Library

Tahap ini memuat pustaka (library) R yang diperlukan untuk pengolahan data deret waktu, eksplorasi statistik, visualisasi, hingga pemodelan ARDL.

# REVISI: library yang tidak dipakai dibuang (nardl, ggcorrplot, corrplot, urca,
# openxlsx, MLmetrics, ARDL, MASS). e1071 juga dibuang karena menimpa
# skewness() dan kurtosis() milik moments sehingga hasilnya membingungkan.
library(readxl)
library(dplyr)
library(ggplot2)
library(moments)
library(forecast)
library(tseries)
library(dynlm)
library(dLagM)
library(lmtest)
library(car)
library(tidyr)
library(patchwork)
library(writexl)
library(dynamac)
library(ivreg)
library(sandwich)

Pustaka-pustaka di atas diaktifkan untuk mendukung seluruh tahapan analisis dari eksplorasi hingga evaluasi model.

Import Data

Tahap ini mengimpor dataset harian dari file Excel lokal ke dalam lingkungan R dan membatasi variabel menjadi peubah terikat tunggal yaitu Tinggi Muka Air (\(\text{TMA}\)) sebagai \(Y\) dan peubah bebas tunggal yaitu Curah Hujan (\(\text{PC}\)) sebagai \(X\).

dataset <- read_excel("D:/IPB UNIVERSITY/KULIAH/SEMESTER 5/MPDW/Metode Peramalan Deret Waktu/Pertemuan 2/Data_Harian_2022_2024.xlsx")
data_sub <- dataset %>% 
  dplyr::select(Date, TMA, PC) %>%
  rename(Yt = TMA, Xt = PC)
str(data_sub)

## tibble [1,096 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Date: POSIXct[1:1096], format: "2022-01-01" "2022-01-02" ...
##  $ Yt  : num [1:1096] 595 588 580 576 574 ...
##  $ Xt  : num [1:1096] 0.85 1.79 0.4 0.22 5.7 ...

Data harian periode 1 Januari 2022 - 31 Desember 2024 (\(1096\) hari) dengan fokus pada pemodelan bivariat antara Tinggi Muka Air (\(\text{TMA}\)) dan Curah Hujan (\(\text{PC}\)). Data Januari - Maret 2025 dipakai terpisah sebagai data evaluasi tambahan. Sumber Data: Data yang digunakan pada analisis ini diperoleh dari link berikut: https://github.com/fajryanti/ARDL-skripsi/tree/main.

Eksplorasi Data

Bagian ini mengeksplorasi karakteristik data melalui ringkasan statistik deskriptif, ukuran penyebaran, kemencengan, kurtosis, visualisasi deret waktu, serta korelasi antara \(\text{TMA}\) dan \(\text{PC}\).

Statistika Deskriptif

Statistika Deskriptif Dasar

Tahap ini melihat ukuran pemusatan dan penyebaran dasar dari peubah \(\text{TMA}\) dan \(\text{PC}\).

summary(data_sub[, c("Yt", "Xt")])

##        Yt              Xt          
##  Min.   :555.4   Min.   :  0.0000  
##  1st Qu.:597.7   1st Qu.:  0.7375  
##  Median :620.0   Median :  4.0000  
##  Mean   :625.3   Mean   :  6.5700  
##  3rd Qu.:650.1   3rd Qu.:  9.1325  
##  Max.   :781.2   Max.   :145.0300

Standard Deviation

Tahap ini menghitung simpangan baku untuk mengetahui tingkat keragaman data harian pada peubah \(\text{TMA}\) dan \(\text{PC}\).

sd_tma <- sd(data_sub$Yt, na.rm = TRUE)
sd_pc <- sd(data_sub$Xt, na.rm = TRUE)
cat("Standard Deviation TMA:", sd_tma, "\n")

## Standard Deviation TMA: 39.44962

cat("Standard Deviation PC:", sd_pc, "\n")

## Standard Deviation PC: 8.838259

Tinggi Muka Air (\(\text{TMA}\)) memiliki simpangan baku sebesar \(39.45\), menunjukkan sebaran data respon yang cukup bervariasi.
Curah Hujan (\(\text{PC}\)) memiliki simpangan baku sebesar \(8.84\), menunjukkan karakteristik fluktuasi harian curah hujan.

Skewness

Tahap ini menghitung koefisien kemencengan (skewness) untuk mengevaluasi tingkat kesimetrisan distribusi data peubah \(\text{TMA}\) dan \(\text{PC}\).

skew_tma <- skewness(data_sub$Yt, na.rm = TRUE)
skew_pc <- skewness(data_sub$Xt, na.rm = TRUE)
cat("Skewness TMA:", skew_tma, "\n")

## Skewness TMA: 0.5304754

cat("Skewness PC:", skew_pc, "\n")

## Skewness PC: 5.070871

Jika skewness \(\approx 0\), berarti distribusi simetris.
\(\text{TMA}\) memiliki skewness \(0.53\), yaitu miring ke kanan secara ringan.
Curah Hujan (\(\text{PC}\)) memiliki skewness \(5.07\), yang tinggi dan positif, menandakan distribusi data sangat miring ke kanan (right-skewed) akibat adanya hari-hari dengan curah hujan ekstrem.

Kurtosis

Tahap ini menghitung nilai kurtosis guna mengetahui tingkat keruncingan distribusi peubah \(\text{TMA}\) dan \(\text{PC}\).

# REVISI: moments::kurtosis() menghasilkan kurtosis biasa (normal = 3),
# jadi dikurangi 3 untuk mendapat excess kurtosis (normal = 0)
kurt_tma <- moments::kurtosis(data_sub$Yt, na.rm = TRUE) - 3
kurt_pc <- moments::kurtosis(data_sub$Xt, na.rm = TRUE) - 3
cat("Excess Kurtosis TMA:", kurt_tma, "\n")

## Excess Kurtosis TMA: -0.01559565

cat("Excess Kurtosis PC:", kurt_pc, "\n")

## Excess Kurtosis PC: 59.19679

Nilai yang ditampilkan adalah excess kurtosis (acuan distribusi normal = 0). Excess kurtosis \(\text{TMA}\) sebesar \(-0.016\) menunjukkan distribusi \(\text{TMA}\) sangat dekat dengan distribusi normal (mesokurtik), sedangkan \(\text{PC}\) sebesar \(59.2\) menunjukkan ekor yang sangat berat (leptokurtik) akibat kejadian hujan ekstrem.

Plot Time Series TMA dan PC

Tahap ini menyajikan grafik garis deret waktu untuk melihat pola pergerakan harian dari peubah \(\text{TMA}\) dan \(\text{PC}\) secara visual.

data_long <- data_sub %>%
  pivot_longer(cols = c(Yt, Xt),
               names_to = "Variabel",
               values_to = "Nilai") %>%
  mutate(Variabel = ifelse(Variabel == "Yt", "TMA (Respon)", "PC (Curah Hujan)"))

ggplot(data_long, aes(x = Date, y = Nilai)) +
  geom_line(color = "steelblue") +
  facet_wrap(~ Variabel, scales = "free_y", ncol = 1) + 
  labs(title = "Time Series Plot TMA dan PC",
       x = "Date",
       y = "Value") +
  theme_minimal()

Berdasarkan plot time series di atas, peubah Curah Hujan (\(\text{PC}\)) menunjukkan fluktuasi yang konstan di sekitar nilai rendah dengan beberapa lonjakan ekstrem yang tajam secara berkala.
Peubah Tinggi Muka Air (\(\text{TMA}\)) memperlihatkan pola naik turun tanpa tren linier yang kaku, mengindikasikan adanya dinamika merespon kondisi musim atau curah hujan.
Pada periode sekitar akhir 2023 hingga awal 2024, \(\text{TMA}\) tampak datar di nilai terendahnya. Pola ini perlu dicek apakah benar kondisi lapangan atau nilai yang terbatas/terisi otomatis oleh alat ukur, karena dapat memengaruhi uji kestasioneran dan kualitas model.

Korelasi

Pemeriksaan korelasi linear sederhana dilakukan untuk melihat sejauh mana hubungan awal antara peubah bebas dan peubah terikat.

cor_val <- cor(data_sub$Yt, data_sub$Xt, use = "complete.obs")
cat("Koefisien Korelasi antara TMA dan PC:", cor_val, "\n")

