library(dLagM)
## Warning: package 'dLagM' was built under R version 4.5.3
## Loading required package: nardl
## Warning: package 'nardl' was built under R version 4.5.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Loading required package: dynlm
## Warning: package 'dynlm' was built under R version 4.5.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.5.2
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(dynlm)
library(MLmetrics)
## Warning: package 'MLmetrics' was built under R version 4.5.3
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:dLagM':
## 
##     MAPE
## The following object is masked from 'package:base':
## 
##     Recall
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.5.2
library(car)
## Warning: package 'car' was built under R version 4.5.2
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.5.2
library(readxl)

Data

Analisis memakai Tinggi Muka Air (TMA) sebagai variabel Y dan Curah Hujan (PC) sebagai variabel X. Data harian 2022–2024 berjumlah 1.096 observasi, dibagi berurutan menjadi 876 data training (80%) dan 220 data testing (20%).

# Impor data
raw <- as.data.frame(read_excel("Data_Harian_2022_2024.xlsx"))
raw <- raw[order(raw$Date), ]

# Bentuk data dengan nama Yt dan Xt (sesuai panduan)
data <- data.frame(Yt = raw$TMA, Xt = raw$PC)

# Membagi data: 80% training, 20% testing (urut waktu)
n_train <- floor(0.8 * nrow(data))
train <- data[1:n_train, ]
test  <- data[(n_train + 1):nrow(data), ]
h     <- nrow(test)

# Mengubah format menjadi objek Time Series (ts)
train.ts <- ts(train)
test.ts  <- ts(test)
data.ts  <- ts(data)

Buat Plot Data

# Plot keseluruhan data (tanpa pembeda train/test)
tgl <- as.Date(raw$Date)

par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))

# Y: TMA
plot(tgl, data$Yt, type = "l", col = "blue",
     xlab = "", ylab = "TMA", main = "Y: TMA (Tinggi Muka Air)")
grid()

# X: PC
plot(tgl, data$Xt, type = "l", col = "darkgreen",
     xlab = "Tanggal", ylab = "PC", main = "X: PC (Curah Hujan)")
grid()

par(mfrow = c(1, 1))

Kedua variabel berfluktuasi harian. TMA memiliki banyak lonjakan tajam dengan kisaran sekitar 550–780. Curah hujan sebagian besar bernilai rendah dengan sesekali lonjakan tinggi. Pada sekitar Agustus–November 2023 (musim kemarau) curah hujan mendekati nol dan TMA bertahan datar di level terendah. Pola ini menunjukkan TMA bergerak searah dengan curah hujan, sehingga X layak dipakai untuk menjelaskan Y.

Autoregressive Distributed Lag (ARDL)

Notasi pada ardlDlm: p = panjang lag X, q = panjang lag Y.

Pemodelan (p = 1, q = 1)

# Model ARDL dengan lag X = 1 dan lag Y = 1
model.ardl <- ardlDlm(formula = Yt ~ Xt, data = train, p = 1, q = 1)
summary(model.ardl)
## 
## Time series regression with "ts" data:
## Start = 2, End = 876
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -85.126 -10.938  -3.332  10.384 119.054 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 188.8437    12.6415  14.938  < 2e-16 ***
## Xt.t          0.9565     0.1227   7.795 1.84e-14 ***
## Xt.1          0.9620     0.1309   7.347 4.67e-13 ***
## Yt.1          0.6789     0.0207  32.790  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.88 on 871 degrees of freedom
## Multiple R-squared:  0.7041, Adjusted R-squared:  0.7031 
## F-statistic:   691 on 3 and 871 DF,  p-value: < 2.2e-16

Ketiga peubah berpengaruh signifikan terhadap TMA (semua \(p-value < 0,05\)):

  • Curah hujan hari ini (Xt) meningkatkan TMA sebesar \(0,9565\) satuan untuk setiap kenaikan 1 satuan curah hujan.

  • Curah hujan kemarin (Xt-1) menambah pengaruh \(0,9620\), sehingga hujan punya efek tertunda satu hari.

  • TMA kemarin (Yt-1) berpengaruh \(0,6789\), yang menunjukkan TMA sangat bergantung pada kondisinya sendiri di hari sebelumnya.

