1 Pendahuluan

Dokumen ini menganalisis hubungan jangka pendek dan jangka panjang antara curah hujan dan kelembapan udara Kota Bogor terhadap tinggi muka air (TMA) Pintu Air Manggarai menggunakan model ARDL (Autoregressive Distributed Lag). Data yang dipakai adalah data harian Januari 2022 - Desember 2024, lalu hasil modelnya diuji untuk meramal Januari - Maret 2025.

Peubah Keterangan Peran
TMA Tinggi muka air Pintu Air Manggarai (cm) Peubah respon (Y)
PC Curah hujan Kota Bogor (mm) Peubah bebas 1 (X1)
RH Kelembapan relatif Kota Bogor (%) Peubah bebas 2 (X2)

2 Library

library(readxl)     
library(dplyr)      
library(tidyr)    
library(ggplot2)    
library(moments)    
library(corrplot)   
library(forecast)   
library(tseries)    
library(ARDL)       
library(lmtest)     
library(car)        
library(patchwork)  

f <- function(x, d = 2) format(round(as.numeric(x), d), nsmall = d, decimal.mark = ",")

3 Import Data

vars <- c("TMA", "PC", "RH")

dataset <- read_excel("Data_Harian_2022_2024.xlsx") %>%
  mutate(Date = as.Date(Date)) %>%
  select(Date, all_of(vars)) %>%
  as.data.frame()

str(dataset)
## 'data.frame':    1096 obs. of  4 variables:
##  $ Date: Date, format: "2022-01-01" "2022-01-02" ...
##  $ TMA : num  595 588 580 576 574 ...
##  $ PC  : num  0.85 1.79 0.4 0.22 5.7 ...
##  $ RH  : num  86.4 88.1 83.2 80.3 83.2 ...

Interpretasi. Data terdiri dari 1.096 pengamatan harian (1 Januari 2022 sampai 31 Desember 2024) tanpa tanggal yang terlewat, sehingga cocok untuk analisis deret waktu.

4 Eksplorasi Data

4.1 Statistika Deskriptif

ringkasan <- data.frame(
  Min      = sapply(dataset[vars], min),
  Median   = sapply(dataset[vars], median),
  Rata2    = sapply(dataset[vars], mean),
  Maks     = sapply(dataset[vars], max),
  SD       = sapply(dataset[vars], sd),
  Skewness = sapply(dataset[vars], skewness),
  Kurtosis = sapply(dataset[vars], kurtosis)
)
round(ringkasan, 2)
##        Min Median  Rata2   Maks    SD Skewness Kurtosis
## TMA 555.42 620.00 625.35 781.25 39.45     0.53     2.98
## PC    0.00   4.00   6.57 145.03  8.84     5.07    62.20
## RH   67.03  86.29  85.29  94.31  4.71    -0.99     4.03

Interpretasi.

  • TMA rata-ratanya 625,35 cm (rentang 555,42 sampai 781,25 cm) dengan simpangan baku 39,45. Skewness 0,53 dan kurtosis 2,98 (mendekati 3) menunjukkan sebarannya cukup simetris dan mendekati normal.
  • PC (curah hujan) punya skewness 5,07 dan kurtosis 62,20: sangat miring ke kanan dengan ekor sangat berat. Artinya hampir semua hari hujannya kecil (median hanya 4,00 mm), tetapi sesekali terjadi hujan ekstrem (maksimum 145,03 mm).
  • RH (kelembapan) punya skewness -0,99 (miring ke kiri): kelembapan umumnya tinggi (median 86,29%) dan hanya sesekali turun jauh.

Catatan: moments::kurtosis memakai skala di mana distribusi normal bernilai 3.

4.2 Plot Deret Waktu

dataset %>%
  pivot_longer(cols = all_of(vars), names_to = "Peubah", values_to = "Nilai") %>%
  mutate(Peubah = factor(Peubah, levels = vars)) %>%
  ggplot(aes(x = Date, y = Nilai)) +
  geom_line(color = "#1f3b73", linewidth = 0.4) +
  facet_wrap(~ Peubah, scales = "free_y", ncol = 1) +
  labs(title = "Plot Deret Waktu TMA, PC, dan RH", x = "Tanggal", y = NULL) +
  theme_minimal()

Interpretasi. Ketiga peubah tidak membentuk pola musiman yang tegas. TMA naik turun tidak beraturan dengan puncak-puncak tajam (tertinggi pada akhir Februari 2023). Curah hujan menunjukkan lonjakan ekstrem sesekali (puncaknya 2 Juli 2024), sedangkan kelembapan relatif stabil di kisaran tinggi. Lonjakan TMA tampak sering berdekatan dengan lonjakan curah hujan, dugaan awal bahwa keduanya berhubungan.

