Peubah yang digunakan: - TMA = Tinggi Muka Air (peubah respon) - RH = Kelembapan relatif Kota Bogor (peubah penjelas)

1 Library

library(readxl)
## Warning: package 'readxl' was built under R version 4.5.2
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
library(moments)
## Warning: package 'moments' was built under R version 4.5.2
library(ggcorrplot)
## Warning: package 'ggcorrplot' was built under R version 4.5.3
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.5.2
## corrplot 0.95 loaded
library(forecast)
## Warning: package 'forecast' was built under R version 4.5.3
library(urca)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(dynlm)
## Warning: package 'dynlm' was built under R version 4.5.3
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(dLagM)
## Warning: package 'dLagM' was built under R version 4.5.3
## Loading required package: nardl
## Warning: package 'nardl' was built under R version 4.5.3
## 
## Attaching package: 'dLagM'
## The following object is masked from 'package:forecast':
## 
##     forecast
library(lmtest)
library(car)
## Warning: package 'car' was built under R version 4.5.2
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(ARDL)
## Warning: package 'ARDL' was built under R version 4.5.3
## To cite the ARDL package in publications:
## 
## Use this reference to refer to the validity of the ARDL package.
## 
##   Natsiopoulos, Kleanthis, and Tzeremes, Nickolaos G. (2022). ARDL
##   bounds test for cointegration: Replicating the Pesaran et al. (2001)
##   results for the UK earnings equation using R. Journal of Applied
##   Econometrics, 37(5), 1079-1090. https://doi.org/10.1002/jae.2919
## 
## Use this reference to cite this specific version of the ARDL package.
## 
##   Kleanthis Natsiopoulos, Nickolaos Tzeremes and Daniel Finnan (2026).
##   ARDL: ARDL, ECM and Bounds-Test for Cointegration. R package version
##   0.2.5. https://CRAN.R-project.org/package=ARDL
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.5.3
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.5.2
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
## 
##     area
library(dynlm)
library(nardl)
library(writexl)
## Warning: package 'writexl' was built under R version 4.5.2
library(openxlsx2)
## Warning: package 'openxlsx2' was built under R version 4.5.3
## 
## Attaching package: 'openxlsx2'
## The following object is masked from 'package:writexl':
## 
##     write_xlsx
## The following object is masked from 'package:readxl':
## 
##     read_xlsx
library(dynamac)
## Warning: package 'dynamac' was built under R version 4.5.3
## 
## Attaching package: 'dynamac'
## The following object is masked from 'package:nardl':
## 
##     pssbounds
library(e1071)
## Warning: package 'e1071' was built under R version 4.5.3
## 
## Attaching package: 'e1071'
## The following objects are masked from 'package:moments':
## 
##     kurtosis, moment, skewness
## The following object is masked from 'package:ggplot2':
## 
##     element

2 Import Data

2.1 Data

dataset <- read_excel("C:\\Users\\ASUS\\Downloads\\Data_Harian_2022_2024.xlsx")
str(dataset)
## tibble [1,096 × 6] (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 ...
##  $ TM  : num [1:1096] 24.1 24.2 24.8 24.9 25.5 ...
##  $ PC  : num [1:1096] 0.85 1.79 0.4 0.22 5.7 ...
##  $ RH  : num [1:1096] 86.4 88.1 83.2 80.3 83.2 ...
##  $ WS  : num [1:1096] 2.08 2.25 1.56 1.11 0.79 1.38 1.92 1.35 1.82 1.41 ...

3 Eksplorasi Data

3.1 Statistika Deskriptif

3.1.1 Statistik Deskriptif Dasar

summary(dataset)
##       Date                          TMA              TM       
##  Min.   :2022-01-01 00:00:00   Min.   :555.4   Min.   :22.85  
##  1st Qu.:2022-10-01 18:00:00   1st Qu.:597.7   1st Qu.:24.57  
##  Median :2023-07-02 12:00:00   Median :620.0   Median :25.10  
##  Mean   :2023-07-02 12:00:00   Mean   :625.3   Mean   :25.18  
##  3rd Qu.:2024-04-01 06:00:00   3rd Qu.:650.1   3rd Qu.:25.78  
##  Max.   :2024-12-31 00:00:00   Max.   :781.2   Max.   :27.97  
##        PC                 RH              WS       
##  Min.   :  0.0000   Min.   :67.03   Min.   :0.460  
##  1st Qu.:  0.7375   1st Qu.:82.92   1st Qu.:1.010  
##  Median :  4.0000   Median :86.29   Median :1.230  
##  Mean   :  6.5700   Mean   :85.29   Mean   :1.377  
##  3rd Qu.:  9.1325   3rd Qu.:88.53   3rd Qu.:1.583  
##  Max.   :145.0300   Max.   :94.31   Max.   :3.700