Persamaannya: \(Ŷt = 188,8437 + 0,9565 Xt + 0,9620 Xt-1 + 0,6789 Yt-1\). Nilai \(R²\) sebesar \(0,7041\) berarti sekitar 70,4% keragaman TMA dapat dijelaskan model. Uji F signifikan (p < 2,2e-16), sehingga model layak secara simultan. Efek jangka panjang curah hujan sebesar \((0,9565 + 0,9620) / (1 − 0,6789) ≈ 5,98\), artinya kenaikan 1 satuan curah hujan yang bertahan menaikkan TMA sekitar 5,98 satuan pada akhirnya.

Lag Optimum ARDL

# Mencari kombinasi p dan q terbaik berdasarkan AIC
model.ardl.opt <- ardlBoundOrders(data = data.frame(train), ic = "AIC",
                                  formula = Yt ~ Xt, max.p = 7, max.q = 7)

# Menarik kombinasi lag dengan AIC terkecil
min_p <- c()
for(i in 1:length(model.ardl.opt$Stat.table)){
  min_p[i] <- min(model.ardl.opt$Stat.table[[i]], na.rm = TRUE)
}
q_opt <- which(min_p == min(min_p, na.rm = TRUE))
p_opt <- which(model.ardl.opt$Stat.table[[q_opt]] ==
                 min(model.ardl.opt$Stat.table[[q_opt]], na.rm = TRUE))

data.frame("Lag_Y_opt(q)" = q_opt, "Lag_X_opt(p)" = p_opt,
           "AIC" = model.ardl.opt$min.Stat)
##   Lag_Y_opt.q. Lag_X_opt.p.      AIC
## 1            7            1 7731.751

Berdasarkan AIC terkecil (7731,751), kombinasi terbaik adalah q = 7 dan p = 1, yaitu TMA dipengaruhi oleh tujuh hari TMA sebelumnya dan curah hujan hari ini serta kemarin. Perhatikan bahwa q optimum jatuh tepat di batas pencarian (max.q = 7), sehingga lag optimum sebenarnya mungkin lebih tinggi. Hal ini bisa dicek dengan menaikkan max.q.

Membangun ulang ARDL dengan lag optimum

model.ardl.best <- ardlDlm(formula = Yt ~ Xt, data = train, p = p_opt, q = q_opt)
summary(model.ardl.best)
## 
## Time series regression with "ts" data:
## Start = 8, End = 876
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -57.799 -11.686  -3.027   7.113 127.256 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 120.57304   12.87234   9.367  < 2e-16 ***
## Xt.t          0.83802    0.11347   7.385  3.6e-13 ***
## Xt.1          1.26147    0.12258  10.291  < 2e-16 ***
## Yt.1          0.35048    0.03221  10.880  < 2e-16 ***
## Yt.2          0.13037    0.03383   3.853 0.000125 ***
## Yt.3          0.09836    0.03388   2.904 0.003784 ** 
## Yt.4          0.08064    0.03403   2.370 0.018031 *  
## Yt.5          0.00119    0.03383   0.035 0.971945    
## Yt.6          0.04512    0.03361   1.343 0.179787    
## Yt.7          0.08004    0.02976   2.690 0.007284 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.94 on 859 degrees of freedom
## Multiple R-squared:  0.7533, Adjusted R-squared:  0.7507 
## F-statistic: 291.5 on 9 and 859 DF,  p-value: < 2.2e-16

Model dengan lag optimum lebih baik pada data training daripada ARDL (1,1): R² naik dari 0,7041 menjadi 0,7533 dan standar error sisaan turun dari 22,88 menjadi 20,94. Semua peubah X signifikan (Xt = 0,838 dan Xt-1 = 1,261). Lag Y yang signifikan adalah Yt-1 sampai Yt-4 dan Yt-7, sedangkan Yt-5 (p = 0,972) dan Yt-6 (p = 0,180) tidak signifikan. Pengaruh TMA masa lalu menurun seiring jauhnya lag (0,350 pada lag 1 menjadi 0,130 pada lag 2), meski muncul kembali kecil pada lag 7, yang mungkin mencerminkan pola mingguan.