## Koefisien Korelasi antara TMA dan PC: 0.3945932

Berdasarkan perhitungan di atas, diperoleh koefisien korelasi antara Tinggi Muka Air (\(\text{TMA}\)) dan Curah Hujan (\(\text{PC}\)) sebesar \(0.3946\).
Nilai korelasi positif menunjukkan hubungan searah yang lemah hingga sedang antara curah hujan harian dan tinggi muka air.

Data Splitting

Bagian ini membagi keseluruhan dataset menjadi data pelatihan (train) dan data pengujian (test) untuk keperluan evaluasi model peramalan.

total_rows <- nrow(data_sub)
test_size <- 153

train_index <- 1:(total_rows - test_size)
test_index  <- (total_rows - test_size + 1):total_rows

train <- data_sub[train_index, ]
test  <- data_sub[test_index, ]

cat("Jumlah data train:", nrow(train), "\n")

## Jumlah data train: 943

cat("Jumlah data test:", nrow(test), "\n")

## Jumlah data test: 153

Data dibagi menjadi data train sebesar \(86\%\) dari total data atau sebanyak \(943\) data dan data test sebesar \(14\%\) dari total data atau sebanyak \(153\) data (1 Agustus - 31 Desember 2024).

Time Series Data Split Plot Plot Time Series Gabungan (TMA dan PC) Train vs Test Tahap ini memvisualisasikan pembagian data latih dan data uji pada peubah \(\text{TMA}\) dan \(\text{PC}\).

combined_data <- bind_rows(
  train %>% mutate(DataType = "Train"),
  test %>% mutate(DataType = "Test")
)

data_long_split <- combined_data %>%
  dplyr::select(Date, Yt, Xt, DataType) %>%
  pivot_longer(
    cols = c(Yt, Xt),
    names_to = "Variabel",
    values_to = "Nilai"
  ) %>%
  mutate(Variabel = ifelse(Variabel == "Yt", "TMA (Respon)", "PC (Curah Hujan)"))

ggplot(data_long_split, aes(x = Date, y = Nilai, color = DataType)) +
  geom_line() +
  facet_wrap(~ Variabel, scales = "free_y", ncol = 1) +
  scale_color_manual(
    values = c("Train" = "#00BFC4", "Test" = "#F8766D")
  ) +
  labs(
    title = "Time Series Plot (Train vs. Test)",
    x = "Date",
    y = "Value",
    color = "Data Type"
  ) +
  theme_minimal()

Kestasioneran Data dan CCF

Bagian ini menganalisis kestasioneran data melalui plot ACF, PACF, uji formal Augmented Dickey-Fuller (ADF) dan KPSS, serta hubungan lag melalui Cross-Correlation Function (CCF).

Plot ACF TMA dan PC

Tahap ini menampilkan plot Autocorrelation Function (ACF) untuk peubah \(\text{TMA}\) dan \(\text{PC}\) pada data latih.

acf_TMA <- ggAcf(train$Yt) + labs(title = "ACF TMA")
acf_PC <- ggAcf(train$Xt) + labs(title = "ACF PC")

combined_acf <- acf_TMA / acf_PC
combined_acf

Plot ACF \(\text{TMA}\) meluruh lambat, artinya data sangat persisten (nilai hari ini sangat dipengaruhi nilai hari sebelumnya). ACF \(\text{PC}\) meluruh lebih cepat. Pola peluruhan lambat ini belum otomatis berarti data tidak stasioner, sehingga kestasioneran diputuskan lewat uji formal (ADF dan KPSS) di bawah.

Plot PACF TMA dan PC

Tahap ini menampilkan plot Partial Autocorrelation Function (PACF) untuk peubah \(\text{TMA}\) dan \(\text{PC}\).

pacf_TMA <- ggPacf(train$Yt) + labs(title = "PACF TMA")
pacf_PC <- ggPacf(train$Xt) + labs(title = "PACF PC")

combined_pacf <- pacf_TMA / pacf_PC
combined_pacf

Plot PACF \(\text{TMA}\) menunjukkan lag pertama yang sangat dominan disusul beberapa lag awal yang kecil tetapi signifikan. Pola ACF meluruh dan PACF terpotong di lag-lag awal ini adalah karakteristik proses autoregresif (AR) yang stasioner dan persisten, bukan tanda ketidakstasioneran.

Uji ADF dan KPSS

Tahap ini melakukan pengujian formal untuk menguji kestasioneran data. ADF menguji \(H_0\): ada akar unit (tidak stasioner), sedangkan KPSS menguji \(H_0\): data stasioner. Keduanya dipakai bersama karena ADF kurang kuat pada data yang sangat persisten.

TMA

adf.test(train$Yt)

## Warning in adf.test(train$Yt): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train$Yt
## Dickey-Fuller = -4.0051, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

kpss.test(train$Yt, null = "Level")

## 
##  KPSS Test for Level Stationarity
## 
## data:  train$Yt
## KPSS Level = 0.68864, Truncation lag parameter = 7, p-value = 0.01458

PC

adf.test(train$Xt)

## Warning in adf.test(train$Xt): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train$Xt
## Dickey-Fuller = -7.5359, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

kpss.test(train$Xt, null = "Level")

## Warning in kpss.test(train$Xt, null = "Level"): p-value greater than printed
## p-value

## 
##  KPSS Test for Level Stationarity
## 
## data:  train$Xt
## KPSS Level = 0.24939, Truncation lag parameter = 7, p-value = 0.1

ADF pada \(\text{TMA}\) dan \(\text{PC}\) menghasilkan \(p\text{-value} < 0.01\) (R memberi peringatan bahwa nilai sebenarnya lebih kecil dari yang dicetak), sehingga tolak \(H_0\): kedua peubah stasioner pada level \(I(0)\).
KPSS pada \(\text{PC}\) tidak menolak \(H_0\) (\(p > 0.10\)), sehingga ADF dan KPSS sepakat bahwa \(\text{PC}\) stasioner. Untuk \(\text{TMA}\), KPSS menolak \(H_0\) (\(p = 0.0146 < 0.05\)) sementara ADF menolak akar unit, sehingga buktinya bertentangan. Hal ini wajar pada deret yang sangat persisten (ACF meluruh lambat): \(\text{TMA}\) paling tepat dipandang sebagai deret stasioner dengan persistensi tinggi, dan keterbatasan ini perlu disebutkan.
Model ARDL utama tetap dibangun pada level, bukan pada data hasil differencing. Alasannya: ADF menolak akar unit pada kedua peubah, tidak ada peubah yang \(I(2)\), dan kerangka ARDL dengan uji batas (bounds test) tetap valid baik untuk peubah \(I(0)\) maupun \(I(1)\), sehingga ketidakpastian pada \(\text{TMA}\) tidak membatalkan pemodelan pada level.

Plot CCF antara TMA dan PC

Tahap ini melihat korelasi silang antara peubah \(\text{TMA}\) dan \(\text{PC}\) pada berbagai lag waktu.

ccf_result <- ccf(train$Yt, train$Xt, lag.max = 10, plot = TRUE, main = "CCF antara TMA dan PC")

Pada ccf(y, x), lag positif berarti \(X\) mendahului \(Y\). Korelasi tertinggi ada di lag \(+1\), sehingga curah hujan hari sebelumnya berkaitan paling kuat dengan \(\text{TMA}\) hari ini. Perlu diingat bahwa CCF pada data level dipengaruhi autokorelasi masing-masing deret, sehingga hampir semua lag tampak signifikan. Penentuan lag optimum tetap sebaiknya dikonfirmasi dengan kriteria informasi (AIC/BIC).

First Difference

Differencing

Karena ADF menunjukkan data sudah stasioner pada level, differencing tidak wajib. Bagian ini tetap dipertahankan untuk model DLM dan Koyck sebagai pembanding pada skala perubahan harian, sedangkan ARDL utama memakai data level.

DTMA <- diff(train$Yt, differences = 1)
DPC <- diff(train$Xt, differences = 1)
df <- data.frame(DTMA, DPC)
df.ts <- ts(df)
head(df)

##         DTMA   DPC
## 1  -7.708333  0.94
## 2  -7.500000 -1.39
## 3  -4.375000 -0.18
## 4  -1.458333  5.48
## 5   9.583333  3.75
## 6 -10.208333  6.72

Plot ACF

DTMA