4.3 Matriks Korelasi

korelasi <- cor(dataset[vars])
round(korelasi, 2)
##      TMA   PC   RH
## TMA 1.00 0.39 0.53
## PC  0.39 1.00 0.42
## RH  0.53 0.42 1.00
corrplot(korelasi, method = "color", addCoef.col = "black", tl.col = "black",
         number.cex = 0.9, mar = c(0, 0, 1, 0))

Interpretasi. TMA berkorelasi positif dengan curah hujan (0,39) dan kelembapan (0,53): semakin tinggi hujan atau kelembapan, TMA cenderung ikut naik. Korelasi antara kedua peubah bebas (PC dan RH) sebesar 0,42, tergolong sedang sehingga belum mengkhawatirkan sebagai multikolinearitas (dicek lagi lewat VIF di bagian uji asumsi).

5 Pembagian Data Latih dan Uji

n_test  <- 160
n_total <- nrow(dataset)

train <- dataset[1:(n_total - n_test), ]
test  <- dataset[(n_total - n_test + 1):n_total, ]

cat("Data train:", nrow(train), "hari (", format(min(train$Date)), "s.d.", format(max(train$Date)), ")\n")
## Data train: 936 hari ( 2022-01-01 s.d. 2024-07-24 )
cat("Data test :", nrow(test),  "hari (", format(min(test$Date)),  "s.d.", format(max(test$Date)),  ")\n")
## Data test : 160 hari ( 2024-07-25 s.d. 2024-12-31 )
bind_rows(mutate(train, Data = "Train"), mutate(test, Data = "Test")) %>%
  ggplot(aes(x = Date, y = TMA, color = Data)) +
  geom_line(linewidth = 0.4) +
  scale_color_manual(values = c(Train = "#00A5A8", Test = "#E8574D")) +
  labs(title = "TMA: Data Train vs Test", x = "Tanggal", y = "TMA (cm)") +
  theme_minimal()

Interpretasi. Data dibagi secara berurutan (bukan acak, karena ini deret waktu): 936 hari (85,4%) untuk melatih model dan 160 hari (14,6%) untuk menguji. Model hanya “melihat” data latih, sedangkan data uji dipakai untuk menilai kemampuan meramal pada data baru.

6 Kestasioneran Data

Model ARDL mensyaratkan tidak ada peubah yang terintegrasi orde 2, I(2). Karena itu perlu dicek kestasioneran tiap peubah pada data latih, secara visual (ACF) dan formal (uji ADF).

6.1 Plot ACF

p1_acf <- ggAcf(train$TMA) + labs(title = "ACF TMA")
p2_acf <- ggAcf(train$PC)  + labs(title = "ACF PC")
p3_acf <- ggAcf(train$RH)  + labs(title = "ACF RH")
p1_acf + p2_acf + p3_acf + plot_layout(ncol = 3)

Interpretasi. ACF TMA dan RH turun sangat lambat (korelasi masih cukup besar sampai lag 30), tanda data belum stasioner atau punya ingatan panjang. ACF PC turun cepat mendekati nol, tanda data cenderung stasioner.

6.2 Uji ADF (Level)

Hipotesis: H0 = data tidak stasioner (ada akar unit), H1 = data stasioner. Tolak H0 jika p-value < 0,05.

uji_adf <- function(x) {
  h <- adf.test(x)
  data.frame(DF = round(unname(h$statistic), 3),
             Lag = unname(h$parameter),
             p_value = round(h$p.value, 3))
}

adf_level <- do.call(rbind, lapply(train[vars], uji_adf))
adf_level$Keputusan <- ifelse(adf_level$p_value < 0.05, "Stasioner", "Tidak stasioner")
adf_level
##         DF Lag p_value       Keputusan
## TMA -4.001   9   0.010       Stasioner
## PC  -7.559   9   0.010       Stasioner
## RH  -3.356   9   0.061 Tidak stasioner

Interpretasi. Pada taraf 5%:

  • TMA (DF = -4,00; p = 0,010): stasioner.
  • PC (DF = -7,56; p = 0,010): stasioner.
  • RH (DF = -3,36; p = 0,061): tidak stasioner.