Peramalan

# Peramalan ARDL awal (p=1, q=1)
fore.ardl <- forecast(model = model.ardl, x = test$Xt, h = h)
mape.ardl <- MAPE(fore.ardl$forecasts, test$Yt)
cat("MAPE Model ARDL (1,1):", mape.ardl, "\n")
## MAPE Model ARDL (1,1): 0.03320921
# Peramalan ARDL optimum
fore.ardl.best <- forecast(model = model.ardl.best, x = test$Xt, h = h)
mape.ardl.best <- MAPE(fore.ardl.best$forecasts, test$Yt)
cat("MAPE Model ARDL Optimum (p=", p_opt, ", q=", q_opt, "):", mape.ardl.best, "\n")
## MAPE Model ARDL Optimum (p= 1 , q= 7 ): 0.04029873

Nilai MAPE ARDL (1,1) sebesar 0,0332 (3,32%) dan ARDL optimum (1,7) sebesar 0,0403 (4,03%). Kedua model sangat akurat karena error di bawah 10%. Keakuratan yang tinggi ini juga dipengaruhi level TMA yang besar (550–780), sehingga error relatif menjadi kecil. Model optimum lebih baik pada data training tetapi sedikit lebih buruk pada data testing. Ini mungkin terjadi karena data testing memuat curah hujan ekstrem yang tidak ada pada data training.

Pendekatan dengan Library dynlm dan Uji Asumsi

# Sama dengan DLM (q=1)
cons_lm1 <- dynlm(Yt ~ Xt + L(Xt), data = train.ts)

# Setara dengan Model Autoregressive Murni (hanya Xt dan Yt-1)
cons_lm2 <- dynlm(Yt ~ Xt + L(Yt), data = train.ts)

# Sama dengan ARDL (p=1, q=1)
cons_lm3 <- dynlm(Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)

# Sama dengan ARDL optimum (p_opt, q_opt)
cons_best <- dynlm(Yt ~ Xt + L(Xt, 1:p_opt) + L(Yt, 1:q_opt), data = train.ts)
summary(cons_best)
## 
## Time series regression with "ts" data:
## Start = 8, End = 876
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt, 1:p_opt) + L(Yt, 1:q_opt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -57.799 -11.686  -3.027   7.113 127.256 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     120.57304   12.87234   9.367  < 2e-16 ***
## Xt                0.83802    0.11347   7.385  3.6e-13 ***
## L(Xt, 1:p_opt)    1.26147    0.12258  10.291  < 2e-16 ***
## L(Yt, 1:q_opt)1   0.35048    0.03221  10.880  < 2e-16 ***
## L(Yt, 1:q_opt)2   0.13037    0.03383   3.853 0.000125 ***
## L(Yt, 1:q_opt)3   0.09836    0.03388   2.904 0.003784 ** 
## L(Yt, 1:q_opt)4   0.08064    0.03403   2.370 0.018031 *  
## L(Yt, 1:q_opt)5   0.00119    0.03383   0.035 0.971945    
## L(Yt, 1:q_opt)6   0.04512    0.03361   1.343 0.179787    
## L(Yt, 1:q_opt)7   0.08004    0.02976   2.690 0.007284 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.94 on 859 degrees of freedom
## Multiple R-squared:  0.7533, Adjusted R-squared:  0.7507 
## F-statistic: 291.5 on 9 and 859 DF,  p-value: < 2.2e-16

Model cons_best menghasilkan koefisien, \(R²\) (0,7533), dan standar error sisaan (20,94) yang identik dengan ardlDlm (1,7). Ini memastikan kedua pendekatan konsisten, dan cons_best dipakai untuk uji asumsi karena bisa langsung dipanggil di fungsi uji sisaan.

Uji asumsi model (pada model optimum cons_best)

  • Uji autokorelasi (H0: tidak ada autokorelasi)
dwtest(cons_best)
## 
##  Durbin-Watson test
## 
## data:  cons_best
## DW = 1.9189, p-value = 0.1115
## alternative hypothesis: true autocorrelation is greater than 0

DW = 1,9189, \(p-value = 0,1115 > 0,05\), sehingga H0 tidak ditolak (tidak ada autokorelasi).

  • Uji heteroskedastisitas (H0: ragam sisaan homogen)
bptest(cons_best)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_best
## BP = 77.44, df = 9, p-value = 5.204e-13

Nilai BP = 77,44 dengan \(p-value = 5,204e-13 < 0,05\), sehingga H0 ditolak. Ragam sisaan tidak homogen (terdapat heteroskedastisitas), sehingga asumsi homoskedastisitas tidak terpenuhi.

  • Uji normalitas (H0: sisaan berdistribusi normal)
shapiro.test(residuals(cons_best))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_best)
## W = 0.90493, p-value < 2.2e-16

