1 Library

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork)
library(moments)
library(corrplot)
library(forecast)
library(tseries)
library(urca)
library(dLagM)
library(dynamac)
library(lmtest)
library(writexl)

2 Import Data

dataset <- read_excel("C:/Users/hp/Downloads/Data_Harian_2022_2024_mpdw.xlsx")
dataset <- dataset %>% arrange(Date)
str(dataset)

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

2.1 Data

head(dataset, 10)

## # A tibble: 10 × 3
##    Date                  TMA    PC
##    <dttm>              <dbl> <dbl>
##  1 2022-01-01 00:00:00  595.  0.85
##  2 2022-01-02 00:00:00  588.  1.79
##  3 2022-01-03 00:00:00  580.  0.4 
##  4 2022-01-04 00:00:00  576.  0.22
##  5 2022-01-05 00:00:00  574.  5.7 
##  6 2022-01-06 00:00:00  584.  9.45
##  7 2022-01-07 00:00:00  574. 16.2 
##  8 2022-01-08 00:00:00  586.  3.37
##  9 2022-01-09 00:00:00  570.  4.79
## 10 2022-01-10 00:00:00  568.  9.14

Dataset memuat tiga kolom: Date (tanggal pengamatan harian), TMA (Tinggi Muka Air, dalam cm), dan PC (curah hujan/precipitation, dalam mm). Data mencakup 1096 amatan harian dari 01 January 2022 hingga 31 December 2024 tanpa nilai hilang, sehingga seluruh data dapat langsung digunakan tanpa proses imputasi.

3 Eksplorasi Data

3.1 Statistika Deskriptif

3.1.1 Statistik Deskriptif Dasar

summary(dataset[, c("TMA", "PC")])

##       TMA              PC          
##  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

Secara rata-rata, TMA berada di sekitar 625.35 cm dengan nilai minimum 555.42 cm dan maksimum 781.25 cm, menunjukkan adanya fluktuasi tinggi muka air yang cukup besar sepanjang periode data. Curah hujan (PC) memiliki rata-rata 6.57 mm/hari, namun dengan nilai maksimum harian mencapai 145.03 mm – jauh di atas rata-ratanya – yang mengindikasikan adanya hari-hari dengan hujan ekstrem yang berpotensi memicu kenaikan TMA secara tiba-tiba.

3.1.2 Standard Deviation