Jadi peubahnya berupa campuran, ada yang sudah stasioner di level (I(0)) dan ada yang belum (I(1)). Ini tidak masalah untuk ARDL, justru salah satu kelebihan ARDL adalah boleh mencampur peubah I(0) dan I(1), asalkan tidak ada yang I(2). Kita pastikan itu lewat differencing berikut. Namun perlu dicatat bahwa uji bounds pada bagian berikutnya mengasumsikan peubah respon (TMA) bersifat I(1), sedangkan TMA di sini stasioner di level; implikasinya dibahas pada bagian uji bounds.

7 First Difference

Setiap peubah dibeda satu kali (selisih hari ini dengan kemarin), lalu diuji ulang.

d_train <- data.frame(DTMA = diff(train$TMA),
                      DPC  = diff(train$PC),
                      DRH  = diff(train$RH))

adf_diff <- do.call(rbind, lapply(d_train, uji_adf))
adf_diff$Keputusan <- ifelse(adf_diff$p_value < 0.05, "Stasioner", "Tidak stasioner")
adf_diff
##           DF Lag p_value Keputusan
## DTMA -12.518   9    0.01 Stasioner
## DPC  -14.984   9    0.01 Stasioner
## DRH  -13.910   9    0.01 Stasioner
(ggAcf(d_train$DTMA) + labs(title = "ACF DTMA")) +
(ggAcf(d_train$DPC)  + labs(title = "ACF DPC")) +
(ggAcf(d_train$DRH)  + labs(title = "ACF DRH")) + plot_layout(ncol = 3)

Interpretasi. Setelah first difference, ketiga peubah stasioner (p-value < 0,05 dan ACF tidak lagi turun lambat). Berarti tidak ada peubah yang I(2), sehingga syarat penggunaan ARDL terpenuhi. Sebagian peubah stasioner di level (I(0)) dan sebagian stasioner setelah dibeda satu kali (I(1)).

8 Korelasi Silang (CCF)

CCF membantu melihat pada jeda (lag) berapa curah hujan dan kelembapan paling berkaitan dengan perubahan TMA. Lag positif berarti peubah bebas hari sebelumnya berkaitan dengan TMA hari ini.

(ggCcf(d_train$DTMA, d_train$DPC, lag.max = 9) + labs(title = "CCF antara DTMA dan DPC")) +
(ggCcf(d_train$DTMA, d_train$DRH, lag.max = 9) + labs(title = "CCF antara DTMA dan DRH")) +
  plot_layout(ncol = 2)

# nilai CCF pada lag 0 sampai 3 (lag 0 ada di urutan ke-10 dari 19 nilai)
ccf_tabel <- data.frame(
  Lag = 0:3,
  DPC = round(drop(ccf(d_train$DTMA, d_train$DPC, lag.max = 9, plot = FALSE)$acf)[10:13], 3),
  DRH = round(drop(ccf(d_train$DTMA, d_train$DRH, lag.max = 9, plot = FALSE)$acf)[10:13], 3)
)
ccf_tabel
##   Lag    DPC    DRH
## 1   0  0.109  0.094
## 2   1  0.258  0.096
## 3   2 -0.131 -0.064
## 4   3 -0.012  0.046

Interpretasi. Korelasi silang DTMA dengan DPC paling kuat pada lag 1 (0,258), lebih besar daripada lag 0 (0,109). Artinya kenaikan TMA paling berkaitan dengan hujan satu hari sebelumnya, masuk akal karena air hujan dari hulu (Bogor) butuh waktu mengalir ke Pintu Air Manggarai. Korelasi DTMA dengan DRH jauh lebih lemah (lag 1 = 0,096; lag 0 = 0,094). Temuan ini menjadi dasar memberi ruang lag bagi PC dan RH pada pemilihan orde ARDL berikut.

9 Model ARDL

9.1 Pemilihan Orde

Notasi ARDL(q, p_PC, p_RH): TMA dipengaruhi oleh q nilai TMA sebelumnya, serta PC dan RH pada hari yang sama sampai p hari sebelumnya. Orde terbaik dipilih dengan AIC terkecil dari semua kombinasi q = 1 sampai 9, p_PC = 0 sampai 2, dan p_RH = 0 sampai 2 (81 model).

train.ts <- ts(train[, vars])