Tahap ini menampilkan plot ACF untuk peubah \(\text{DTMA}\) hasil differencing pertama.

# REVISI: chunk duplikat (versi tanpa dan dengan ukuran font) digabung jadi satu
acf_DTMA <- ggAcf(DTMA) +
  labs(title = "ACF DTMA") +
  theme(
    axis.text.x = element_text(size = 14),
    axis.text.y = element_text(size = 14)
  )
acf_DTMA

Plot ACF \(\text{DTMA}\) tidak meluruh perlahan sehingga stasioner dalam rata-rata. Namun autokorelasi lag-1 bernilai negatif cukup besar (sekitar \(-0.34\)), yang merupakan ciri khas over-differencing pada data yang sebenarnya sudah stasioner.

DPC

Tahap ini menampilkan plot ACF untuk peubah \(\text{DPC}\) hasil differencing pertama.

acf_DPC <- ggAcf(DPC) +
  labs(title = "ACF DPC") +
  theme(
    axis.text.x = element_text(size = 14),
    axis.text.y = element_text(size = 14)
  )
acf_DPC

Plot ACF \(\text{DPC}\) tidak meluruh perlahan sehingga stasioner dalam rata-rata, dengan autokorelasi lag-1 negatif besar (sekitar \(-0.39\)) yang juga mengindikasikan over-differencing.

Plot PACF

DTMA

Tahap ini menampilkan plot PACF untuk peubah \(\text{DTMA}\) hasil differencing pertama.

pacf_DTMA <- ggPacf(DTMA) +
  labs(title = "PACF DTMA") +
  theme(
    axis.text.x = element_text(size = 14),
    axis.text.y = element_text(size = 14)
  )
pacf_DTMA

PACF \(\text{DTMA}\) bernilai negatif dan signifikan berturut-turut sampai sekitar lag 8, lalu mengecil bertahap (tails off), bukan terpotong tegas. Karena itu lag maksimum tidak ditentukan dari “cut-off” saja, tetapi dari AIC/BIC pada bagian pemilihan orde ARDL.

DPC

Tahap ini menampilkan plot PACF untuk peubah \(\text{DPC}\) hasil differencing pertama.

pacf_DPC <- ggPacf(DPC) + 
  labs(title = "PACF DPC") +
  theme(
    axis.text.x = element_text(size = 14),
    axis.text.y = element_text(size = 14)
  )
pacf_DPC

PACF \(\text{DPC}\) juga bernilai negatif dan signifikan pada banyak lag awal, mengecil bertahap, konsisten dengan pola over-differencing.

Uji ADF

DTMA

adf.test(DTMA)

## Warning in adf.test(DTMA): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  DTMA
## Dickey-Fuller = -12.566, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

Nilai \(p\text{-value} < 0.01\) (lebih kecil dari \(\alpha = 0.05\)), sehingga tolak \(H_0\). \(\text{DTMA}\) stasioner secara uji formal.

DPC

adf.test(DPC)

## Warning in adf.test(DPC): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  DPC
## Dickey-Fuller = -15.041, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

Nilai \(p\text{-value} < 0.01\) (lebih kecil dari \(\alpha = 0.05\)), sehingga tolak \(H_0\). \(\text{DPC}\) stasioner secara uji formal.

Plot CCF (Lag Optimum Peubah Penjelas)

Tahap ini menampilkan plot Cross-Correlation Function (CCF) antara peubah \(\text{DTMA}\) dan \(\text{DPC}\) untuk menentukan lag optimum dari peubah penjelas pada model DLM.

ccf_result1 <- ggCcf(DTMA, DPC, lag.max = 9) +
  ggtitle("CCF antara DTMA dan DPC") +
  theme(
    axis.text.x = element_text(size = 14), 
    axis.text.y = element_text(size = 14)  
  )
ccf_result1

Nilai Korelasi

# REVISI: plot = FALSE lalu dicetak, supaya nilai korelasinya benar-benar tampil
ccf_nilai <- ccf(DTMA, DPC, lag.max = 9, plot = FALSE)
ccf_nilai

## 
## Autocorrelations of series 'X', by lag
## 
##     -9     -8     -7     -6     -5     -4     -3     -2     -1      0      1 
##  0.053 -0.055 -0.008  0.027 -0.026  0.010 -0.099  0.106 -0.207  0.109  0.258 
##      2      3      4      5      6      7      8      9 
## -0.131 -0.012  0.038  0.010 -0.078  0.049 -0.019  0.008

Korelasi silang terbesar ada di lag \(+1\) (\(\approx 0.26\)), yaitu \(\text{DPC}\) hari sebelumnya berkaitan dengan \(\text{DTMA}\) hari ini, sehingga dipilih lag optimum \(q_1 = 1\). Lag lain (0, 2, dan beberapa lag negatif) juga melewati batas signifikansi, tetapi lag \(+1\) paling dominan dan tetap masuk akal secara hidrologi (air hujan butuh waktu mengalir ke titik pengukuran).

Model Distributed Lag (DLM)

Pemodelan DLM

Tahap ini memodelkan hubungan keterlambatan waktu menggunakan Distributed Lag Model (DLM). Peubah respon (perubahan Tinggi Muka Air / \(\text{DTMA}\)) dipengaruhi oleh peubah penjelas (perubahan Curah Hujan / \(\text{DPC}\)) pada periode saat ini (\(t\)) serta periode sebelumnya (\(t-q\)) dengan \(q = 1\).

model_dlm <- dlm(x = as.vector(DPC), y = as.vector(DTMA), q = 1)
summary(model_dlm$model)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -108.387   -9.656   -0.435    7.632  138.664 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.01530    0.80867   0.019    0.985    
## x.t          0.71077    0.09552   7.441 2.25e-13 ***
## x.1          1.01849    0.09552  10.663  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.81 on 938 degrees of freedom
## Multiple R-squared:  0.1186, Adjusted R-squared:  0.1168 
## F-statistic: 63.14 on 2 and 938 DF,  p-value: < 2.2e-16

Persamaan Model:

\[\widehat{\text{DTMA}}_t = 0.01530 + 0.71077\,\text{DPC}_t + 1.01849\,\text{DPC}_{t-1}\]

Pengujian Koefisien Individu (Uji t): Semua koefisien peubah penjelas signifikan pada \(\alpha = 5\%\), sedangkan intersep tidak signifikan (\(p = 0.985\)).
- \(\text{DPC}_t\) (\(0.71077\)): kenaikan perubahan curah hujan harian sebesar \(1\) satuan berkaitan dengan kenaikan perubahan \(\text{TMA}\) sebesar \(0.71077\) satuan pada hari yang sama.
- \(\text{DPC}_{t-1}\) (\(1.01849\)): kenaikan perubahan curah hujan sebesar \(1\) satuan kemarin berkaitan dengan kenaikan perubahan \(\text{TMA}\) hari ini sebesar \(1.01849\) satuan. Dampak hari sebelumnya lebih besar daripada hari yang sama.
- Karena data sudah di-difference, koefisien dibaca sebagai hubungan antar perubahan harian, bukan antar level curah hujan dan level \(\text{TMA}\).
Uji Simultan (Uji F): \(F = 63.14\) dengan \(p\text{-value} < 2.2 \times 10^{-16}\), sehingga tolak \(H_0\): \(\text{DPC}_t\) dan \(\text{DPC}_{t-1}\) bersama-sama berpengaruh signifikan terhadap \(\text{DTMA}_t\).
Koefisien Determinasi (\(R^2\)): Multiple R-squared \(0.1186\) (Adjusted \(0.1168\)), artinya sekitar \(11.86\%\) keragaman perubahan \(\text{TMA}\) dijelaskan oleh curah hujan saat ini dan \(1\) hari sebelumnya. Sisanya dipengaruhi faktor lain di luar model.

Evaluasi Model DLM

Tahap ini mengevaluasi kriteria kebaikan model (Goodness of Fit) DLM melalui AIC, BIC, serta metrik kesalahan in-sample pada data latih.

cat("AIC DLM:", AIC(model_dlm$model), "\n")

## AIC DLM: 8718.745

cat("BIC DLM:", BIC(model_dlm$model), "\n")

## BIC DLM: 8738.133

GoF(model_dlm)

##             n      MAE MPE MAPE    sMAPE      MASE      MSE      MRAE    GMRAE
## model_dlm 941 15.77835 NaN  Inf 1.535554 0.5764154 613.4048 951195863 5.180658