sd_tma <- sd(dataset$TMA, na.rm = TRUE)
sd_pc  <- sd(dataset$PC, 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

Simpangan baku PC (8.84) relatif besar dibandingkan rata-ratanya (6.57), menegaskan bahwa curah hujan bersifat sangat bervariasi (tidak menentu) dari hari ke hari – karakteristik umum data curah hujan harian di wilayah tropis seperti Bogor.

3.1.3 Skewness

skew_tma <- skewness(dataset$TMA, na.rm = TRUE)
skew_pc  <- skewness(dataset$PC, na.rm = TRUE)

cat("Skewness TMA:", skew_tma, "\n")

## Skewness TMA: 0.5304754

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

## Skewness PC : 5.070871

Nilai skewness TMA sebesar 0.53 (menjulur ke kanan/positive skew), sedangkan PC memiliki skewness 5.071 yang menjulur kuat ke kanan. Hal ini wajar terjadi pada data curah hujan karena sebagian besar hari memiliki curah hujan rendah/tanpa hujan, namun terdapat beberapa hari dengan curah hujan sangat tinggi (nilai ekstrem) yang menarik ekor distribusi ke kanan.

3.1.4 Kurtosis

kurt_tma <- kurtosis(dataset$TMA, na.rm = TRUE)
kurt_pc  <- kurtosis(dataset$PC, na.rm = TRUE)

cat("Kurtosis TMA:", kurt_tma, "\n")

## Kurtosis TMA: 2.984404

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

## Kurtosis PC : 62.19679

Nilai kurtosis PC sebesar 62.2 (leptokurtik/lebih runcing dari distribusi normal) mengonfirmasi adanya nilai-nilai ekstrem (outlier) pada curah hujan harian, sejalan dengan hasil skewness sebelumnya.

3.2 Time Series

3.2.1 Plot Time Series TMA

ggplot(dataset, aes(x = Date, y = TMA)) +
  geom_line(color = "#00688B") +
  labs(title = "Plot Time Series Tinggi Muka Air (TMA)", x = "Tanggal", y = "TMA (cm)") +
  theme_minimal()

3.2.2 Plot Time Series PC

ggplot(dataset, aes(x = Date, y = PC)) +
  geom_line(color = "#B22222") +
  labs(title = "Plot Time Series Curah Hujan (PC)", x = "Tanggal", y = "PC (mm)") +
  theme_minimal()

3.2.3 Plot Time Series Gabungan

data_long <- dataset %>%
  pivot_longer(cols = c(TMA, PC), names_to = "Variabel", values_to = "Nilai")

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

Secara visual, pola TMA menunjukkan komponen tren/level yang bergerak naik-turun dengan cukup persisten dari waktu ke waktu (mengindikasikan kemungkinan data tidak stasioner di level), sedangkan PC berpola fluktuatif tajam antar hari (khas data curah hujan) tanpa tren yang jelas.

3.3 Matriks Korelasi

cor_tma_pc <- cor(dataset$TMA, dataset$PC, use = "complete.obs")
cor_tma_pc

## [1] 0.3945932

cor_matrix <- matrix(
  c(1, cor_tma_pc, cor_tma_pc, 1),
  nrow = 2,
  dimnames = list(c("TMA", "PC"), c("TMA", "PC"))
)

corrplot(cor_matrix, method = "color", addCoef.col = "black",
         tl.col = "black", tl.srt = 0, number.cex = 1,
         title = "Matriks Korelasi TMA dan PC",
         mar = c(0, 0, 1, 0))

Koefisien korelasi Pearson antara TMA dan PC sebesar 0.395 menunjukkan hubungan linear yang sedang dan searah (positif) secara kontemporer (pada hari yang sama). Nilai korelasi yang tidak terlalu tinggi ini wajar terjadi karena air hujan memerlukan waktu (lag) untuk terakumulasi dan memengaruhi tinggi muka air sungai/saluran air, sehingga hubungan sesungguhnya antara PC dan TMA lebih tepat dieksplorasi melalui Cross-Correlation Function (CCF) pada bagian selanjutnya, bukan hanya korelasi kontemporer.

4 Data Splitting

total_rows <- nrow(dataset)
test_size  <- 160

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

train_df <- dataset[train_index, ]
test_df  <- dataset[test_index, ]

str(train_df)

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

str(test_df)

## tibble [160 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Date: POSIXct[1:160], format: "2024-07-25" "2024-07-26" ...
##  $ TMA : num [1:160] 600 600 600 600 600 ...
##  $ PC  : num [1:160] 0.01 0.01 0 0 0 ...

train_df <- train_df %>% mutate(DataType = "Train")
test_df  <- test_df  %>% mutate(DataType = "Test")

combined_data <- bind_rows(train_df, test_df)

Data dibagi menjadi data latih (train) sebesar 85% dari total data atau sebanyak 936 amatan dan data uji (test) sebesar 15% dari total data atau sebanyak 160 amatan terakhir, mengikuti prinsip time series split di mana data uji selalu berada setelah data latih (tidak diacak) agar mencerminkan kondisi peramalan riil ke depan.

4.1 Time series data split

4.1.1 Plot Time Series TMA

data_long <- combined_data %>%
  pivot_longer(cols = c(TMA), names_to = "Variabel", values_to = "Nilai")

ggplot(data_long, aes(x = Date, y = Nilai, color = DataType)) +
  geom_line() +
  facet_wrap(~ Variabel, scales = "free_y") +
  scale_color_manual(values = c("Train" = "#00BFC4", "Test" = "#F8766D")) +
  labs(title = "Time Series TMA (Train vs. Test)", x = "Tanggal", y = "Tinggi Muka Air (cm)", color = "Jenis Data") +
  theme_minimal()

4.1.2 Plot Time Series PC

data_long <- combined_data %>%
  pivot_longer(cols = c(PC), names_to = "Variabel", values_to = "Nilai")

ggplot(data_long, aes(x = Date, y = Nilai, color = DataType)) +
  geom_line() +
  facet_wrap(~ Variabel, scales = "free_y") +
  scale_color_manual(values = c("Train" = "#00BFC4", "Test" = "#F8766D")) +
  labs(title = "Time Series PC (Train vs. Test)", x = "Tanggal", y = "Curah Hujan (mm)", color = "Jenis Data") +
  theme_minimal()

4.1.3 Plot Time Series Gabungan

data_long <- combined_data %>%
  pivot_longer(cols = c(TMA, PC), names_to = "Variabel", values_to = "Nilai")

ggplot(data_long, 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 TMA dan PC (Train vs. Test)", x = "Tanggal", y = "Nilai", color = "Jenis Data") +
  theme_minimal()

5 Kestasioneran Data

5.1 Plot ACF

5.1.1 Plot ACF TMA

acf_tma <- ggAcf(train_df$TMA) + labs(title = "ACF TMA")
acf_tma

5.1.2 Plot ACF PC

acf_pc <- ggAcf(train_df$PC) + labs(title = "ACF PC")
acf_pc

5.1.3 Plot ACF Gabungan

acf_tma + acf_pc

Plot ACF TMA dan PC meluruh secara lambat (slowly decaying), yang merupakan indikasi awal bahwa kedua data pada level tidak stasioner. Akan tetapi, ACF PC cenderung meluruh lebih cepat menuju nol, mengindikasikan PC relatif lebih stasioner dibandingkan TMA. Karakteristik ini lazim terjadi pada data hidrologi (TMA memiliki memory/persistensi tinggi, curah hujan bersifat lebih acak harian).

5.2 Plot PACF

5.2.1 Plot PACF TMA

pacf_tma <- ggPacf(train_df$TMA) + labs(title = "PACF TMA")
pacf_tma

5.2.2 Plot PACF PC

pacf_pc <- ggPacf(train_df$PC) + labs(title = "PACF PC")
pacf_pc

5.2.3 Plot PACF Gabungan

pacf_tma + pacf_pc

Berdasarkan grafik PACF, pada data TMA spike signifikan terlihat hingga sekitar lag-4 sehingga kandidat orde autoregressive (\(p\)) berada pada rentang 1-4. Sementara itu, pada data PC, spike signifikan terlihat hingga lag-3, sehingga kandidat orde autoregressive (\(p\)) berada pada rentang 1-3.

5.3 Uji ADF

5.3.1 Uji ADF TMA

adf_tma <- adf.test(train_df$TMA)
adf_tma

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_df$TMA
## Dickey-Fuller = -4.0015, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan Augmented Dickey-Fuller (ADF) test, diperoleh p-value sebesar 0.01. Jika nilai ini lebih kecil dari taraf signifikansi 5%, maka \(H_0\) ditolak sehingga data TMA dinyatakan stasioner pada tingkat level – meskipun secara visual polanya tampak persisten.

5.3.2 Uji ADF PC

adf_pc <- adf.test(train_df$PC)
adf_pc

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_df$PC
## Dickey-Fuller = -7.5589, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan Augmented Dickey-Fuller (ADF) test, diperoleh p-value sebesar 0.01, dengan interpretasi yang sama seperti TMA di atas.

5.4 Plot CCF

5.4.1 Plot CCF TMA dan PC

ccf_result <- ggCcf(train_df$TMA, train_df$PC, lag.max = 9) +
  ggtitle("CCF antara TMA dan PC")
ccf_result

5.4.2 Nilai Korelasi CCF

ccf_raw <- ccf(train_df$TMA, train_df$PC, lag.max = 9, plot = FALSE)
ccf_df <- data.frame(Lag = ccf_raw$lag[, 1, 1], Correlation = ccf_raw$acf[, 1, 1])
ccf_df

##    Lag Correlation
## 1   -9  0.12078700
## 2   -8  0.09656928
## 3   -7  0.11018666
## 4   -6  0.12779790
## 5   -5  0.12737960
## 6   -4  0.14512834
## 7   -3  0.15687309
## 8   -2  0.23595661
## 9   -1  0.24265445
## 10   0  0.39085581
## 11   1  0.46370382
## 12   2  0.35997618
## 13   3  0.34560526
## 14   4  0.33940062
## 15   5  0.30725905
## 16   6  0.26824765
## 17   7  0.28237555
## 18   8  0.26322418
## 19   9  0.25697559

lag_optimum <- ccf_df$Lag[which.max(abs(ccf_df$Correlation))]
cat("Lag dengan korelasi silang absolut tertinggi:", lag_optimum, "\n")

## Lag dengan korelasi silang absolut tertinggi: 1

Berdasarkan hasil analisis korelasi silang, diperoleh nilai korelasi silang absolut tertinggi pada lag 1. Hal ini menunjukkan bahwa perubahan curah hujan memiliki hubungan paling kuat dengan perubahan tinggi muka air pada lag tersebut, sehingga informasi ini dipakai sebagai dasar penentuan struktur lag pada pemodelan ARDL berikutnya.

6 First Difference

6.1 Differencing

full_diff <- data.frame(
  Date = dataset$Date[-1],
  DTMA = diff(dataset$TMA, differences = 1),
  DPC  = diff(dataset$PC,  differences = 1)
)

n_train <- length(train_index)
train_diff_rows <- 1:(n_train - 1)
test_diff_rows  <- n_train:(n_train + length(test_index) - 1)

df      <- full_diff[train_diff_rows, c("DTMA", "DPC")]
df_test <- full_diff[test_diff_rows,  c("DTMA", "DPC")]

DTMA <- df$DTMA
DPC  <- df$DPC
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

6.2 Plot ACF

6.2.1 Plot ACF DTMA

acf_dtma <- ggAcf(DTMA) + labs(title = "ACF DTMA")
acf_dtma

6.2.2 Plot ACF DPC

acf_dpc <- ggAcf(DPC) + labs(title = "ACF DPC")
acf_dpc

6.2.3 Plot ACF Gabungan

acf_dtma + acf_dpc

Berdasarkan grafik ACF setelah differencing, yaitu DTMA dan DPC, terlihat bahwa keduanya masih memiliki autokorelasi negatif yang signifikan pada lag-1, sedangkan sebagian besar lag berikutnya berada dalam batas signifikansi. Grafik tersebut dapat digunakan sebagai informasi tambahan mengenai pola autokorelasi.

6.3 Plot PACF

6.3.1 Plot PACF DTMA (Lag optimum peubah respon)

pacf_dtma <- ggPacf(DTMA) + labs(title = "PACF DTMA")
pacf_dtma

6.3.2 Plot PACF DPC

pacf_dpc <- ggPacf(DPC) + labs(title = "PACF DPC")
pacf_dpc

6.3.3 Plot PACF Gabungan

pacf_dtma + pacf_dpc

Berdasarkan grafik PACF, DTMA menunjukkan spike signifikan pada beberapa lag awal hingga sekitar lag-9, sedangkan DPC menunjukkan spike signifikan pada beberapa lag awal. Hasil ini mendasari pemilihan struktur lag pada model ARDL (lag DTMA 1:9, lag DPC 1) pada bagian selanjutnya.

6.4 Uji ADF

6.4.1 Uji ADF DTMA

adf_dtma <- adf.test(DTMA)
adf_dtma

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

6.4.2 Uji ADF DPC

adf_dpc <- adf.test(DPC)
adf_dpc

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

Setelah differencing satu kali, p-value ADF untuk DTMA sebesar 0.01 dan untuk DPC sebesar 0.01, mengonfirmasi bahwa kedua peubah stasioner setelah pembedaan pertama – konsisten dengan hasil uji ADF pada level yang juga sudah menyatakan stasioner, sehingga differencing di sini berfungsi sebagai analisis tambahan, bukan syarat mutlak.

6.5 Plot CCF (Lag optimum peubah penjelas)

6.5.1 Plot CCF DTMA dan DPC

ccf_d_result <- ggCcf(DTMA, DPC, lag.max = 9) +
  ggtitle("CCF antara DTMA dan DPC")
ccf_d_result

6.5.2 Nilai Korelasi CCF

ccf_d_raw <- ccf(DTMA, DPC, lag.max = 9, plot = FALSE)
ccf_d_df <- data.frame(Lag = ccf_d_raw$lag[, 1, 1], Correlation = ccf_d_raw$acf[, 1, 1])
ccf_d_df

##    Lag  Correlation
## 1   -9  0.052958627
## 2   -8 -0.055055240
## 3   -7 -0.007564919
## 4   -6  0.027277549
## 5   -5 -0.026068771
## 6   -4  0.009520092
## 7   -3 -0.098523738
## 8   -2  0.105694759
## 9   -1 -0.206857132
## 10   0  0.108686920
## 11   1  0.258098583
## 12   2 -0.130903673
## 13   3 -0.012055675
## 14   4  0.037902380
## 15   5  0.010291092
## 16   6 -0.077987191
## 17   7  0.048995088
## 18   8 -0.019289336
## 19   9  0.008364891

lag_optimum_d <- ccf_d_df$Lag[which.max(abs(ccf_d_df$Correlation))]
cat("Lag optimum (data differencing):", lag_optimum_d, "\n")

## Lag optimum (data differencing): 1

Nilai korelasi silang pada data hasil differencing memperkuat hasil sebelumnya, dengan lag ke-1 menunjukkan hubungan paling kuat antara perubahan curah hujan (DPC) dengan perubahan tinggi muka air (DTMA).

7 Pemodelan ARDL (dynardl/dynamac)

7.1 Uji Kointegrasi

Uji kointegrasi dilakukan dengan pendekatan ARDL Bound Test (Pesaran, Shin, & Smith) untuk menguji ada tidaknya hubungan jangka panjang (kointegrasi) antara TMA dan PC, tanpa mensyaratkan kedua peubah berordo integrasi sama.

modelbound <- ardlBound(
  data = train_df,
  formula = TMA ~ PC,
  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.20482, df1 = 1, df2 = 911, p-value = 0.651
## 
## ------------------------------------------------------ 
## 
##  Ljung-Box Test for the autocorrelation in residuals:
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 0.0026547, df = 1, p-value = 0.9589
## 
## ------------------------------------------------------ 
## 
##  Breusch-Pagan Test for the homoskedasticity of residuals:
## 
##  studentized Breusch-Pagan test
## 
## data:  modelFull$model
## BP = 201.72, df = 13, p-value < 2.2e-16
## 
## The p-value of Breusch-Pagan test for the homoskedasticity of residuals:  6.053769e-36 < 0.05!
## ------------------------------------------------------ 
## 
##  Shapiro-Wilk test of normality of residuals:
## 
##  Shapiro-Wilk normality test
## 
## data:  modelFull$model$residual
## W = 0.8827, p-value < 2.2e-16
## 
## The p-value of Shapiro-Wilk test normality of residuals:  7.210048e-26 < 0.05!
## ------------------------------------------------------ 
## 
##  PESARAN, SHIN AND SMITH (2001) COINTEGRATION TEST 
## 
##  Observations: 935 
##  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.4329166719859 
##   
##  ------------------------------------------------------ 
##  F-statistic note: Asymptotic critical values used. 
##  
## ------------------------------------------------------ 
## 
##  Ramsey's RESET Test for model specification:
## 
##  RESET test
## 
## data:  modelECM$model
## RESET = 18.987, df1 = 1, df2 = 913, p-value = 1.466e-05
## 
## the p-value of RESET test:  1.465705e-05 < 0.05!
## ------------------------------------------------------

## ------------------------------------------------------ 
## Error Correction Model Output: 
## 
## Time series regression with "ts" data:
## Start = 10, End = 935
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -131.870  -11.401   -2.585    6.663  135.093 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## ec.1   -0.15599    0.01398 -11.161  < 2e-16 ***
## dPC.t   0.75402    0.08992   8.386  < 2e-16 ***
## dPC.1   0.33707    0.10054   3.353 0.000833 ***
## dTMA.1 -0.43910    0.03012 -14.577  < 2e-16 ***
## dTMA.2 -0.28628    0.03345  -8.557  < 2e-16 ***
## dTMA.3 -0.21494    0.03413  -6.299 4.66e-10 ***
## dTMA.4 -0.13859    0.03459  -4.006 6.67e-05 ***
## dTMA.5 -0.15805    0.03409  -4.636 4.07e-06 ***
## dTMA.6 -0.12899    0.03407  -3.786 0.000163 ***
## dTMA.7 -0.10075    0.03369  -2.990 0.002861 ** 
## dTMA.8 -0.04581    0.03262  -1.404 0.160526    
## dTMA.9  0.02087    0.02945   0.709 0.478674    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.31 on 914 degrees of freedom
## Multiple R-squared:  0.3656, Adjusted R-squared:  0.3573 
## F-statistic: 43.89 on 12 and 914 DF,  p-value: < 2.2e-16
## 
## ------------------------------------------------------ 
## Long-run coefficients: 
##      TMA.1       PC.1 
## -0.1559866  1.2969018 
##

f_stat_bound <- 41.4329166719859   # sesuai "F-statistic =" 
i1_5persen   <- 4.16               # sesuai baris "5% critical value" kolom I(1)

Statistik-F yang diperoleh sebesar 41.433, jauh lebih besar dibandingkan batas atas I(1) pada taraf nyata 5% sebesar 4.16. Karena statistik-F melampaui batas atas I(1), maka \(H_0\) (tidak ada hubungan jangka panjang) ditolak, sehingga disimpulkan terdapat hubungan jangka panjang (kointegrasi) antara TMA dan PC.

7.2 Model

7.2.1 Jangka Pendek

Model dinamis jangka pendek (error correction form) diduga menggunakan fungsi dynardl() dari paket dynamac, dengan lag autoregresif DTMA hingga lag ke-9 (mengikuti pola PACF DTMA) dan lag DPC sebesar 1 (mengikuti hasil CCF data differencing).

modelardl <- dynardl(
  DTMA ~ DPC,
  data = df,
  lags = list(
    "DTMA" = 1:9,
    "DPC"  = 1
  ),
  ec = TRUE
)

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

summary(modelardl)

## 
## Call:
## lm(formula = as.formula(paste(paste(dvnamelist), "~", paste(colnames(IVs), 
##     collapse = "+"), collapse = " ")))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -89.431 -12.296  -3.093   6.997 141.229 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.07623    0.75510   0.101 0.919613    
## l.1.DTMA    -1.48033    0.03166 -46.754  < 2e-16 ***
## l.2.DTMA    -0.25904    0.03557  -7.284 7.01e-13 ***
## l.3.DTMA    -0.20871    0.03622  -5.763 1.13e-08 ***
## l.4.DTMA    -0.12523    0.03686  -3.398 0.000709 ***
## l.5.DTMA    -0.16258    0.03655  -4.448 9.73e-06 ***
## l.6.DTMA    -0.12678    0.03663  -3.461 0.000563 ***
## l.7.DTMA    -0.10015    0.03622  -2.765 0.005813 ** 
## l.8.DTMA    -0.05345    0.03499  -1.528 0.126901    
## l.9.DTMA     0.02171    0.03175   0.684 0.494307    
## l.1.DPC      0.76884    0.08443   9.106  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.98 on 915 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.7247, Adjusted R-squared:  0.7217 
## F-statistic: 240.8 on 10 and 915 DF,  p-value: < 2.2e-16

AIC(modelardl$model)

## [1] 8445.934

coefs <- summary(modelardl)$coefficients
b1_dpc <- coefs["l.1.DPC", "Estimate"]
p_dpc  <- coefs["l.1.DPC", "Pr(>|t|)"]

Koefisien l.1.DPC sebesar 0.76884 dengan p-value <2e-16. Koefisien yang bernilai positif menunjukkan bahwa perubahan PC satu periode sebelumnya berhubungan positif dengan perubahan TMA pada periode sekarang.

7.2.2 Jangka Panjang

Model jangka panjang diduga menggunakan level (bukan bentuk differencing) untuk menangkap hubungan keseimbangan (equilibrium) antara TMA dan PC dalam jangka panjang.

modelardll <- dynardl(
  formula = TMA ~ PC,
  data = train_df,
  levels = c("PC"),
  ec = FALSE
)

## [1] "Dependent variable to be run in levels."

summary(modelardll)

## 
## Call:
## lm(formula = as.formula(paste(paste(dvnamelist), "~", paste(colnames(IVs), 
##     collapse = "+"), collapse = " ")))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -91.775 -10.860  -3.639   9.701 135.091 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 155.38860   12.04126   12.90   <2e-16 ***
## l.1.TMA       0.74133    0.01940   38.21   <2e-16 ***
## PC            0.99491    0.09169   10.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.65 on 932 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6698, Adjusted R-squared:  0.6691 
## F-statistic: 945.2 on 2 and 932 DF,  p-value: < 2.2e-16

AIC(modelardll$model)

## [1] 8573.629

Model jangka pendek menjelaskan dinamika perubahan TMA dari waktu ke waktu, sedangkan model jangka panjang menggambarkan hubungan keseimbangan antara level TMA dan PC. Kedua model ini tidak dapat dibandingkan langsung menggunakan AIC karena peubah dependennya berbeda (DTMA vs TMA) – pemilihan model perlu disesuaikan dengan tujuan analisis, bukan berdasarkan kriteria informasi.

coefs_ll   <- summary(modelardll)$coefficients
b_PC       <- coefs_ll["PC", "Estimate"]
b_TMA_lag1 <- coefs_ll["l.1.TMA", "Estimate"]
longrun_PC_level_model <- b_PC / (1 - b_TMA_lag1)

cat("Long-run coefficient PC (dari model level modelardll):", longrun_PC_level_model, "\n")

## Long-run coefficient PC (dari model level modelardll): 3.846219

Berdasarkan hasil estimasi model level (modelardll), diperoleh koefisien PC kontemporer sebesar 0.99491 dan koefisien l.1.TMA sebesar 0.74133. Koefisien jangka panjang PC dari model ini dihitung sebagai koefisien PC dibagi (1 - koefisien l.1.TMA), yaitu sebesar 3.8462.

7.3 Uji Asumsi

7.3.1 Uji Autokorelasi Residual (Box-Ljung Test)

res1 <- residuals(modelardl$model)
lb_test <- Box.test(res1, lag = 20, type = "Ljung-Box")
lb_test

## 
##  Box-Ljung test
## 
## data:  res1
## X-squared = 19.245, df = 20, p-value = 0.506

Dengan p-value sebesar 0.506, residual model tidak menunjukkan autokorelasi yang signifikan (gagal tolak H0), sehingga asumsi non-autokorelasi terpenuhi.

7.3.2 Uji Normalitas (Jarque-Bera Test)

jb_test <- jarque.bera.test(res1)
jb_test

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

P-value uji Jarque-Bera sebesar <2e-16 mengindikasikan bahwa residual model tidak menyebar normal (tolak H0), meskipun pada data berukuran besar pelanggaran normalitas umumnya tidak terlalu memengaruhi validitas estimasi berkat sifat asimtotik OLS.

7.3.3 Uji Heteroskedastisitas (Breusch-Pagan Test)

fitted_val1 <- fitted(modelardl$model)
bp_test <- bptest(res1 ~ fitted_val1)
bp_test

## 
##  studentized Breusch-Pagan test
## 
## data:  res1 ~ fitted_val1
## BP = 1.377, df = 1, p-value = 0.2406

P-value uji Breusch-Pagan sebesar 0.2406 menunjukkan bahwa ragam residual bersifat homogen (homoskedastis), sehingga asumsi kesamaan ragam terpenuhi.

7.4 Forecasting

coefs <- summary(modelardl)$coefficients
intercept <- coefs["(Intercept)", "Estimate"]
phi       <- coefs[paste0("l.", 1:9, ".DTMA"), "Estimate"]
beta_dpc  <- coefs["l.1.DPC", "Estimate"]

7.4.0.1 Pemeriksaan Kestabilan Model AR(9)

p_ar <- length(phi)
comp_matrix <- matrix(0, p_ar, p_ar)
comp_matrix[1, ] <- phi
if (p_ar > 1) {
  for (i in 2:p_ar) comp_matrix[i, i - 1] <- 1
}
akar_karakteristik <- eigen(comp_matrix, only.values = TRUE)$values
modulus_akar <- Mod(akar_karakteristik)

data.frame(akar = akar_karakteristik, modulus = round(modulus_akar, 4))

##                     akar modulus
## 1 -1.38132343+0.0000000i  1.3813
## 2  0.55701157+0.4175298i  0.6961
## 3  0.55701157-0.4175298i  0.6961
## 4 -0.43806137+0.5200372i  0.6800
## 5 -0.43806137-0.5200372i  0.6800
## 6 -0.66625350+0.0000000i  0.6663
## 7  0.04232685+0.6545533i  0.6559
## 8  0.04232685-0.6545533i  0.6559
## 9  0.24468799+0.0000000i  0.2447

cat("Modulus akar terbesar:", round(max(modulus_akar), 4), "\n")

## Modulus akar terbesar: 1.3813

if (max(modulus_akar) > 1) {
  cat("Modulus akar terbesar >1, artinya struktur AR(9) pada model ini **tidak stabil secara dinamis (explosive)**. Konsekuensinya, forecasting rekursif murni (memakai hasil prediksi sebelumnya sebagai lag untuk memprediksi langkah berikutnya, berulang-ulang) akan menghasilkan nilai yang membesar tanpa batas seiring bertambahnya horizon peramalan, meskipun model ini menunjukkan kecocokan (fit) yang baik untuk prediksi satu langkah ke depan menggunakan nilai lag aktual.\n\nKarena TMA pada praktiknya diukur/dicatat setiap hari, forecasting pada data latih, data uji, maupun prediksi 2025 di bawah ini dilakukan dengan pendekatan **rolling one-step-ahead**: setiap prediksi TMA hari ke-$t$ dihitung menggunakan nilai TMA/DTMA **aktual** 9 hari sebelumnya (bukan hasil prediksi hari-hari sebelumnya), sehingga kesalahan tidak terakumulasi dari hari ke hari. Pendekatan ini merupakan evaluasi yang realistis untuk model peramalan harian, meskipun tidak mencerminkan simulasi dinamis murni beberapa hari ke depan tanpa pembaruan data.\n\n")
} else {
  cat("Modulus seluruh akar karakteristik <= 1, sehingga struktur AR(9) pada model ini stabil secara dinamis, dan forecasting rekursif murni dapat digunakan tanpa risiko nilai yang membesar tanpa batas.\n\n")
}

Modulus akar terbesar >1, artinya struktur AR(9) pada model ini tidak stabil secara dinamis (explosive). Konsekuensinya, forecasting rekursif murni (memakai hasil prediksi sebelumnya sebagai lag untuk memprediksi langkah berikutnya, berulang-ulang) akan menghasilkan nilai yang membesar tanpa batas seiring bertambahnya horizon peramalan, meskipun model ini menunjukkan kecocokan (fit) yang baik untuk prediksi satu langkah ke depan menggunakan nilai lag aktual.

Karena TMA pada praktiknya diukur/dicatat setiap hari, forecasting pada data latih, data uji, maupun prediksi 2025 di bawah ini dilakukan dengan pendekatan rolling one-step-ahead: setiap prediksi TMA hari ke-\(t\) dihitung menggunakan nilai TMA/DTMA aktual 9 hari sebelumnya (bukan hasil prediksi hari-hari sebelumnya), sehingga kesalahan tidak terakumulasi dari hari ke hari. Pendekatan ini merupakan evaluasi yang realistis untuk model peramalan harian, meskipun tidak mencerminkan simulasi dinamis murni beberapa hari ke depan tanpa pembaruan data.

forecast_dtma_one_step <- function(idx_target, dtma_series, dpc_series,
                                    intercept, phi, beta_dpc) {
  lag_vals <- rev(dtma_series[(idx_target - 9):(idx_target - 1)])
  dpc_l1   <- dpc_series[idx_target - 1]
  intercept + sum(phi * lag_vals) + beta_dpc * dpc_l1
}

7.4.1 Data Training

DTMA_fitted_train <- fitted(modelardl$model)   # 9 obs pertama hilang karena lag

TMA_actual_lag  <- train_df$TMA[(1 + 9):(nrow(train_df) - 1)]   # TMA aktual t-1
TMA_actual_train <- train_df$TMA[(1 + 9 + 1):nrow(train_df)]    # TMA aktual t
TMA_pred_train   <- TMA_actual_lag + DTMA_fitted_train           # TMA_pred[t] = TMA_aktual[t-1] + DTMA_fitted[t]

plot(TMA_actual_train, type = "l", col = "black",
     main = "Data Latih: Aktual vs Prediksi TMA (one-step-ahead, anchor ke TMA aktual)")
lines(TMA_pred_train, col = "blue")
legend("topleft", legend = c("Aktual", "Prediksi"), col = c("black", "blue"), lty = 1)

MAE_train  <- mean(abs(TMA_pred_train - TMA_actual_train), na.rm = TRUE)
RMSE_train <- sqrt(mean((TMA_pred_train - TMA_actual_train)^2, na.rm = TRUE))
MAPE_train <- mean(abs((TMA_actual_train - TMA_pred_train) / TMA_actual_train), na.rm = TRUE) * 100

cat("MAE :", MAE_train, "\n")

## MAE : 23.91399

cat("RMSE:", RMSE_train, "\n")

## RMSE: 35.04561

cat("MAPE:", MAPE_train, "%\n")

## MAPE: 3.723197 %

Pada data latih, model menghasilkan MAPE sebesar 3.72% untuk prediksi satu hari ke depan menggunakan TMA aktual hari sebelumnya sebagai acuan.

7.4.2 Data testing

idx_test_start <- n_train

DTMA_pred_test <- sapply(seq_along(df_test$DPC), function(t) {
  forecast_dtma_one_step(
    idx_target  = idx_test_start + t - 1,
    dtma_series = full_diff$DTMA,
    dpc_series  = full_diff$DPC,
    intercept = intercept, phi = phi, beta_dpc = beta_dpc
  )
})

TMA_actual_lag_test <- c(tail(train_df$TMA, 1), test_df$TMA[-nrow(test_df)])   # TMA aktual t-1 untuk tiap target di data uji
TMA_pred_test  <- TMA_actual_lag_test + DTMA_pred_test
TMA_actual_test <- test_df$TMA
date_seq        <- test_df$Date

plot(date_seq, TMA_actual_test, type = "l", col = "black",
     main = "Data Uji: Aktual vs Prediksi TMA (rolling one-step-ahead)",
     xlab = "Tanggal", ylab = "TMA (cm)")
lines(date_seq, TMA_pred_test, col = "blue")
legend("topleft", legend = c("Aktual", "Prediksi"), col = c("black", "blue"), lty = 1)

MAE_test  <- mean(abs(TMA_pred_test - TMA_actual_test), na.rm = TRUE)
RMSE_test <- sqrt(mean((TMA_pred_test - TMA_actual_test)^2, na.rm = TRUE))
MAPE_test <- mean(abs((TMA_actual_test - TMA_pred_test) / TMA_actual_test), na.rm = TRUE) * 100

cat("MAE :", MAE_test, "\n")

## MAE : 16.09976

cat("RMSE:", RMSE_test, "\n")

## RMSE: 26.43804

cat("MAPE:", MAPE_test, "%\n")

## MAPE: 2.528116 %

Pada data uji, dengan pendekatan rolling one-step-ahead (memakai TMA aktual hari sebelumnya sebagai acuan setiap hari, bukan hasil prediksi berantai), model menghasilkan MAPE sebesar 2.53%.

7.4.3 Prediksi 2025

Selain dievaluasi pada data uji internal, model juga diuji pada data benar-benar baru yang belum pernah dilihat sama sekali oleh model saat pendugaan, yaitu data harian periode Januari-Maret 2025.

data_prediksi_25 <- read_excel("C:/Users/hp/Downloads/Data_Jan25_Mar25_mpdw.xlsx")
data_prediksi_25$Date <- as.Date(data_prediksi_25$Date)
data_prediksi_25 <- data_prediksi_25 %>% arrange(Date)
str(data_prediksi_25)

## tibble [90 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Date: Date[1:90], format: "2025-01-01" "2025-01-02" ...
##  $ TMA : num [1:90] 641 640 640 640 640 ...
##  $ PC  : num [1:90] 5.27 2.45 6.22 6.03 2.63 ...

last_row_2024 <- dataset[nrow(dataset), c("Date", "TMA", "PC")]
data_25_ext   <- bind_rows(last_row_2024, data_prediksi_25[, c("Date", "TMA", "PC")])

full_diff_25 <- data.frame(
  Date = data_25_ext$Date[-1],
  DTMA = diff(data_25_ext$TMA, differences = 1),
  DPC  = diff(data_25_ext$PC,  differences = 1)
)

combined_diff_25 <- data.frame(
  Date = c(full_diff$Date, full_diff_25$Date),
  DTMA = c(full_diff$DTMA, full_diff_25$DTMA),
  DPC  = c(full_diff$DPC,  full_diff_25$DPC)
)
idx_25_start <- nrow(full_diff) + 1   
DTMA_pred_25 <- sapply(seq_len(nrow(full_diff_25)), function(t) {
  forecast_dtma_one_step(
    idx_target  = idx_25_start + t - 1,
    dtma_series = combined_diff_25$DTMA,
    dpc_series  = combined_diff_25$DPC,
    intercept = intercept, phi = phi, beta_dpc = beta_dpc
  )
})

TMA_actual_lag_25 <- c(last_row_2024$TMA, data_prediksi_25$TMA[-nrow(data_prediksi_25)])
TMA_pred_25 <- TMA_actual_lag_25 + DTMA_pred_25

hasil_prediksi <- data.frame(
  date       = data_prediksi_25$Date,
  TMA_actual = data_prediksi_25$TMA,
  TMA_pred   = TMA_pred_25
)
hasil_prediksi

##          date TMA_actual TMA_pred
## 1  2025-01-01   640.8000 637.6159
## 2  2025-01-02   640.0000 642.2061
## 3  2025-01-03   640.0000 638.8851
## 4  2025-01-04   640.0000 643.0150
## 5  2025-01-05   640.0000 639.9969
## 6  2025-01-06   640.0000 637.4323
## 7  2025-01-07   640.0000 656.9963
## 8  2025-01-08   641.4583 629.7566
## 9  2025-01-09   666.2500 649.1697
## 10 2025-01-10   642.5000 626.1948
## 11 2025-01-11   637.5000 660.8647
## 12 2025-01-12   634.5833 651.1169
## 13 2025-01-13   636.6667 634.0375
## 14 2025-01-14   636.6667 637.3765
## 15 2025-01-15   647.2917 650.9729
## 16 2025-01-16   642.9167 619.5552
## 17 2025-01-17   640.0000 645.4124
## 18 2025-01-18   640.0000 660.2725
## 19 2025-01-19   640.4167 631.5599
## 20 2025-01-20   645.0000 643.8032
## 21 2025-01-21   642.0833 644.8423
## 22 2025-01-22   645.8333 644.7517
## 23 2025-01-23   640.0000 627.9360
## 24 2025-01-24   640.0000 665.2013
## 25 2025-01-25   642.5000 636.5349
## 26 2025-01-26   630.0000 626.3179
## 27 2025-01-27   630.0000 663.4112
## 28 2025-01-28   646.2500 612.8474
## 29 2025-01-29   692.2917 646.0369
## 30 2025-01-30   715.4167 646.4065
## 31 2025-01-31   666.6304 629.6007
## 32 2025-02-01   652.1739 716.9412
## 33 2025-02-02   650.8333 676.7328
## 34 2025-02-03   648.7500 651.2497
## 35 2025-02-04   640.0000 649.0273
## 36 2025-02-05   640.0000 656.4172
## 37 2025-02-06   640.0000 649.5095
## 38 2025-02-07   643.3333 650.3855
## 39 2025-02-08   655.0000 644.8293
## 40 2025-02-09   635.8333 638.1617
## 41 2025-02-10   634.5833 665.8154
## 42 2025-02-11   630.0000 631.9928
## 43 2025-02-12   630.0000 637.5518
## 44 2025-02-13   630.0000 630.7541
## 45 2025-02-14   630.0000 633.2076
## 46 2025-02-15   630.0000 634.0666
## 47 2025-02-16   630.0000 632.6095
## 48 2025-02-17   630.0000 641.8015
## 49 2025-02-18   637.0833 632.6463
## 50 2025-02-19   640.0000 623.0321
## 51 2025-02-20   642.2917 639.8135
## 52 2025-02-21   642.5000 622.0106
## 53 2025-02-22   640.0000 641.8083
## 54 2025-02-23   640.0000 640.6439
## 55 2025-02-24   640.0000 647.5168
## 56 2025-02-25   669.7917 638.3898
## 57 2025-02-26   641.2500 619.7407
## 58 2025-02-27   645.4167 675.2089
## 59 2025-02-28   676.7391 646.6346
## 60 2025-03-01   707.9167 636.7807
## 61 2025-03-02   652.5000 647.1427
## 62 2025-03-03   746.4583 734.4395
## 63 2025-03-04   796.0417 625.5365
## 64 2025-03-05   706.2500 670.8162
## 65 2025-03-06   670.0000 798.6682
## 66 2025-03-07   681.6667 727.7606
## 67 2025-03-08   658.7500 680.7333
## 68 2025-03-09   683.9583 684.9245
## 69 2025-03-10   672.9167 656.9007
## 70 2025-03-11   660.0000 691.7401
## 71 2025-03-12   660.0000 689.6858
## 72 2025-03-13   660.0000 675.2839
## 73 2025-03-14   660.0000 664.1578
## 74 2025-03-15   660.0000 667.4359
## 75 2025-03-16   660.0000 660.2925
## 76 2025-03-17   660.8333 654.7164
## 77 2025-03-18   730.8333 667.8425
## 78 2025-03-19   665.0000 624.0074
## 79 2025-03-20   660.0000 743.2292
## 80 2025-03-21   660.0000 667.6190
## 81 2025-03-22   660.0000 665.7944
## 82 2025-03-23   660.0000 664.0360
## 83 2025-03-24   660.0000 659.0338
## 84 2025-03-25   653.7500 660.9430
## 85 2025-03-26   650.0000 665.7358
## 86 2025-03-27   651.2500 663.4623
## 87 2025-03-28   654.1667 650.8821
## 88 2025-03-29   650.0000 652.9496
## 89 2025-03-30   650.0000 653.7997
## 90 2025-03-31   650.2174 652.6227

nama_file_excel4 <- "hasil_peramalan_2025.xlsx"
write_xlsx(hasil_prediksi, path = nama_file_excel4)

ggplot(hasil_prediksi, aes(x = date, y = TMA_pred)) +
  geom_line(color = "blue", linewidth = 1) +
  labs(title = "Prediksi Tinggi Muka Air (TMA) Januari-Maret 2025",
       x = "Tanggal", y = "TMA Prediksi (cm)") +
  theme_minimal()

plot(hasil_prediksi$date, hasil_prediksi$TMA_actual, type = "l", col = "black",
     ylab = "TMA (cm)", xlab = "Tanggal",
     main = "TMA Aktual vs Prediksi Januari-Maret 2025", cex.axis = 1.2)
lines(hasil_prediksi$date, hasil_prediksi$TMA_pred, col = "blue")
legend("topleft", legend = c("Aktual", "Prediksi"), col = c("black", "blue"), lty = 1)

MAE_25  <- mean(abs(hasil_prediksi$TMA_actual - hasil_prediksi$TMA_pred))
RMSE_25 <- sqrt(mean((hasil_prediksi$TMA_actual - hasil_prediksi$TMA_pred)^2))
MAPE_25 <- mean(abs((hasil_prediksi$TMA_actual - hasil_prediksi$TMA_pred) / hasil_prediksi$TMA_actual)) * 100

cat("MAE :", MAE_25, "\n")

## MAE : 18.71324

cat("RMSE:", RMSE_25, "\n")

## RMSE: 32.55252

cat("MAPE:", MAPE_25, "%\n")

## MAPE: 2.769723 %

Pada data benar-benar baru periode Januari-Maret 2025, dengan pendekatan rolling one-step-ahead yang sama, model menghasilkan MAPE sebesar 2.77%.

Keterbatasan:

Model AR(9) pada DTMA terbukti tidak stabil secara dinamis (lihat pemeriksaan akar karakteristik di atas), sehingga forecasting genuinely multi-hari ke depan tanpa pembaruan data (simulasi dinamis murni) tidak dapat diandalkan dengan spesifikasi lag saat ini. Evaluasi di atas menggunakan pendekatan rolling one-step-ahead sebagai alternatif yang realistis, dengan asumsi TMA diperbarui/diukur setiap hari.

Forecasting pada bagian data uji dan prediksi 2025 di atas masih memakai nilai curah hujan (PC) aktual sebagai input, bukan curah hujan yang diramalkan – konsisten dengan asumsi TMA yang juga memakai data aktual hari sebelumnya pada pendekatan rolling one-step-ahead ini.

Analisis ARDL: Pengaruh Curah Hujan terhadap Tinggi Muka Air Kota Bogor

Fadly Syabani

09 September 2026