grid <- expand.grid(q = 1:9, p_PC = 0:2, p_RH = 0:2)
grid$AIC <- NA
for (i in seq_len(nrow(grid))) {
  m <- ardl(TMA ~ PC + RH, data = train.ts,
            order = c(grid$q[i], grid$p_PC[i], grid$p_RH[i]))
  grid$AIC[i] <- AIC(m)
}
grid <- grid[order(grid$AIC), ]
head(grid, 5)
##    q p_PC p_RH      AIC
## 81 9    2    2 8313.011
## 27 9    2    0 8314.567
## 54 9    2    1 8315.290
## 72 9    1    2 8319.151
## 80 8    2    2 8322.825
orde_opt <- as.numeric(grid[1, c("q", "p_PC", "p_RH")])
ord_q <- orde_opt[1]; p_pc <- orde_opt[2]; p_rh <- orde_opt[3]

Interpretasi. Model terbaik adalah ARDL(9, 2, 2). Artinya TMA hari ini dipengaruhi oleh TMA 9 hari terakhir (air sungai punya “ingatan” yang panjang), curah hujan hari ini sampai 2 hari sebelumnya, dan kelembapan hari ini sampai 2 hari sebelumnya. Orde PC sampai lag 2 sejalan dengan temuan CCF bahwa efek hujan terasa dengan jeda.

9.2 Model ARDL (Level)

model_ardl <- ardl(TMA ~ PC + RH, data = train.ts, order = orde_opt)
summary(model_ardl)
## 
## Time series regression with "ts" data:
## Start = 10, End = 936
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -125.798  -11.573   -2.591    6.303  134.954 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 73.37831   16.56558   4.430 1.06e-05 ***
## L(TMA, 1)    0.39650    0.03285  12.071  < 2e-16 ***
## L(TMA, 2)    0.15393    0.03420   4.500 7.66e-06 ***
## L(TMA, 3)    0.07621    0.03389   2.249  0.02477 *  
## L(TMA, 4)    0.07606    0.03364   2.261  0.02398 *  
## L(TMA, 5)   -0.02272    0.03360  -0.676  0.49906    
## L(TMA, 6)    0.02263    0.03345   0.676  0.49894    
## L(TMA, 7)    0.02280    0.03336   0.683  0.49453    
## L(TMA, 8)    0.05384    0.03310   1.626  0.10421    
## L(TMA, 9)    0.05468    0.02966   1.843  0.06560 .  
## PC           0.65993    0.09618   6.861 1.26e-11 ***
## L(PC, 1)     0.85793    0.10367   8.276 4.50e-16 ***
## L(PC, 2)    -0.29453    0.10391  -2.835  0.00469 ** 
## RH           0.84351    0.29902   2.821  0.00489 ** 
## L(RH, 1)     0.03773    0.34337   0.110  0.91252    
## L(RH, 2)    -0.61716    0.30060  -2.053  0.04035 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.23 on 911 degrees of freedom
## Multiple R-squared:  0.7369, Adjusted R-squared:  0.7325 
## F-statistic: 170.1 on 15 and 911 DF,  p-value: < 2.2e-16

Interpretasi. Model ARDL(9, 2, 2) mampu menjelaskan sekitar 0,737 (R-squared) keragaman TMA, tetapi sebagian besar berasal dari lag TMA sendiri (TMA sangat persisten), sehingga angka ini tidak boleh dibaca sebagai kekuatan PC dan RH dalam menjelaskan TMA. Koefisien pada model level ini sulit dibaca satu per satu karena ada banyak suku lag yang saling berkaitan, sehingga makna jangka pendek dan jangka panjangnya kita lihat lewat bentuk error correction dan multiplier di bagian berikut.

9.3 Uji Hubungan Jangka Panjang (Bounds Test)

Uji ini menjawab: apakah TMA, PC, dan RH punya hubungan jangka panjang? Model ARDL ditulis ulang dalam bentuk Unrestricted Error Correction Model (UECM), lalu dilakukan uji F Pesaran-Shin-Smith (case 2: intercept terbatas, tanpa tren).

  • H0: tidak ada hubungan jangka panjang
  • H1: ada hubungan jangka panjang
  • Jika F hitung di atas batas atas I(1) pada nilai kritis, maka tolak H0.
model_uecm <- uecm(model_ardl)

uji_bound <- bounds_f_test(model_uecm, case = 2)
uji_bound
## 
##  Bounds F-test (Wald) for no cointegration
## 
## data:  d(TMA) ~ L(TMA, 1) + L(PC, 1) + L(RH, 1) + d(L(TMA, 1)) + d(L(TMA,     2)) + d(L(TMA, 3)) + d(L(TMA, 4)) + d(L(TMA, 5)) + d(L(TMA,     6)) + d(L(TMA, 7)) + d(L(TMA, 8)) + d(PC) + d(L(PC, 1)) +     d(RH) + d(L(RH, 1))
## F = 30.489, p-value = 1e-06
## alternative hypothesis: Possible cointegration
## null values:
##    k    T 
##    2 1000
F_bound <- as.numeric(uji_bound$statistic)
p_bound <- if (is.null(uji_bound$p.value)) NA else max(uji_bound$p.value)
kes_bound <- ifelse(is.na(p_bound), "(lihat nilai kritis pada output)",
             ifelse(p_bound < 0.05, "terdapat hubungan jangka panjang yang stabil",
                                    "tidak terdapat hubungan jangka panjang"))