3.1.2 Standar Deviasi

sd_rh <- sd(dataset$RH, na.rm = TRUE)
sd_rh
## [1] 4.710888

3.1.3 Skewness

skew_rh <- skewness(dataset$RH, na.rm = TRUE)
skew_rh
## [1] -0.986807

3.1.4 Kurtosis

kurt_rh <- kurtosis(dataset$RH, na.rm = TRUE)
kurt_rh
## [1] 1.017691

3.2 Time Series

ggplot(data = dataset, aes(x = Date, y = RH)) +
  geom_line() +
  labs(title = "Time Series Plot of RH", x = "Date", y = "RH") +
  theme_minimal()

4 Data Splitting

total_rows <- 1096
test_size <- 160

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

train_2225 <- dataset[train_index, ]
test_2225  <- dataset[test_index, ]
str(train_2225)
## tibble [936 × 6] (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 ...
##  $ TM  : num [1:936] 24.1 24.2 24.8 24.9 25.5 ...
##  $ PC  : num [1:936] 0.85 1.79 0.4 0.22 5.7 ...
##  $ RH  : num [1:936] 86.4 88.1 83.2 80.3 83.2 ...
##  $ WS  : num [1:936] 2.08 2.25 1.56 1.11 0.79 1.38 1.92 1.35 1.82 1.41 ...

Data dibagi menjadi data train sebesar 86% dari total data (936 data) dan data test sebesar 14% (160 data).

str(test_2225)
## tibble [160 × 6] (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 ...
##  $ TM  : num [1:160] 24.4 24.9 24.9 24.7 24.5 ...
##  $ PC  : num [1:160] 0.01 0.01 0 0 0 ...
##  $ RH  : num [1:160] 85.3 84.4 82.3 82.3 81.6 ...
##  $ WS  : num [1:160] 1.32 1.13 1.24 1.01 1.16 1.35 1.25 1.01 0.98 0.99 ...
train_2225.ts<-ts(train_2225)
test_2225.ts<-ts(test_2225)

4.1 Time Series Data Split

train_2225 <- train_2225 %>%
  mutate(DataType = "Train")

test_2225 <- test_2225 %>%
  mutate(DataType = "Test")
combined_data <- bind_rows(train_2225, test_2225)
data_long <- combined_data %>%
  pivot_longer(cols = c(RH), 
               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 = 2) +
  scale_color_manual(values = c("Train" = "#00BFC4", "Test" = "#F8766D")) +
  labs(title = "Time Series of RH (Train vs. Test)",
       x = "Date",
       y = "Kelembapan Relatif (persen)",
       color = "Data Type") +
  theme_minimal()

5 Kestasioneran Data

5.1 Plot ACF

acf_RH <- ggAcf(train_2225$RH) + labs(title = "ACF RH")
acf_RH

plot ACF yang menurun secara perlahan mengindikasikan data tidak stasioner.

acf_RH <- ggAcf(train_2225$RH) +
  labs(title = "ACF RH") +
  theme(
    axis.text.x = element_text(size = 14),
    axis.text.y = element_text(size = 14)
  )
acf_RH

pacf_RH <- ggPacf(train_2225$RH) + labs(title = "PACF RH")
pacf_RH

5.3 Uji ADF

adf.test(train_2225$RH)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_2225$RH
## Dickey-Fuller = -3.3561, Lag order = 9, p-value = 0.06101
## alternative hypothesis: stationary

5.4 Plot CCF

ccf_result3 <- ggCcf(train_2225$TMA, train_2225$RH, lag.max=9, plot = TRUE)+ggtitle("CCF antara TMA dan RH")
ccf_result3

6 Firs Difference

6.1 Differencing

DTMA <- diff(train_2225$TMA, differences = 1)
DPC <- diff(train_2225$PC, differences = 1)
DTM <- diff(train_2225$TM, differences = 1)
DRH <- diff(train_2225$RH, differences = 1)
DWS <- diff(train_2225$WS, differences = 1)

# Gabungkan dalam dataframe
df <- data.frame(DTMA, DPC, DTM, DRH, DWS)
df.ts<-ts(df)
head(df)
##         DTMA   DPC   DTM   DRH   DWS
## 1  -7.708333  0.94  0.12  1.75  0.17
## 2  -7.500000 -1.39  0.54 -4.93 -0.69
## 3  -4.375000 -0.18  0.16 -2.91 -0.45
## 4  -1.458333  5.48  0.57  2.87 -0.32
## 5   9.583333  3.75 -0.58  2.39  0.59
## 6 -10.208333  6.72 -0.80  4.06  0.54

6.2 Plot ACF

acf_DRH <- ggAcf(DRH) + labs(title = "ACF DRH")
acf_DRH

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

6.3 Plot PACF