Kriteria Informasi: AIC \(= 8718.745\) dan BIC \(= 8738.133\).
Tingkat Kesalahan (Goodness of Fit): \(n = 941\), MAE \(= 15.78\), MSE \(= 613.40\), MASE \(= 0.576\).
- MAPE dan MPE bernilai Inf/NaN karena \(\text{DTMA}\) banyak bernilai nol atau mendekati nol, sehingga tidak dilaporkan.
- Metrik di atas berada pada skala perubahan harian \(\text{TMA}\). Galat prediksi satu hari ke depan pada level \(Y_t\) sama besarnya dengan galat pada \(\Delta Y_t\), sehingga nilainya dapat dibandingkan langsung dengan model lain pada tabel akhir.

Model Autoregressive Lag (ARL/Koyck)

Pemodelan ARL

Tahap ini memodelkan hubungan dinamik menggunakan Model Koyck. Berbeda dengan DLM, model ini memasukkan lag peubah respon (\(\text{DTMA}_{t-1}\)) sebagai penjelas di samping \(\text{DPC}_t\), dengan asumsi pengaruh peubah bebas meluruh secara eksponensial. Model diestimasi dengan Instrumental Variable (IV) melalui koyckDlm().

model_koyck <- koyckDlm(x = as.vector(DPC), y = as.vector(DTMA))
summary(model_koyck$model)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -232.759  -16.173   -1.249   10.972  310.149 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.01127    1.07451   0.010    0.992    
## Y.1         -0.50644    0.04662 -10.864  < 2e-16 ***
## X.t         -2.29141    0.31430  -7.291 6.55e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.96 on 938 degrees of freedom
## Multiple R-Squared: -0.5561, Adjusted R-squared: -0.5594 
## Wald test: 61.65 on 2 and 938 DF,  p-value: < 2.2e-16

Persamaan Model:

\[\widehat{\text{DTMA}}_t = 0.01127 - 0.50644\,\text{DTMA}_{t-1} - 2.29141\,\text{DPC}_t\]

Pengujian Koefisien Individu (Uji t):
- Intersep (\(0.01127\)): tidak signifikan (\(p = 0.992\)).
- Lag respon (\(-0.50644\), \(p < 2 \times 10^{-16}\)): ini adalah parameter peluruhan Koyck (\(\lambda\)). Nilainya negatif, padahal asumsi Koyck mensyaratkan \(0 < \lambda < 1\).
- \(\text{DPC}_t\) (\(-2.29141\), \(p = 6.55 \times 10^{-13}\)): signifikan, tetapi bertanda negatif, berlawanan dengan model DLM dan logika hidrologi (hujan naik seharusnya menaikkan \(\text{TMA}\)). Tanda yang terbalik ini mengindikasikan estimasi IV tidak stabil (instrumen lemah), sehingga koefisien tersebut tidak dapat ditafsirkan secara substantif.
Uji Simultan (Wald): \(61.65\) dengan \(p\text{-value} < 2.2 \times 10^{-16}\), sehingga tolak \(H_0\).
Kelayakan Model: Nilai \(\lambda\) negatif melanggar syarat \(0 < \lambda < 1\) sehingga asumsi peluruhan eksponensial Koyck tidak sesuai untuk data ini. \(R^2\) negatif (\(-0.5561\)) juga muncul, tetapi hal itu memang lazim pada estimasi IV sehingga bukan alasan utama; alasan utamanya adalah \(\lambda\) di luar rentang dan tanda koefisien yang tidak masuk akal.

Evaluasi Model ARL

GoF(model_koyck)

##               n      MAE MPE MAPE    sMAPE      MASE      MSE       MRAE
## model_koyck 941 21.38385 NaN  Inf 1.538932 0.7811958 1082.998 2323802165
##               GMRAE
## model_koyck 7.99342

\(n = 941\), MAE \(= 21.38\), MSE \(= 1082.998\), MASE \(= 0.781\).
Dibandingkan DLM (MAE \(15.78\), MSE \(613.40\)), Koyck menghasilkan kesalahan yang lebih besar. Restriksi peluruhan eksponensial ala Koyck kurang mampu menangkap dinamika curah hujan dan tinggi muka air, sehingga pemodelan dilanjutkan ke ARDL yang lebih fleksibel.

ARDL (Level) dan Dynardl (dynamac)

Pada bagian ini ARDL dibangun pada data level karena \(\text{TMA}\) dan \(\text{PC}\) sudah stasioner. Bentuk umumnya:

\[Y_t = c + \sum_{i=1}^{p}\phi_i Y_{t-i} + \sum_{j=0}^{q}\beta_j X_{t-j} + \varepsilon_t\]

Uji Kointegrasi

modelbound <- ardlBound(data = train, formula = Yt ~ Xt, case = 2, max.p = 1, max.q = 9)

##   
## Orders being calculated with max.p = 1 and max.q = 9 ...
## 
## Autoregressive order: 10 and p-orders: 2 
## ------------------------------------------------------ 
## 
##  Breusch-Godfrey Test for the autocorrelation in residuals:
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelFull$model
## LM test = 0.22164, df1 = 1, df2 = 918, p-value = 0.6379
## 
## ------------------------------------------------------ 
## 
##  Ljung-Box Test for the autocorrelation in residuals:
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 0.0029075, df = 1, p-value = 0.957
## 
## ------------------------------------------------------ 
## 
##  Breusch-Pagan Test for the homoskedasticity of residuals:
## 
##  studentized Breusch-Pagan test
## 
## data:  modelFull$model
## BP = 203.43, df = 13, p-value < 2.2e-16
## 
## The p-value of Breusch-Pagan test for the homoskedasticity of residuals:  2.700262e-36 < 0.05!
## ------------------------------------------------------ 
## 
##  Shapiro-Wilk test of normality of residuals:
## 
##  Shapiro-Wilk normality test
## 
## data:  modelFull$model$residual
## W = 0.88265, p-value < 2.2e-16
## 
## The p-value of Shapiro-Wilk test normality of residuals:  5.754731e-26 < 0.05!
## ------------------------------------------------------ 
## 
##  PESARAN, SHIN AND SMITH (2001) COINTEGRATION TEST 
## 
##  Observations: 942 
##  Number of Regressors (k): 1 
##  Case: 2 
## 
##  ------------------------------------------------------ 
##  -                       F-test                       - 
##  ------------------------------------------------------ 
##                  <------- I(0) ------------ I(1) -----> 
##  10% critical value       3.02            3.51 
##  5% critical value        3.62            4.16 
##  1% critical value        4.94            5.58 
##  
## 
##  F-statistic = 41.6397821191344 
##   
##  ------------------------------------------------------ 
##  F-statistic note: Asymptotic critical values used. 
##  
## ------------------------------------------------------ 
## 
##  Ramsey's RESET Test for model specification:
## 
##  RESET test
## 
## data:  modelECM$model
## RESET = 18.937, df1 = 1, df2 = 920, p-value = 1.502e-05
## 
## the p-value of RESET test:  1.502293e-05 < 0.05!
## ------------------------------------------------------

## ------------------------------------------------------ 
## Error Correction Model Output: 
## 
## Time series regression with "ts" data:
## Start = 10, End = 942
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -131.649  -10.959   -2.548    6.419  135.043 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## ec.1  -0.15640    0.01398 -11.189  < 2e-16 ***
## dXt.t  0.75245    0.08956   8.402  < 2e-16 ***
## dXt.1  0.33910    0.10011   3.387 0.000736 ***
## dYt.1 -0.43833    0.03002 -14.604  < 2e-16 ***
## dYt.2 -0.28542    0.03332  -8.565  < 2e-16 ***
## dYt.3 -0.21422    0.03400  -6.301 4.57e-10 ***
## dYt.4 -0.13794    0.03446  -4.003 6.77e-05 ***
## dYt.5 -0.15754    0.03397  -4.638 4.02e-06 ***
## dYt.6 -0.12855    0.03395  -3.787 0.000163 ***
## dYt.7 -0.10040    0.03357  -2.991 0.002856 ** 
## dYt.8 -0.04558    0.03250  -1.403 0.161080    
## dYt.9  0.02097    0.02934   0.715 0.474982    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.23 on 921 degrees of freedom
## Multiple R-squared:  0.3654, Adjusted R-squared:  0.3571 
## F-statistic: 44.19 on 12 and 921 DF,  p-value: < 2.2e-16
## 
## ------------------------------------------------------ 
## Long-run coefficients: 
##       Yt.1       Xt.1 
## -0.1564025  1.2918485 
##