Interpretasi. Nilai F hitung = 30,49, jauh di atas batas atas nilai kritis I(1) bahkan pada taraf 1%. Kesimpulannya terdapat hubungan jangka panjang yang stabil: artinya TMA, curah hujan, dan kelembapan memiliki hubungan keseimbangan yang stabil dalam jangka panjang. Catatan: uji bounds mengasumsikan peubah respon I(1), sedangkan TMA terdeteksi stasioner di level (uji ADF). Karena itu penolakan H0 tidak dapat dibaca sebagai kointegrasi dalam arti ketat, melainkan sebagai hubungan jangka panjang yang stabil (mean-reverting). Hubungan ini tetap dapat ditafsirkan lewat koefisien penyesuaian dan pengali jangka panjang pada bagian berikut.

10 Hubungan Jangka Pendek dan Jangka Panjang

10.1 Jangka Pendek (Error Correction Model)

Dalam bentuk UECM, d(...) berarti perubahan harian dan L(TMA, 1) adalah TMA satu hari sebelumnya. Koefisien L(TMA, 1) adalah kecepatan penyesuaian (error correction term).

summary(model_uecm)
## 
## Time series regression with "ts" data:
## Start = 10, End = 936
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -125.798  -11.573   -2.591    6.303  134.954 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  73.37831   16.56558   4.430 1.06e-05 ***
## L(TMA, 1)    -0.16608    0.02485  -6.684 4.06e-11 ***
## L(PC, 1)      1.22333    0.13047   9.376  < 2e-16 ***
## L(RH, 1)      0.26408    0.20974   1.259 0.208325    
## d(L(TMA, 1)) -0.43742    0.03572 -12.246  < 2e-16 ***
## d(L(TMA, 2)) -0.28350    0.03750  -7.559 9.86e-14 ***
## d(L(TMA, 3)) -0.20729    0.03690  -5.618 2.57e-08 ***
## d(L(TMA, 4)) -0.13123    0.03646  -3.600 0.000336 ***
## d(L(TMA, 5)) -0.15395    0.03539  -4.350 1.51e-05 ***
## d(L(TMA, 6)) -0.13132    0.03462  -3.793 0.000159 ***
## d(L(TMA, 7)) -0.10852    0.03304  -3.284 0.001061 ** 
## d(L(TMA, 8)) -0.05468    0.02966  -1.843 0.065601 .  
## d(PC)         0.65993    0.09618   6.861 1.26e-11 ***
## d(L(PC, 1))   0.29453    0.10391   2.835 0.004691 ** 
## d(RH)         0.84351    0.29902   2.821 0.004893 ** 
## d(L(RH, 1))   0.61716    0.30060   2.053 0.040348 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.23 on 911 degrees of freedom
## Multiple R-squared:  0.3726, Adjusted R-squared:  0.3623 
## F-statistic: 36.07 on 15 and 911 DF,  p-value: < 2.2e-16
sr <- as.data.frame(summary(model_uecm)$coefficients)
colnames(sr) <- c("b", "se", "t", "p")
ect <- unlist(sr["L(TMA, 1)", c("b", "t", "p")])

Interpretasi jangka pendek.

  • Kecepatan penyesuaian (ECT). Koefisien L(TMA, 1) = -0,166 (t = -6,68; p < 0,001). Tandanya negatif dan signifikan, sesuai syarat penyesuaian menuju keseimbangan jangka panjang. Artinya sekitar 16,6% penyimpangan TMA dari keseimbangan jangka panjangnya dikoreksi setiap hari, atau separuh penyimpangan hilang dalam kira-kira 3,8 hari.
  • Curah hujan (PC). Koefisien d(PC) = 0,66 (p < 0,001): kenaikan curah hujan 1 mm pada hari yang sama diikuti kenaikan TMA sekitar 0,66 cm pada hari itu juga. Hujan hari sebelumnya, d(L(PC, 1)) = 0,29 (p = 0,005), masih menambah pengaruh positif.
  • Kelembapan (RH). Koefisien d(RH) = 0,84 (p = 0,005) dan d(L(RH, 1)) = 0,62 (p = 0,040): kenaikan kelembapan diikuti kenaikan TMA dalam jangka pendek.