Nilai W = 0,905 dengan \(p-value < 2,2e-16 < 0,05\), sehingga H0 ditolak. Sisaan tidak berdistribusi normal.

  • Uji multikolinearitas (VIF > 10 berarti bermasalah)
vif(cons_best)
##                    GVIF Df GVIF^(1/(2*Df))
## Xt             1.345213  1        1.159833
## L(Xt, 1:p_opt) 1.571457  1        1.253578
## L(Yt, 1:q_opt) 1.318218  7        1.019930

Semua nilai GVIF berada di bawah 10. Tidak ada multikolinearitas antar peubah bebas, sehingga asumsi ini terpenuhi.

Perbandingan Performa Model

akurasi <- matrix(c(mape.ardl, mape.ardl.best))
row.names(akurasi) <- c("ARDL (1,1)", paste0("ARDL Optimum (", p_opt, ",", q_opt, ")"))
colnames(akurasi) <- "MAPE"
print(akurasi)
##                          MAPE
## ARDL (1,1)         0.03320921
## ARDL Optimum (1,7) 0.04029873
# Visualisasi aktual vs peramalan
plot(test$Yt, type = "b", col = "black", lwd = 2,
     ylim = range(c(test$Yt, fore.ardl$forecasts, fore.ardl.best$forecasts)),
     ylab = "TMA", xlab = "Periode testing (hari)",
     main = "Perbandingan Model Peramalan ARDL")
lines(fore.ardl$forecasts,      col = "green", type = "b", pch = 18)
lines(fore.ardl.best$forecasts, col = "red",   type = "b", pch = 16)
legend("topleft",
       legend = c("Aktual", "ARDL (1,1)", paste0("ARDL Optimum (", p_opt, ",", q_opt, ")")),
       col = c("black", "green", "red"), lty = 1, pch = c(1, 18, 16),
       cex = 0.8, inset = 0.02)

  • Tabel MAPE menunjukkan ARDL (1,1) (0,0332) lebih kecil daripada ARDL optimum (1,7) (0,0403). Untuk tujuan peramalan pada data ini, model yang lebih sederhana bekerja sedikit lebih baik. Sebaliknya, ARDL (1,7) lebih baik dari sisi kecocokan pada data training (AIC, \(R²\), dan standar error sisaan).

  • Garis peramalan kedua model mengikuti pola umum data aktual. Sesuaikan kalimat ini dengan plot Anda, terutama apakah lonjakan tajam TMA ikut tertangkap atau justru terpangkas.

Kesimpulan

Curah hujan hari ini dan kemarin berpengaruh signifikan terhadap TMA, dan TMA juga dipengaruhi nilainya di hari-hari sebelumnya. Kedua model memberi akurasi peramalan yang baik (MAPE 3–4%), dengan ARDL (1,1) sedikit lebih unggul pada data testing. Namun uji asumsi menunjukkan sisaan mengandung heteroskedastisitas dan tidak normal, sehingga uji signifikansi (nilai-p) sebaiknya ditafsirkan dengan hati-hati.