Hasil Uji Kointegrasi Pesaran, Shin and Smith (2001) dengan \(942\) pengamatan, \(k = 1\), dan kasus \(2\): F-statistik \(= 41.64\). Nilai kritis batas atas \(I(1)\) adalah \(3.51\) (10%), \(4.16\) (5%), dan \(5.58\) (1%). Karena \(41.64\) jauh melampaui batas atas, tolak \(H_0\) dan disimpulkan ada hubungan jangka panjang antara \(\text{TMA}\) dan \(\text{PC}\).

Catatan penting:

Karena kedua peubah sudah \(I(0)\), uji batas (bounds test) memang hampir pasti menolak \(H_0\), sehingga hasil ini tidak terlalu informatif secara konsep. Yang lebih penting adalah kecepatan penyesuaian pada model ECM di bawah.
Diagnostik pada output uji ini menunjukkan masalah: Breusch-Pagan menolak homoskedastisitas (\(p < 2.2 \times 10^{-16}\)), Shapiro-Wilk menolak normalitas residual, dan Ramsey RESET menolak spesifikasi model linear (\(p = 1.5 \times 10^{-5}\)). Ketiganya dibahas ulang pada uji asumsi.

Pemilihan Orde ARDL

Orde \(p\) (lag \(Y\)) dan \(q\) (lag \(X\)) dipilih dengan AIC pada semua kombinasi. Semua kombinasi dihitung pada sampel yang sama (baris awal yang kosong akibat lag dibuang) agar AIC dapat dibandingkan.

# fungsi untuk membuat kolom lag Y dan X
buat_lag <- function(dat, p, q) {
  out <- dat
  for (i in 1:p) out[[paste0("Y_lag", i)]] <- dplyr::lag(dat$Yt, i)
  if (q > 0) {
    for (j in 1:q) out[[paste0("X_lag", j)]] <- dplyr::lag(dat$Xt, j)
  }
  out
}

# fungsi untuk membuat formula ARDL(p, q)
buat_formula <- function(p, q) {
  rhs <- c(paste0("Y_lag", 1:p), "Xt")
  if (q > 0) rhs <- c(rhs, paste0("X_lag", 1:q))
  as.formula(paste("Yt ~", paste(rhs, collapse = " + ")))
}

p_max <- 10
q_max <- 3

train_lag_max <- buat_lag(train, p_max, q_max)[-(1:p_max), ]

kandidat <- expand.grid(p = 1:p_max, q = 1:q_max)
kandidat$AIC <- NA
kandidat$BIC <- NA

for (i in 1:nrow(kandidat)) {
  m <- lm(buat_formula(kandidat$p[i], kandidat$q[i]), data = train_lag_max)
  kandidat$AIC[i] <- AIC(m)
  kandidat$BIC[i] <- BIC(m)
}

head(kandidat[order(kandidat$AIC), ], 5)

##     p q      AIC      BIC
## 19  9 2 8365.736 8433.474
## 20 10 2 8367.228 8439.804
## 18  8 2 8367.347 8430.246
## 29  9 3 8367.593 8440.169
## 30 10 3 8369.046 8446.461

p_opt <- kandidat$p[which.min(kandidat$AIC)]
q_opt <- kandidat$q[which.min(kandidat$AIC)]
cat("Orde terpilih (AIC): p =", p_opt, ", q =", q_opt, "\n")

## Orde terpilih (AIC): p = 9 , q = 2

Berdasarkan AIC terkecil, dipilih ARDL(\(9\), \(2\)). Tiga model teratas berselisih AIC kurang dari 2 (8365.7; 8367.2; 8367.3), sehingga ARDL(9, 2) tidak jauh lebih baik daripada ARDL(8, 2). Pada model terpilih, koefisien \(Y_{t-5}\) sampai \(Y_{t-8}\) juga tidak signifikan. Sebagai pembanding, berikut tiga model teratas menurut BIC yang lebih menghukum banyaknya parameter:

kandidat[order(kandidat$BIC), ][1:3, ]

##    p q      AIC      BIC
## 14 4 2 8383.436 8426.982
## 18 8 2 8367.347 8430.246
## 16 6 2 8377.310 8430.532

BIC memilih ARDL(4, 2) (BIC \(8426.98\), AIC \(8383.44\)), yaitu model yang jauh lebih hemat parameter daripada ARDL(9, 2) (BIC \(8433.47\)). Selisih BIC-nya kecil, tetapi AIC-nya lebih besar sekitar \(17.7\), sehingga kedua model layak. ARDL(9, 2) dipertahankan sebagai model utama karena AIC-nya terkecil, dan ARDL(4, 2) dapat dijadikan analisis sensitivitas: ubah p_opt <- 4 lalu knit ulang untuk melihat apakah kesimpulan berubah.

Model

Jangka Pendek (ECM dengan dynardl)

Model koreksi galat (Error Correction Model, ECM) ekuivalen dengan ARDL(\(9\), \(2\)) pada level. Data yang dimasukkan adalah data level (\(\text{Yt}\) dan \(\text{Xt}\)), dan dynardl sendiri yang membentuk selisihnya.

train_df <- as.data.frame(train[, c("Yt", "Xt")])

lagdiffs_list <- list()
if (p_opt > 1) lagdiffs_list[["Yt"]] <- 1:(p_opt - 1)
if (q_opt > 1) lagdiffs_list[["Xt"]] <- 1:(q_opt - 1)
if (length(lagdiffs_list) == 0) lagdiffs_list <- NULL

modelardl <- dynardl(
  Yt ~ Xt,
  data = train_df,
  lags = list("Yt" = 1, "Xt" = 1),
  diffs = c("Xt"),
  lagdiffs = lagdiffs_list,
  ec = TRUE,
  simulate = FALSE
)

## [1] "Error correction (EC) specified; dependent variable to be run in differences."

summary(modelardl$model)

## 
## Call:
## lm(formula = as.formula(paste(paste(dvnamelist), "~", paste(colnames(IVs), 
##     collapse = "+"), collapse = " ")))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -131.177  -11.452   -2.489    6.528  135.374 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 88.08571   13.24058   6.653 4.93e-11 ***
## l.1.Yt      -0.15416    0.02154  -7.158 1.67e-12 ***
## ld.1.Yt     -0.44144    0.03412 -12.937  < 2e-16 ***
## ld.2.Yt     -0.29015    0.03658  -7.932 6.24e-15 ***
## ld.3.Yt     -0.21879    0.03628  -6.030 2.37e-09 ***
## ld.4.Yt     -0.14362    0.03597  -3.993 7.04e-05 ***
## ld.5.Yt     -0.16226    0.03503  -4.633 4.13e-06 ***
## ld.6.Yt     -0.13388    0.03438  -3.894 0.000106 ***
## ld.7.Yt     -0.10745    0.03290  -3.266 0.001133 ** 
## ld.8.Yt     -0.05568    0.02963  -1.879 0.060530 .  
## d.1.Xt       0.75176    0.09109   8.253 5.31e-16 ***
## l.1.Xt       1.29517    0.12258  10.566  < 2e-16 ***
## ld.1.Xt      0.33221    0.10064   3.301 0.001000 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.25 on 921 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.3645, Adjusted R-squared:  0.3562 
## F-statistic: 44.03 on 12 and 921 DF,  p-value: < 2.2e-16

Cara membaca hasil:

Koefisien l.1.Yt adalah kecepatan penyesuaian (speed of adjustment) sebesar \(-0.154\). Nilainya negatif dan berada di antara \(-1\) dan \(0\), sehingga sistem stabil dan kembali ke keseimbangan. Sekitar \(15\%\) ketidakseimbangan terkoreksi setiap hari (waktu paruh sekitar \(4\) hari), dan nilai ini konsisten dengan ECM dari ardlBound (\(-0.156\)).
Koefisien d.1.Xt sebesar \(0.752\) adalah dampak jangka pendek: kenaikan curah hujan \(1\) satuan pada hari yang sama langsung menaikkan \(\text{TMA}\) sekitar \(0.75\) satuan. Efek totalnya berkembang bertahap menuju pengganda jangka panjang di bawah.
\(R^2\) pada model ECM dihitung untuk respon \(\Delta Y_t\), jadi tidak boleh dibandingkan langsung dengan \(R^2\) model level.