Semua pengaruh di atas dibaca ceteris paribus (peubah lain dianggap tetap). Suku d(L(TMA, k)) cenderung bertanda negatif karena bentuk penulisan ulang (reparameterisasi) model ARDL ke UECM, bukan karena TMA berbalik arah setelah lonjakan. Dinamika sebenarnya dilihat dari koefisien lag TMA yang positif pada model level, yaitu TMA sangat persisten namun tetap kembali ke keseimbangan.

10.2 Jangka Panjang

Koefisien jangka panjang (long-run multiplier) adalah pengaruh total peubah bebas terhadap TMA setelah semua penyesuaian selesai, dihitung dari model ARDL dengan delta method.

lr <- multipliers(model_ardl)
lr
##          Term   Estimate Std. Error  t value     Pr(>|t|)
## 1 (Intercept) 441.822663  95.541683 4.624397 4.299790e-06
## 2          PC   7.365892   1.271358 5.793718 9.477868e-09
## 3          RH   1.590065   1.162977 1.367237 1.718884e-01
ambil_lr <- function(nama, kolom) lr[lr$Term == nama, kolom]

Interpretasi jangka panjang.

  • Curah hujan (PC): koefisien jangka panjang 7,37 (p < 0,001), signifikan. Jika curah hujan naik 1 mm dan bertahan, TMA akan naik sekitar 7,37 cm dalam jangka panjang, jauh lebih besar dibanding efek jangka pendeknya karena pengaruh hujan terakumulasi lewat beberapa hari.
  • Kelembapan (RH): koefisien jangka panjang 1,59 dengan p = 0,172, sehingga tidak signifikan pada taraf 5%. Kelembapan hanya berpengaruh nyata pada jangka pendek, tetapi tidak terbukti berpengaruh jangka panjang setelah curah hujan dikendalikan.

Jadi, curah hujan adalah penggerak utama TMA baik jangka pendek maupun panjang, sedangkan kelembapan berperan sebagai pelengkap jangka pendek.

11 Uji Asumsi Residual

Untuk uji asumsi dan peramalan, model yang sama (ARDL dalam bentuk level) dimodelkan ulang dengan lm() agar mudah dipakai memprediksi. Hasilnya identik dengan ardl().

buat_lag <- function(df, ord_q, p_pc, p_rh) {
  hasil <- data.frame(TMA = df$TMA, PC = df$PC, RH = df$RH)
  for (i in 1:ord_q)       hasil[[paste0("TMA_L", i)]] <- dplyr::lag(df$TMA, i)
  for (i in seq_len(p_pc)) hasil[[paste0("PC_L", i)]]  <- dplyr::lag(df$PC, i)
  for (i in seq_len(p_rh)) hasil[[paste0("RH_L", i)]]  <- dplyr::lag(df$RH, i)
  hasil
}

lag_train <- na.omit(buat_lag(train, ord_q, p_pc, p_rh))
model_lm  <- lm(TMA ~ ., data = lag_train)

# cek: koefisien sama dengan model ardl()
all.equal(sort(unname(coef(model_lm))), sort(unname(coef(model_ardl))))
## [1] TRUE

11.1 Autokorelasi Residual (Ljung-Box)

H0: tidak ada autokorelasi residual.

res <- residuals(model_lm)
uji_lb <- Box.test(res, lag = 20, type = "Ljung-Box")
uji_lb
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 8.7487, df = 20, p-value = 0.9856

Interpretasi. p-value = 0,986 (> 0,05), maka H0 tidak ditolak: tidak ada autokorelasi pada residual. Ini hasil yang baik, artinya struktur lag pada ARDL sudah cukup menangkap ketergantungan antarwaktu.

11.2 Normalitas Residual (Jarque-Bera)

H0: residual berdistribusi normal.

uji_jb <- jarque.bera.test(res)
uji_jb
## 
##  Jarque Bera Test
## 
## data:  res
## X-squared = 1574.2, df = 2, p-value < 2.2e-16

Interpretasi. p-value < 0,05 (statistik JB = 1574,2), maka H0 ditolak: residual tidak berdistribusi normal. Penyebabnya wajar, yaitu adanya lonjakan TMA dan hujan ekstrem (banjir) yang menghasilkan galat besar dan jarang terjadi.