pacf_DRH <- ggPacf(DRH) + labs(title = "PACF DRH")
pacf_DRH

Uji ADF

adf.test(DRH)
## Warning in adf.test(DRH): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  DRH
## Dickey-Fuller = -13.91, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

6.5 Plot CCF (Lag optimum peubah penjelas)

ccf_result3 <- ggCcf(DTMA, DRH, lag.max=9, plot = TRUE)+ggtitle("CCF antara DTMA dan DRH")
ccf_result3

ccf_result1 <- ggCcf(DTMA, DRH, lag.max = 9) +
  ggtitle("CCF antara DTMA dan DRH") +
  theme(
    axis.text.x = element_text(size = 14),  # Besarkan angka di sumbu X
    axis.text.y = element_text(size = 14)   # Besarkan angka di sumbu Y
  )

ccf_result1

Dynardl (dynamac)

Uji Kointegrasi (Bounds Test)

modelbound <- ardlBound(
  data    = as.data.frame(train_2225[, c("TMA", "RH")]),
  formula = TMA ~ RH,
  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.72021, df1 = 1, df2 = 911, p-value = 0.3963
## 
## ------------------------------------------------------ 
## 
##  Ljung-Box Test for the autocorrelation in residuals:
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 0.011403, df = 1, p-value = 0.915
## 
## ------------------------------------------------------ 
## 
##  Breusch-Pagan Test for the homoskedasticity of residuals:
## 
##  studentized Breusch-Pagan test
## 
## data:  modelFull$model
## BP = 28.001, df = 13, p-value = 0.009047
## 
## The p-value of Breusch-Pagan test for the homoskedasticity of residuals:  0.009046565 < 0.05!
## ------------------------------------------------------ 
## 
##  Shapiro-Wilk test of normality of residuals:
## 
##  Shapiro-Wilk normality test
## 
## data:  modelFull$model$residual
## W = 0.88414, p-value < 2.2e-16
## 
## The p-value of Shapiro-Wilk test normality of residuals:  9.937552e-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 = 12.5612281845871 
##   
##  ------------------------------------------------------ 
##  F-statistic note: Asymptotic critical values used. 
##  
## ------------------------------------------------------ 
## 
##  Ramsey's RESET Test for model specification:
## 
##  RESET test
## 
## data:  modelECM$model
## RESET = 0.1483, df1 = 1, df2 = 913, p-value = 0.7003
## 
## ------------------------------------------------------
## ------------------------------------------------------ 
## 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 
## -47.265 -13.466  -3.586   7.973 146.539 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## ec.1   -0.15913    0.02589  -6.145 1.19e-09 ***
## dRH.t   1.51125    0.30519   4.952 8.75e-07 ***
## dRH.1   1.16648    0.31250   3.733 0.000201 ***
## dTMA.1 -0.36909    0.03683 -10.022  < 2e-16 ***
## dTMA.2 -0.22712    0.03788  -5.996 2.91e-09 ***
## dTMA.3 -0.12968    0.03848  -3.370 0.000783 ***
## dTMA.4 -0.08047    0.03843  -2.094 0.036546 *  
## dTMA.5 -0.11085    0.03803  -2.915 0.003649 ** 
## dTMA.6 -0.09482    0.03769  -2.516 0.012051 *  
## dTMA.7 -0.07525    0.03701  -2.033 0.042304 *  
## dTMA.8 -0.03551    0.03548  -1.001 0.317201    
## dTMA.9  0.02205    0.03196   0.690 0.490471    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.1 on 914 degrees of freedom
## Multiple R-squared:  0.2548, Adjusted R-squared:  0.245 
## F-statistic: 26.04 on 12 and 914 DF,  p-value: < 2.2e-16
## 
## ------------------------------------------------------ 
## Long-run coefficients: 
##      TMA.1       RH.1 
## -0.1591344  1.0241813 
##

Kointegrasi terjadi jika F-statistic > batas atas (I(1)) critical value.

Model

Jangka Pendek dan Jangka Panjang

modelardl <- dynardl(
  DTMA ~ DRH,
  data = df,
  lags = list(
    "DTMA" = 1:9,
    "DRH"  = 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 
## -50.489 -12.860  -3.271   6.958 148.496 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.082860   0.780800   0.106 0.915509    
## l.1.DTMA    -1.485002   0.032839 -45.220  < 2e-16 ***
## l.2.DTMA    -0.317888   0.036136  -8.797  < 2e-16 ***
## l.3.DTMA    -0.217508   0.037437  -5.810 8.62e-09 ***
## l.4.DTMA    -0.156467   0.037898  -4.129 3.98e-05 ***
## l.5.DTMA    -0.178432   0.037740  -4.728 2.62e-06 ***
## l.6.DTMA    -0.145611   0.037801  -3.852 0.000125 ***
## l.7.DTMA    -0.115795   0.037439  -3.093 0.002042 ** 
## l.8.DTMA    -0.070884   0.036129  -1.962 0.050068 .  
## l.9.DTMA     0.003252   0.032739   0.099 0.920900    
## l.1.DRH      1.270363   0.296897   4.279 2.08e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.76 on 915 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.7056, Adjusted R-squared:  0.7024 
## F-statistic: 219.3 on 10 and 915 DF,  p-value: < 2.2e-16
AIC(modelardl$model)
## [1] 8507.914

Uji Asumsi

Uji autokorelasi residual (Box-Ljung Test)

res1 <- residuals(modelardl$model)
Box.test(res1, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  res1
## X-squared = 18.846, df = 20, p-value = 0.5318

H0: tidak ada autokorelasi dalam residual; H1: terdapat autokorelasi dalam residual. Jika p-value > 0.05, tak tolak H0.

Uji Normalitas (Jarque-Bera Test)

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

Uji Heteroskedastisitas (Breusch-Pagan Test)

bptest(modelardl$model)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelardl$model
## BP = 18.654, df = 10, p-value = 0.04489

Multikolinearitas

Dengan satu peubah penjelas (RH) di model, multikolinearitas antar peubah penjelas tidak perlu diuji. Lag DRH dan lag DTMA dalam model tetap dapat dicek dengan VIF bila diperlukan.

vif(modelardl$model)
## l.1.DTMA l.2.DTMA l.3.DTMA l.4.DTMA l.5.DTMA l.6.DTMA l.7.DTMA l.8.DTMA 
## 1.249778 1.513953 1.625276 1.665838 1.652173 1.657539 1.625984 1.514292 
## l.9.DTMA  l.1.DRH 
## 1.243479 1.020193

Forecasting

Data Training

# Differencing
DTMA <- diff(train_2225$TMA, differences = 1)
DRH  <- diff(train_2225$RH,  differences = 1)

train_baru <- data.frame(DRH)

coefs     <- summary(modelardl$model)$coefficients
intercept <- coefs["(Intercept)", "Estimate"]
b1        <- coefs["l.1.DRH", "Estimate"]

train_baru$DTMA_pred <- intercept + b1 * train_baru$DRH

train_baru$TMA_pred <- NA
train_baru$TMA_pred[1] <- train_2225$TMA[1] + train_baru$DTMA_pred[1]
for (i in 2:nrow(train_baru)) {
  train_baru$TMA_pred[i] <- train_baru$TMA_pred[i - 1] + train_baru$DTMA_pred[i]
}

TMA_actual <- train_2225$TMA[-1]  
TMA_pred   <- train_baru$TMA_pred

plot(TMA_actual, type = "l", col = "black", main = "Train: Actual vs Predicted TMA")
lines(TMA_pred, col = "blue")
legend("topleft", legend = c("Actual", "Predicted"), col = c("black", "blue"), lty = 1)

# Evaluasi
MAE  <- mean(abs(TMA_pred - TMA_actual), na.rm = TRUE)
RMSE <- sqrt(mean((TMA_pred - TMA_actual)^2, na.rm = TRUE))
MAPE <- mean(abs((TMA_actual - TMA_pred) / TMA_actual), na.rm = TRUE) * 100

cat("MAE :", MAE, "\n")
## MAE : 31.5902
cat("RMSE:", RMSE, "\n")
## RMSE: 40.80613
cat("MAPE:", MAPE, "%\n")
## MAPE: 5.095374 %

Data Testing

DRH_test <- diff(c(tail(train_2225$RH, 1), test_2225$RH))

test_baru <- data.frame(DRH = DRH_test)
test_baru$DTMA_pred <- intercept + b1 * test_baru$DRH

test_baru$TMA_pred <- tail(train_2225$TMA, 1) + cumsum(test_baru$DTMA_pred)

TMA_actual_test <- test_2225$TMA
TMA_pred_test   <- test_baru$TMA_pred

plot(TMA_actual_test, type = "l", col = "black",
     ylim = range(c(TMA_actual_test, TMA_pred_test)),
     main = "Test: Actual vs Predicted TMA")
lines(TMA_pred_test, col = "red")
legend("topleft", legend = c("Actual", "Predicted"), col = c("black", "red"), lty = 1)

MAE_t  <- mean(abs(TMA_pred_test - TMA_actual_test))
RMSE_t <- sqrt(mean((TMA_pred_test - TMA_actual_test)^2))
MAPE_t <- mean(abs((TMA_actual_test - TMA_pred_test) / TMA_actual_test)) * 100

cat("MAE :", MAE_t, "\n")
## MAE : 19.67632
cat("RMSE:", RMSE_t, "\n")
## RMSE: 27.58166
cat("MAPE:", MAPE_t, "%\n")
## MAPE: 3.056368 %