Jangka Panjang

Pengganda jangka panjang dihitung dari koefisien ARDL level (bukan hanya koefisien \(X_t\)):

\[\text{LRM} = \frac{\sum_j \beta_j}{1-\sum_i \phi_i}\]

# model ARDL level (dipakai juga untuk diagnostik dan peramalan)
n_buang <- max(p_opt, q_opt)
train_lag <- buat_lag(train, p_opt, q_opt)[-(1:n_buang), ]
model_ardl <- lm(buat_formula(p_opt, q_opt), data = train_lag)
summary(model_ardl)

## 
## Call:
## lm(formula = buat_formula(p_opt, q_opt), data = train_lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -131.177  -11.452   -2.489    6.528  135.374 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 88.08571   13.24058   6.653 4.93e-11 ***
## Y_lag1       0.40440    0.03273  12.357  < 2e-16 ***
## Y_lag2       0.15128    0.03413   4.433 1.04e-05 ***
## Y_lag3       0.07136    0.03381   2.111   0.0351 *  
## Y_lag4       0.07517    0.03362   2.236   0.0256 *  
## Y_lag5      -0.01864    0.03360  -0.555   0.5793    
## Y_lag6       0.02838    0.03344   0.849   0.3962    
## Y_lag7       0.02642    0.03334   0.792   0.4283    
## Y_lag8       0.05177    0.03308   1.565   0.1179    
## Y_lag9       0.05568    0.02963   1.879   0.0605 .  
## Xt           0.75176    0.09109   8.253 5.31e-16 ***
## X_lag1       0.87562    0.09995   8.761  < 2e-16 ***
## X_lag2      -0.33221    0.10064  -3.301   0.0010 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.25 on 921 degrees of freedom
## Multiple R-squared:  0.7343, Adjusted R-squared:  0.7309 
## F-statistic: 212.1 on 12 and 921 DF,  p-value: < 2.2e-16

koef <- coef(model_ardl)
phi_total <- sum(koef[grep("^Y_lag", names(koef))])
beta_total <- sum(koef[c("Xt", grep("^X_lag", names(koef), value = TRUE))])

lrm <- beta_total / (1 - phi_total)
kec_penyesuaian <- -(1 - phi_total)

cat("Jumlah koefisien lag Y :", phi_total, "\n")

## Jumlah koefisien lag Y : 0.8458407

cat("Jumlah koefisien X     :", beta_total, "\n")

## Jumlah koefisien X     : 1.295166

cat("Pengganda jangka panjang (LRM):", lrm, "\n")

## Pengganda jangka panjang (LRM): 8.401475

cat("Kecepatan penyesuaian (EC)    :", kec_penyesuaian, "\n")

## Kecepatan penyesuaian (EC)    : -0.1541593

cat("AIC ARDL level:", AIC(model_ardl), "\n")

## AIC ARDL level: 8374.42

Pengganda jangka panjang sebesar \(8.401\): kenaikan tetap curah hujan sebesar \(1\) satuan berkaitan dengan kenaikan tinggi muka air keseimbangan sebesar sekitar \(8.401\) satuan.
Kecepatan penyesuaian model level (\(-0.154\)) sama dengan koefisien l.1.Yt pada ECM (\(-0.154\)), sehingga kedua bentuk konsisten (ECM dan ARDL level adalah penulisan ulang dari model yang sama).
Pengganda jangka panjang adalah rasio, dan karena \(\sum_i \phi_i = 0.846\) mendekati \(1\), nilainya sensitif terhadap perubahan kecil pada koefisien. Angka ini belum disertai selang kepercayaan, sehingga sebaiknya dibaca sebagai perkiraan kasar.
Adjusted \(R^2\) model level sebesar \(0.731\) tinggi terutama karena \(\text{TMA}\) sangat persisten (ada lag \(Y\) di ruas kanan), bukan semata karena curah hujan.

Uji Asumsi

Uji asumsi dilakukan pada residual model ARDL level.

res_ardl <- residuals(model_ardl)

Uji Autokorelasi Residual

Karena ada lag peubah respon di ruas kanan, uji Breusch-Godfrey lebih tepat dipakai daripada uji Durbin-Watson. Ljung-Box tetap ditampilkan sebagai pelengkap.

bgtest(model_ardl, order = 10)

## 
##  Breusch-Godfrey test for serial correlation of order up to 10
## 
## data:  model_ardl
## LM test = 2.5411, df = 10, p-value = 0.9903

Box.test(res_ardl, lag = 20, type = "Ljung-Box", fitdf = p_opt)

## 
##  Box-Ljung test
## 
## data:  res_ardl
## X-squared = 9.0678, df = 11, p-value = 0.6156

Hipotesis: \(H_0\) tidak ada autokorelasi residual, \(H_1\) ada autokorelasi. Uji Breusch-Godfrey menghasilkan \(p = 0.990\) dan Ljung-Box \(p = 0.616\), keduanya lebih besar dari \(0.05\), sehingga gagal tolak \(H_0\): tidak ada bukti autokorelasi residual, dan jumlah lag pada model sudah memadai.

Uji Normalitas (Jarque-Bera)

jarque.bera.test(res_ardl)

## 
##  Jarque Bera Test
## 
## data:  res_ardl
## X-squared = 1706.8, df = 2, p-value < 2.2e-16

\(H_0\): residual berdistribusi normal. Statistik Jarque-Bera \(= 1706.8\) dengan \(p < 2.2 \times 10^{-16}\), sehingga tolak \(H_0\): residual tidak normal (ekor berat akibat hari hujan ekstrem). Uji t dan selang kepercayaan hanya berlaku secara aproksimasi.

Uji Heteroskedastisitas (Breusch-Pagan)

# REVISI: bptest langsung pada model, bukan pada regresi res ~ fitted
bptest(model_ardl)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_ardl
## BP = 200.54, df = 12, p-value < 2.2e-16

\(H_0\): ragam residual konstan. Hasil \(BP = 200.54\) dengan \(p < 2.2 \times 10^{-16}\), sehingga tolak \(H_0\): ragam residual tidak konstan (heteroskedastis), sesuai dengan temuan pada ardlBound. Karena itu inferensi koefisien memakai galat baku yang tahan heteroskedastisitas (HC3):

coeftest(model_ardl, vcov. = vcovHC(model_ardl, type = "HC3"))

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 88.085707  15.912786  5.5355 4.046e-08 ***
## Y_lag1       0.404404   0.050797  7.9612 4.995e-15 ***
## Y_lag2       0.151285   0.044197  3.4229 0.0006468 ***
## Y_lag3       0.071362   0.043955  1.6235 0.1048196    
## Y_lag4       0.075169   0.036203  2.0763 0.0381401 *  
## Y_lag5      -0.018639   0.033868 -0.5503 0.5822224    
## Y_lag6       0.028383   0.036833  0.7706 0.4411510    
## Y_lag7       0.026423   0.035895  0.7361 0.4618383    
## Y_lag8       0.051769   0.039092  1.3243 0.1857337    
## Y_lag9       0.055683   0.034950  1.5932 0.1114477    
## Xt           0.751756   0.244716  3.0720 0.0021891 ** 
## X_lag1       0.875623   0.532310  1.6449 0.1003218    
## X_lag2      -0.332213   0.150298 -2.2104 0.0273248 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dengan galat baku HC3, koefisien \(X_t\) (\(p = 0.002\)) dan \(X_{t-2}\) (\(p = 0.027\)) tetap signifikan, tetapi \(X_{t-1}\) (\(p = 0.100\)) dan \(Y_{t-3}\) (\(p = 0.105\)) tidak lagi signifikan pada \(5\%\). Jadi klaim bahwa curah hujan hari sebelumnya berpengaruh signifikan menjadi lebih lemah setelah heteroskedastisitas diperhitungkan, walaupun jumlah koefisien curah hujan (\(1.295\)) tetap positif.

Uji Spesifikasi (RESET)

resettest(model_ardl, power = 2:3, type = "fitted")

## 
##  RESET test
## 
## data:  model_ardl
## RESET = 25.464, df1 = 2, df2 = 919, p-value = 1.724e-11