11.3 Heteroskedastisitas (Breusch-Pagan)

H0: ragam residual konstan (homoskedastis).

uji_bp <- bptest(model_lm)
uji_bp
## 
##  studentized Breusch-Pagan test
## 
## data:  model_lm
## BP = 189.72, df = 15, p-value < 2.2e-16

Interpretasi. p-value < 0,05 (BP = 189,7; db = 15), maka H0 ditolak: ragam residual tidak konstan (heteroskedastis), galat cenderung lebih besar saat TMA sedang tinggi (korelasi antara galat mutlak dan nilai prediksi sekitar 0,3).

11.4 Multikolinearitas (VIF)

vif(model_lm)
##       PC       RH   TMA_L1   TMA_L2   TMA_L3   TMA_L4   TMA_L5   TMA_L6 
## 1.446758 3.993935 3.741219 4.058723 3.990133 3.934067 3.929661 3.900801 
##   TMA_L7   TMA_L8   TMA_L9    PC_L1    PC_L2    RH_L1    RH_L2 
## 3.881715 3.825161 3.071692 1.679632 1.686699 5.267513 4.035601

Interpretasi. Nilai VIF terbesar adalah 5,27, masih jauh di bawah batas umum 10. Jadi tidak ada multikolinearitas yang serius, walaupun beberapa suku lag wajar saling berkorelasi.

Ringkasan uji asumsi. Residual bebas autokorelasi dan tidak ada multikolinearitas serius, tetapi tidak normal dan heteroskedastis. Dampaknya, koefisien tetap wajar dipakai untuk menggambarkan hubungan dan meramal, tetapi galat baku dan p-value sebaiknya dibaca hati-hati (kurang andal), terutama untuk koefisien yang p-valuenya mepet 0,05 seperti d(L(RH, 1)). Sebagai pengembangan, dapat dipertimbangkan galat baku yang tahan heteroskedastisitas (HAC) atau transformasi data.

12 Peramalan dan Evaluasi

Peramalan dilakukan satu langkah ke depan (one-step-ahead): prediksi TMA hari ini memakai TMA hari-hari sebelumnya yang sudah diamati, ditambah curah hujan dan kelembapan pada hari-hari yang relevan. Koefisien tetap dari model yang dilatih pada data latih. Dua periode diuji: data uji (25 Juli - 31 Desember 2024) dan data baru Januari - Maret 2025.

data_2025 <- read_excel("Data_Jan25_Mar25.xlsx") %>%
  mutate(Date = as.Date(Date)) %>%
  select(Date, all_of(vars)) %>%
  as.data.frame()

# gabungkan agar lag di awal 2025 diambil dari data akhir 2024
data_all <- bind_rows(dataset, data_2025)
urutan   <- seq_len(nrow(data_all))

data_all$Periode <- case_when(urutan <= nrow(train)   ~ "Train",
                              urutan <= nrow(dataset) ~ "Test",
                              TRUE                    ~ "2025")
data_all$Prediksi <- predict(model_lm, newdata = buat_lag(data_all, ord_q, p_pc, p_rh))
data_all$Naive    <- dplyr::lag(data_all$TMA, 1)   # pembanding: TMA kemarin

table(data_all$Periode)
## 
##  2025  Test Train 
##    90   160   936

12.1 Ukuran Akurasi

Ukuran akurasi: MAE (rata-rata galat mutlak, cm), RMSE (akar rata-rata galat kuadrat, cm; lebih peka terhadap galat besar), dan MAPE (rata-rata galat persentase, %). Sebagai pembanding ditambahkan metode naif, yaitu menebak TMA hari ini sama dengan TMA kemarin.

evaluasi <- function(aktual, prediksi) {
  c(MAE  = mean(abs(aktual - prediksi)),
    RMSE = sqrt(mean((aktual - prediksi)^2)),
    MAPE = mean(abs((aktual - prediksi) / aktual)) * 100)
}

tabel_eval <- bind_rows(lapply(c("Train", "Test", "2025"), function(p) {
  d <- data_all[data_all$Periode == p & !is.na(data_all$Prediksi), ]
  data.frame(Periode = p, Metode = c("ARDL", "Naive"),
             rbind(evaluasi(d$TMA, d$Prediksi), evaluasi(d$TMA, d$Naive)))
}))
round_df <- tabel_eval; round_df[, 3:5] <- round(round_df[, 3:5], 2)
round_df
##   Periode Metode   MAE  RMSE MAPE
## 1   Train   ARDL 14.27 21.04 2.21
## 2   Train  Naive 16.77 26.57 2.59
## 3    Test   ARDL  9.88 16.45 1.56
## 4    Test  Naive  9.61 18.77 1.49
## 5    2025   ARDL 11.14 16.80 1.67
## 6    2025  Naive 11.35 22.64 1.66
ev <- function(p, metode, kol) tabel_eval[tabel_eval$Periode == p & tabel_eval$Metode == metode, kol]

Interpretasi.

  • Data uji (Jul - Des 2024): MAE = 9,88 cm, RMSE = 16,45 cm, MAPE = 1,56%.
  • Data baru (Jan - Mar 2025): MAE = 11,14 cm, RMSE = 16,80 cm, MAPE = 1,67%.

Galat pada data baru (2025) setara dengan data uji, sehingga model tampak stabil dan tidak overfitting. Namun MAPE yang kecil ini tidak otomatis menandakan model yang unggul, karena metode naif pun mencapai MAPE yang hampir sama (lihat perbandingan di bawah): TMA sangat persisten sehingga menebak “sama seperti kemarin” saja sudah akurat.

Perbandingan dengan metode naif. Pada MAE dan MAPE, model ARDL hanya setara dengan tebakan “TMA hari ini = kemarin” (MAPE naif data uji 1,49%, data 2025 1,66%), karena TMA sangat bergantung pada nilai kemarin. Keunggulan ARDL tampak pada RMSE (data uji 16,45 vs 18,77 cm; 2025 16,80 vs 22,64 cm): ARDL lebih baik menangkap lonjakan TMA, kemungkinan karena memanfaatkan informasi curah hujan dan kelembapan yang tidak dimiliki metode naif.

12.2 Plot Aktual vs Prediksi

plot_aktual_prediksi <- function(periode, judul) {
  data_all %>%
    filter(Periode == periode) %>%
    ggplot(aes(x = Date)) +
    geom_line(aes(y = TMA, color = "Aktual"), linewidth = 0.5) +
    geom_line(aes(y = Prediksi, color = "Prediksi"), linewidth = 0.5) +
    scale_color_manual(values = c(Aktual = "black", Prediksi = "#2563eb")) +
    labs(title = judul, x = "Tanggal", y = "TMA (cm)", color = NULL) +
    theme_minimal() + theme(legend.position = "top")
}

plot_aktual_prediksi("Test", "Data Uji: Aktual vs Prediksi (25 Jul - 31 Des 2024)")

plot_aktual_prediksi("2025", "Data Baru: Aktual vs Prediksi (Jan - Mar 2025)")

Interpretasi. Garis prediksi (biru) mengikuti pola garis aktual (hitam) dengan baik, termasuk naik turun harian. Kelemahan yang terlihat: puncak lonjakan ekstrem cenderung diprediksi lebih rendah dari kenyataan (underestimate), wajar karena kejadian banjir besar jarang muncul di data latih. Perlu diingat juga bahwa ini peramalan satu hari ke depan yang memakai curah hujan dan kelembapan aktual; untuk meramal beberapa hari ke depan, nilai kedua peubah bebas itu sendiri harus diramal terlebih dahulu.

13 Kesimpulan

  1. Kestasioneran. TMA dan PC stasioner di level, sedangkan RH stasioner setelah first difference. Tidak ada peubah I(2), sehingga ARDL layak dipakai.
  2. Model terbaik adalah ARDL(9, 2, 2) berdasarkan AIC terkecil.
  3. Hubungan jangka panjang. Bounds test (F = 30,49) menunjukkan terdapat hubungan jangka panjang yang stabil antara TMA, curah hujan, dan kelembapan. Karena TMA stasioner di level, ini bukan kointegrasi dalam arti ketat.
  4. Jangka pendek. Koefisien penyesuaian (ECT) = -0,166: sekitar 16,6% penyimpangan dari keseimbangan dikoreksi tiap hari. Curah hujan (0,66 cm per mm pada hari yang sama) dan kelembapan sama-sama berpengaruh positif.
  5. Jangka panjang. Hanya curah hujan yang berpengaruh signifikan (7,37 cm per mm); kelembapan tidak signifikan.
  6. Asumsi. Residual bebas autokorelasi dan tanpa multikolinearitas serius, tetapi tidak normal dan heteroskedastis, sehingga p-value dibaca hati-hati.
  7. Akurasi. MAPE data uji 1,56% dan data 2025 1,67% namun pada MAE dan MAPE hanya setara dengan metode naif; keunggulan ARDL tampak pada RMSE.