\(H_0\): bentuk fungsional linear sudah memadai. Hasil \(RESET = 25.46\) dengan \(p = 1.7 \times 10^{-11}\), sehingga tolak \(H_0\): spesifikasi linear belum memadai. Hubungan curah hujan dan \(\text{TMA}\) kemungkinan tidak linear (mis. efek hujan ekstrem berbeda), sehingga transformasi seperti \(\log(1+X_t)\) layak dicoba pada pengembangan berikutnya.

Multikolinearitas

vif(model_ardl)

##   Y_lag1   Y_lag2   Y_lag3   Y_lag4   Y_lag5   Y_lag6   Y_lag7   Y_lag8 
## 3.718772 4.046942 3.976778 3.933900 3.935929 3.901369 3.882894 3.824497 
##   Y_lag9       Xt   X_lag1   X_lag2 
## 3.068784 1.300768 1.565318 1.586148

Nilai VIF berkisar antara \(1.30\) sampai \(4.05\), jauh di bawah batas kritis \(10\) (bahkan di bawah \(5\)), sehingga tidak ada masalah multikolinearitas serius antar lag.

Forecasting

Peramalan yang dipakai adalah peramalan satu hari ke depan (one-step-ahead): \(\hat{Y}_t\) dihitung dari \(Y_{t-1},\dots,Y_{t-p}\) aktual dan \(X_t, X_{t-1},\dots\) aktual. Model pembanding sederhana (naif) \(\hat{Y}_t = Y_{t-1}\) ditambahkan agar akurasi dapat dinilai secara wajar, karena \(\text{TMA}\) sangat persisten sehingga MAPE kecil saja tidak cukup membuktikan model bagus.

# fungsi ukuran akurasi
ukur_akurasi <- function(aktual, prediksi) {
  data.frame(
    MAE  = mean(abs(aktual - prediksi)),
    RMSE = sqrt(mean((aktual - prediksi)^2)),
    MAPE = mean(abs((aktual - prediksi) / aktual)) * 100
  )
}

Data Training

pred_train <- fitted(model_ardl)
aktual_train <- train_lag$Yt
naif_train <- train_lag$Y_lag1

plot(train_lag$Date, aktual_train, type = "l", col = "black",
     main = "Train: Actual vs Predicted TMA", xlab = "Tanggal", ylab = "TMA")
lines(train_lag$Date, pred_train, col = "blue")
legend("topleft", legend = c("Actual", "Predicted"), col = c("black", "blue"), lty = 1)

tabel_train <- rbind(
  ARDL = ukur_akurasi(aktual_train, pred_train),
  Naif = ukur_akurasi(aktual_train, naif_train)
)
round(tabel_train, 3)

##         MAE   RMSE  MAPE
## ARDL 14.206 21.098 2.204
## Naif 16.647 26.466 2.571

Pada data training, ARDL memiliki MAPE \(2.2\%\) dan RMSE \(21.1\), sedangkan model naif RMSE-nya \(26.47\). RMSE ARDL lebih kecil dibandingkan model naif. RMSE training (\(21.1\)) hampir sama dengan simpangan baku residual model (\(21.25\)), yang menandakan prediksi sudah dihitung konsisten dengan model yang diestimasi. Pada grafik, prediksi tampak melonjak tajam mendekati \(780\) di sekitar hari hujan ekstrem (\(145\) mm), yang menunjukkan model sensitif terhadap pencilan curah hujan.

Data Testing

Data pengujian memiliki \(153\) observasi (1 Agustus - 31 Desember 2024). Lag untuk beberapa hari pertama diambil dari data training yang sebenarnya teramati.

# lag dibuat dari seluruh data agar hari-hari awal data test memakai nilai aktual sebelumnya
full_lag <- buat_lag(data_sub, p_opt, q_opt)
test_lag <- full_lag[test_index, ]

aktual_test <- test_lag$Yt
pred_ardl_test <- predict(model_ardl, newdata = test_lag)
naif_test <- test_lag$Y_lag1

# DLM dan Koyck (model pada skala selisih), dikembalikan ke level TMA
dX <- c(NA, diff(data_sub$Xt))
dY <- c(NA, diff(data_sub$Yt))
Y_prev <- data_sub$Yt[test_index - 1]

b_dlm <- unname(coef(model_dlm$model))      # urutan: intersep, x.t, x.1
pred_dlm_test <- Y_prev + b_dlm[1] + b_dlm[2] * dX[test_index] + b_dlm[3] * dX[test_index - 1]

b_koyck <- unname(coef(model_koyck$model))  # urutan: intersep, Y.1, X.t
pred_koyck_test <- Y_prev + b_koyck[1] + b_koyck[2] * dY[test_index - 1] + b_koyck[3] * dX[test_index]

plot(test_lag$Date, aktual_test, type = "l", col = "black",
     main = "Test: Actual vs Predicted TMA", xlab = "Tanggal", ylab = "TMA")
lines(test_lag$Date, pred_ardl_test, col = "blue")
legend("topleft", legend = c("Actual", "ARDL"), col = c("black", "blue"), lty = 1)

tabel_test <- rbind(
  ARDL  = ukur_akurasi(aktual_test, pred_ardl_test),
  DLM   = ukur_akurasi(aktual_test, pred_dlm_test),
  Koyck = ukur_akurasi(aktual_test, pred_koyck_test),
  Naif  = ukur_akurasi(aktual_test, naif_test)
)
round(tabel_test, 3)

##          MAE   RMSE  MAPE
## ARDL   9.820 16.555 1.545
## DLM   10.655 19.696 1.667
## Koyck 18.705 33.842 2.937
## Naif  10.052 19.192 1.562

Tabel di atas membandingkan semua model pada data testing dengan skala yang sama (galat prediksi satu hari ke depan pada level \(\text{TMA}\)).
RMSE ARDL sebesar \(16.56\), DLM \(19.7\), Koyck \(33.84\), dan model naif \(19.19\).
RMSE ARDL lebih kecil dibandingkan model naif, dan lebih kecil dibandingkan DLM dan Koyck.
Ukuran lain memberi gambaran yang lebih hati-hati. MAE ARDL (\(9.82\)) hanya sedikit lebih kecil dari model naif (\(10.05\)), dan MAPE keduanya hampir sama (\(1.55\%\) vs \(1.56\%\)). Keunggulan ARDL terutama muncul pada RMSE, yaitu pada hari-hari dengan lonjakan besar, bukan pada hari biasa.
Karena level \(\text{TMA}\) berada di sekitar \(600\), MAPE yang kecil mudah dicapai oleh model sederhana sekalipun. Sebutan “sangat akurat” kurang tepat. Kesimpulan yang lebih jujur: ARDL memiliki RMSE lebih kecil, terutama saat lonjakan \(\text{TMA}\), tetapi keunggulannya atas model naif belum signifikan secara statistik (lihat uji Diebold-Mariano di bawah).
DLM tidak lebih baik daripada model naif (RMSE \(19.7\) vs \(19.19\)), dan Koyck jauh lebih buruk.
Pada grafik, ARDL mengikuti lonjakan \(\text{TMA}\) dengan cukup baik, termasuk puncak sekitar \(745\) pada akhir November yang diprediksi sekitar \(720\).

Untuk menguji apakah selisih ARDL dan model naif nyata secara statistik, dipakai uji Diebold-Mariano pada galat kuadrat, dengan \(H_1\): ARDL lebih akurat daripada model naif.

dm.test(aktual_test - pred_ardl_test, aktual_test - naif_test,
        alternative = "less", h = 1, power = 2)

## 
##  Diebold-Mariano Test
## 
## data:  aktual_test - pred_ardl_testaktual_test - naif_test
## DM = -1.5974, Forecast horizon = 1, Loss function power = 2, p-value =
## 0.05612
## alternative hypothesis: less

Hasilnya \(DM = -1.60\) dengan \(p = 0.056\). Nilai ini sedikit di atas \(0.05\), sehingga keunggulan ARDL atas model naif belum signifikan pada taraf 5% (hanya mendekati signifikan pada taraf 10%). Perbaikan RMSE sebesar sekitar \(14\%\) pada data testing belum dapat dipastikan bukan kebetulan.

Prediksi 2025

Model yang sudah dilatih dipakai untuk peramalan satu hari ke depan pada data Januari - Maret 2025 (\(90\) hari). Sebanyak \(9\) hari terakhir 2024 dipakai hanya untuk membentuk lag pada hari-hari pertama.

data_2025 <- read_excel("D:/IPB UNIVERSITY/KULIAH/SEMESTER 5/MPDW/Metode Peramalan Deret Waktu/Pertemuan 2/Data_Jan25_Mar25.xlsx")
data_2025 <- data_2025 %>%
  transmute(Date = as.Date(Date), Yt = TMA, Xt = PC)
cat("Jumlah data 2025:", nrow(data_2025), "\n")

## Jumlah data 2025: 90

# REVISI: tail(..., 120) dihapus karena datanya hanya 90 baris
n_lag <- max(p_opt, q_opt)
gabung <- bind_rows(
  data_sub %>% mutate(Date = as.Date(Date)) %>% tail(n_lag),
  data_2025
)
gabung_lag <- buat_lag(gabung, p_opt, q_opt)[-(1:n_lag), ]

pred_2025 <- predict(model_ardl, newdata = gabung_lag)

hasil_prediksi <- data.frame(
  date = gabung_lag$Date,
  TMA_actual = gabung_lag$Yt,
  TMA_pred = pred_2025,
  TMA_naif = gabung_lag$Y_lag1
)

write_xlsx(hasil_prediksi, path = "hasil_peramalan.xlsx")

plot(hasil_prediksi$date, hasil_prediksi$TMA_actual, type = "l", col = "black",
     ylab = "TMA", xlab = "Tanggal", main = "TMA Aktual vs Prediksi (Jan - Mar 2025)")
lines(hasil_prediksi$date, hasil_prediksi$TMA_pred, col = "blue")
legend("topleft", legend = c("Aktual", "Prediksi ARDL"), col = c("black", "blue"), lty = 1)

tabel_2025 <- rbind(
  ARDL = ukur_akurasi(hasil_prediksi$TMA_actual, hasil_prediksi$TMA_pred),
  Naif = ukur_akurasi(hasil_prediksi$TMA_actual, hasil_prediksi$TMA_naif)
)
round(tabel_2025, 3)

##         MAE  RMSE  MAPE
## ARDL 11.093 16.64 1.659
## Naif 11.354 22.64 1.660

Pada periode Januari - Maret 2025, ARDL memiliki RMSE \(16.64\) dan MAPE \(1.66\%\), sedangkan model naif memiliki RMSE \(22.64\) dan MAPE \(1.66\%\). Seperti pada data testing, MAE dan MAPE ARDL hampir identik dengan model naif (MAE \(11.09\) vs \(11.35\)). Keunggulan RMSE datang dari lonjakan awal Maret, yang tetap diprediksi jauh lebih rendah dari nilai aktualnya (sekitar \(725\) vs \(795\)).

Catatan: ini adalah evaluasi peramalan satu hari ke depan dengan curah hujan aktual sebagai masukan. Untuk peramalan jauh ke depan yang sesungguhnya, nilai curah hujan masa depan juga harus diramalkan atau diberikan sebagai skenario.

Kesimpulan dan Saran

Kesimpulan

Analisis deret waktu bivariat antara Tinggi Muka Air (TMA) dan Curah Hujan (PC) memberikan kesimpulan berikut:

Kestasioneran: Uji ADF menunjukkan \(\text{TMA}\) dan \(\text{PC}\) stasioner pada level. KPSS sepakat untuk \(\text{PC}\), tetapi menolak kestasioneran \(\text{TMA}\) (\(p = 0.0146\)), sehingga \(\text{TMA}\) dipandang sebagai deret sangat persisten dengan bukti yang tidak sepenuhnya seragam. Differencing menimbulkan ciri over-differencing (ACF lag-1 negatif besar), sehingga ARDL utama dibangun pada level, yang tetap valid dalam kerangka uji batas.
Model Distributed Lag (DLM): DLM dengan \(q = 1\) pada data selisih menunjukkan bahwa perubahan \(\text{TMA}\) dipengaruhi perubahan curah hujan hari yang sama (\(0.71077\)) dan hari sebelumnya (\(1.01849\)), dengan \(R^2\) hanya \(11.86\%\), MAE \(= 15.78\), dan MSE \(= 613.40\) pada data latih.
Model Koyck: Estimasi Koyck tidak layak dipakai karena parameter peluruhan \(\lambda = -0.506\) berada di luar rentang \(0 < \lambda < 1\) dan koefisien curah hujan bertanda negatif (tidak masuk akal secara hidrologi). MAE \(= 21.38\) dan MSE \(= 1082.00\) juga lebih buruk dari DLM.
Model ARDL: ARDL level terpilih adalah ARDL(\(9\), \(2\)) dengan kecepatan penyesuaian \(-0.154\) per hari dan pengganda jangka panjang \(8.401\). Uji batas Pesaran-Shin-Smith menunjukkan adanya hubungan jangka panjang, walaupun uji ini kurang informatif bila peubah sudah stasioner.
Kinerja peramalan: Pada data testing, RMSE ARDL sebesar \(16.56\) dibandingkan \(19.19\) untuk model naif, sedangkan DLM dan Koyck lebih buruk. Namun MAE dan MAPE ARDL hanya sedikit lebih baik daripada model naif, dan uji Diebold-Mariano menghasilkan \(p = 0.056\), sehingga keunggulan ARDL atas model naif belum signifikan pada taraf 5%. Perbaikannya bersifat moderat dan terutama muncul saat lonjakan \(\text{TMA}\). Pola yang sama terlihat pada data Januari - Maret 2025. Lonjakan tertinggi masih cenderung diprediksi lebih rendah dari nilai aktualnya.
Keterbatasan: Residual tidak normal (Jarque-Bera dan Shapiro-Wilk), ragam residual tidak konstan (Breusch-Pagan pada ardlBound), dan uji RESET menunjukkan kemungkinan hubungan tidak linear. Karena itu inferensi (uji t, selang kepercayaan) perlu ditafsirkan hati-hati, dan galat baku yang tahan heteroskedastisitas sebaiknya dipakai. Dengan galat baku HC3, pengaruh curah hujan hari sebelumnya (\(X_{t-1}\)) tidak lagi signifikan pada \(5\%\).

Saran

Bandingkan selalu dengan model naif sebelum menyebut model akurat, dan laporkan RMSE/MAE di samping MAPE.
Coba transformasi curah hujan (mis. \(\log(1+X_t)\)) atau efek non-linear untuk menangani hujan ekstrem.
Pertimbangkan peubah penjelas tambahan yang sudah tersedia di data 2025 (suhu, kelembapan, kecepatan angin).
Periksa kualitas data \(\text{TMA}\) pada periode yang tampak datar (akhir 2023 - awal 2024).
Model ini dapat menjadi dasar awal sistem peringatan dini, tetapi belum layak dijadikan satu-satunya dasar keputusan mitigasi banjir sebelum kinerjanya pada lonjakan \(\text{TMA}\) ekstrem dievaluasi dan diperbaiki.

Tugas Kuliah 2 - Analisis ARDL Data Curah Hujan

Arifatu Dzikra

2026-09-29

Library

Import Data

Eksplorasi Data

Statistika Deskriptif

Statistika Deskriptif Dasar

Standard Deviation

Skewness

Kurtosis

Plot Time Series TMA dan PC

Korelasi

Data Splitting

Kestasioneran Data dan CCF

Plot ACF TMA dan PC

Plot PACF TMA dan PC

Uji ADF dan KPSS

TMA

PC

Plot CCF antara TMA dan PC

First Difference

Differencing

Plot ACF

DTMA

DPC

Plot PACF

DTMA

DPC

Uji ADF

DTMA

DPC

Plot CCF (Lag Optimum Peubah Penjelas)

Model Distributed Lag (DLM)

Pemodelan DLM

Evaluasi Model DLM

Model Autoregressive Lag (ARL/Koyck)

Pemodelan ARL

Evaluasi Model ARL

ARDL (Level) dan Dynardl (dynamac)

Uji Kointegrasi

Pemilihan Orde ARDL

Model

Jangka Pendek (ECM dengan dynardl)

Jangka Panjang

Uji Asumsi

Uji Autokorelasi Residual

Uji Normalitas (Jarque-Bera)

Uji Heteroskedastisitas (Breusch-Pagan)

Uji Spesifikasi (RESET)

Multikolinearitas

Forecasting

Data Training

Data Testing

Prediksi 2025

Kesimpulan dan Saran

Kesimpulan

Saran