Tugas Individu MPDW Pertemuan 3

Library

library(readxl)
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)
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.3

## corrplot 0.95 loaded

library(forecast)

## Warning: package 'forecast' was built under R version 4.5.3

library(urca)

## Warning: package 'urca' was built under R version 4.5.3

library(tseries)

## Warning: package 'tseries' was built under R version 4.5.3

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

## Warning: package 'zoo' was built under R version 4.5.3

## 
## 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)

## Warning: package 'lmtest' was built under R version 4.5.3

library(car)

## Warning: package 'car' was built under R version 4.5.3

## 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)
library(MASS)

## Warning: package 'MASS' was built under R version 4.5.3

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(patchwork)

## 
## 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.3

library(openxlsx)

## Warning: package 'openxlsx' was built under R version 4.5.3

## Registered S3 method overwritten by 'openxlsx':
##   method               from         
##   as.character.formula formula.tools

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

Import Data

Data

dataset <- read_excel("C:\\Users\\bagan\\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 ...

Menentukan Variable

Y <- dataset$TMA
X <- dataset$PC 
data_model_bersih <- data.frame(Y = Y, X = X)

Check Missing Value

sum(is.na(data_model_bersih$Y))

## [1] 0

sum(is.na(data_model_bersih$X))

## [1] 0

SPlitting Data

# (80% untuk Train)
batas <- round(0.8 * nrow(data_model_bersih))

# Memecah data
data_train <- data_model_bersih[1:batas, ]
data_test <- data_model_bersih[(batas+1):nrow(data_model_bersih), ]

# Melihat jumlah data latih dan uji
cat("Jumlah Data Latih:", nrow(data_train), "\n")

## Jumlah Data Latih: 877

cat("Jumlah Data Uji:", nrow(data_test), "\n")

## Jumlah Data Uji: 219

Plot ACF

acf(data_train$Y, main = "ACF Y")

Interpretasi: Berdasarkan plot ACF, nilai autokorelasi Y positif dan signifikan hingga lag ke-30, karena seluruh spike berada di luar batas signifikansi. Pola ACF yang menurun secara perlahan menunjukkan bahwa Y memiliki autokorelasi yang kuat dan kemungkinan belum stasioner.

Plot PACF

pacf(data_train$Y, main = "PACF Y")

Interpretasi: Berdasarkan plot PACF, terdapat autokorelasi parsial yang signifikan pada beberapa lag awal, terutama lag 1 hingga lag 4. Setelah lag tersebut, sebagian besar nilai PACF berada dalam batas signifikansi. Hal ini menunjukkan bahwa Y dipengaruhi terutama oleh beberapa periode sebelumnya, dengan pengaruh terbesar pada lag 1.

Model Autoregressive Distributed Lag (ARDL)

# Model ARDL dengan lag X = 1 dan lag Y = 1
model.ardl <- ardlDlm(
  formula = Y ~ X,
  data = data_train,
  p = 1,
  q = 1
)

summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 2, End = 877
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -84.681 -10.972  -3.422  10.407 119.069 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 190.62733   12.61895  15.106  < 2e-16 ***
## X.t           0.96670    0.12273   7.876 1.00e-14 ***
## X.1           0.95862    0.13110   7.312 5.94e-13 ***
## Y.1           0.67589    0.02066  32.709  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.91 on 872 degrees of freedom
## Multiple R-squared:  0.7032, Adjusted R-squared:  0.7022 
## F-statistic: 688.7 on 3 and 872 DF,  p-value: < 2.2e-16

Interpretasi: Berdasarkan hasil pemodelan ARDL dengan \(p=1\) dan \(q=1\), diperoleh nilai \(R^2=0,7032\), yang berarti model mampu menjelaskan 70,32% variasi Y. Variabel \(Y_{t-1}\) berpengaruh positif dan signifikan terhadap \(Y_t\) dengan koefisien 0,67589 dan p-value < 0,05. Sementara itu, \(X_t\) dan \(X_{t-1}\) juga berpengaruh positif dan signifikan terhadap \(Y_t\), dengan koefisien masing-masing 0,96670 dan 0,95862, serta p-value < 0,05.

Mencari Lag Optimum

# Mencari kombinasi lag optimum ARDL
model.ardl.opt <- ardlBoundOrders(
  data = data.frame(data_train),
  ic = "AIC",
  formula = Y ~ X
)

# Menentukan kombinasi lag dengan AIC terkecil
min_p <- c()

for (i in 1:6) {
  min_p[i] <- min(model.ardl.opt$Stat.table[[i]])
}

q_opt <- which(
  min_p == min(min_p, na.rm = TRUE)
)

p_opt <- which(
  model.ardl.opt$Stat.table[[q_opt]] ==
    min(model.ardl.opt$Stat.table[[q_opt]], na.rm = TRUE)
)

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

##   Lag_Y_opt.q. Lag_X_opt.p.     AIC
## 1            3           15 7675.78

Interpretasi: Berdasarkan hasil pemilihan model dengan nilai AIC terkecil di dapat lag optimum Y = 3 dan lag optimum X = 15

Pemodelan menggunakan Lag Optimum

model_ardl_final <- ardlDlm(
  formula = Y ~ X,
  data = data_train,
  p = p_opt,
  q = q_opt
)

summary(model_ardl_final)

## 
## Time series regression with "ts" data:
## Start = 16, End = 877
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -58.068 -11.420  -3.122   7.668 121.110 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 109.69342   16.04472   6.837 1.56e-11 ***
## X.t           0.97959    0.11758   8.331 3.23e-16 ***
## X.1           1.41410    0.12963  10.909  < 2e-16 ***
## X.2          -0.57739    0.13798  -4.185 3.16e-05 ***
## X.3          -0.23856    0.13677  -1.744   0.0815 .  
## X.4          -0.17455    0.13687  -1.275   0.2026    
## X.5           0.06385    0.12703   0.503   0.6153    
## X.6          -0.15061    0.12714  -1.185   0.2365    
## X.7           0.00745    0.12701   0.059   0.9532    
## X.8           0.14307    0.12680   1.128   0.2595    
## X.9           0.13961    0.12644   1.104   0.2699    
## X.10          0.12814    0.12650   1.013   0.3114    
## X.11         -0.21077    0.12652  -1.666   0.0961 .  
## X.12          0.04798    0.12707   0.378   0.7058    
## X.13         -0.04427    0.12596  -0.351   0.7253    
## X.14          0.07535    0.12534   0.601   0.5479    
## X.15         -0.16425    0.11791  -1.393   0.1640    
## Y.1           0.42091    0.03388  12.425  < 2e-16 ***
## Y.2           0.20868    0.03630   5.749 1.25e-08 ***
## Y.3           0.18083    0.03395   5.326 1.29e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.01 on 842 degrees of freedom
## Multiple R-squared:  0.7548, Adjusted R-squared:  0.7493 
## F-statistic: 136.5 on 19 and 842 DF,  p-value: < 2.2e-16

Interpretasi: Berdasarkan hasil pemodelan ARDL dengan lag \(Y=3\) dan \(X=15\), diperoleh \(R^2=0,7548\), yang berarti model mampu menjelaskan 75,48% variasi Y. Variabel \(X_t\) dan \(X_{t-1}\) berpengaruh positif signifikan, sedangkan \(X_{t-2}\) berpengaruh negatif signifikan terhadap \(Y_t\). Variabel \(Y_{t-1}\), \(Y_{t-2}\), dan \(Y_{t-3}\) juga berpengaruh positif dan signifikan terhadap \(Y_t\). Secara simultan, model signifikan dengan p-value < 0,05.

Peramalan dengan Lag Optimum

fore.ardl <- forecast(
  model = model_ardl_final,
  x = data_test$X,
  h = nrow(data_test)
)

fore.ardl

## $forecasts
##   [1] 649.2949 644.3122 636.8464 632.0122 623.2128 622.8804 626.8814 673.2213
##   [9] 701.5060 638.2109 634.1257 628.8342 620.7670 610.1675 614.8382 621.1842
##  [17] 627.9722 622.0932 607.2234 624.7901 633.9142 622.5549 603.7989 607.9755
##  [25] 603.2512 598.5477 595.7151 596.5802 597.0670 603.1723 603.1450 595.1158
##  [33] 598.9655 615.3259 652.8417 666.9728 785.0090 920.0697 763.0500 741.8652
##  [41] 739.2525 749.0699 698.9066 703.9656 707.9724 696.5570 694.8808 652.8104
##  [49] 659.4495 642.3159 645.6927 611.4304 613.9149 608.9883 603.1480 598.1874
##  [57] 596.3610 594.0470 591.2411 590.3104 588.9690 587.6917 586.7060 585.7710
##  [65] 584.9549 584.2870 592.4699 627.2480 651.8130 629.9165 618.7536 610.6533
##  [73] 605.1143 599.2537 598.7050 601.2582 604.3837 602.7883 595.2348 594.0656
##  [81] 591.8748 590.3619 585.3170 584.5223 583.8022 583.1101 582.5628 582.0965
##  [89] 581.6439 581.6634 581.8104 581.2723 581.0969 581.2094 580.7228 580.7401
##  [97] 580.1733 580.1051 580.5033 581.0265 582.3834 585.4289 583.3520 581.5115
## [105] 582.4206 596.7199 623.8874 628.5602 624.7966 631.2340 620.7092 608.2208
## [113] 608.5020 604.2651 603.6026 603.1531 601.1636 598.3762 598.1491 599.9957
## [121] 671.6389 739.7841 663.9010 677.4229 652.5241 642.0539 630.1130 626.5827
## [129] 629.1260 631.3512 642.2905 631.1743 630.6898 617.0494 615.9582 600.1213
## [137] 600.8474 598.3956 594.7110 596.3769 602.4613 606.4195 599.8026 593.4446
## [145] 592.7220 590.0004 588.0839 586.8825 588.1006 593.1694 609.7283 629.3834
## [153] 622.1742 603.5285 607.3665 604.0310 599.6850 598.2664 603.1268 623.4754
## [161] 639.4839 621.2143 627.1047 638.3450 626.3344 608.9186 627.1638 648.8706
## [169] 641.7349 649.7484 626.9287 622.8544 618.5560 625.1147 638.8240 669.8024
## [177] 678.4356 670.0462 682.3008 677.1119 662.8573 677.7900 669.3204 652.5893
## [185] 667.3681 712.5827 726.2993 659.0620 668.3499 671.8605 699.0632 730.4988
## [193] 712.2053 699.8297 694.5050 693.8806 674.2873 691.5600 667.5972 660.7668
## [201] 656.4048 659.5705 656.7666 657.6622 663.3591 642.4810 665.7150 678.5008
## [209] 655.8600 647.2109 659.2361 659.9049 648.9182 650.6700 641.6714 637.4699
## [217] 627.8224 626.7331 626.1894
## 
## $call
## forecast.ardlDlm(model = model_ardl_final, x = data_test$X, h = nrow(data_test))
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

# Menghitung MAPE
mape.ardl <- MLmetrics::MAPE(
  fore.ardl$forecasts,
  data_test$Y
)

cat("MAPE Model ARDL Optimum:", mape.ardl, "\n")

## MAPE Model ARDL Optimum: 0.03560578

Visualisasi

library(ggplot2)

# 1. Buat dataframe khusus visualisasi
df_plot <- data.frame(
  Periode = 1:nrow(data_test), 
  Aktual = data_test$Y,
  Forecast = fore.ardl$forecasts
)

ggplot(df_plot, aes(x = Periode)) +
  geom_line(aes(y = Aktual, color = "Aktual"), linewidth = 1) +
  geom_line(aes(y = Forecast, color = "Forecast (ARDL)"), linewidth = 1, linetype = "dashed") +
  scale_color_manual(values = c("Aktual" = "#1F77B4", "Forecast (ARDL)" = "#D62728")) +
  labs(
    title = "Perbandingan Nilai Aktual vs Forecast Model ARDL",
    subtitle = paste0("MAPE: ", round(mape.ardl * 100, 2), "%"),
    x = "Periode",
    y = "Nilai Y",
    color = "Keterangan"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "top"
  )

Interpretasi: Berdasarkan grafik perbandingan nilai aktual dan forecast model ARDL, diperoleh nilai MAPE sebesar 3,56%. Hal ini menunjukkan bahwa hasil peramalan model ARDL memiliki tingkat kesalahan yang relatif kecil, sehingga nilai forecast cukup mendekati nilai aktual.

Uji Asumsi Model ARDL

# 1. Uji Normalitas
shapiro.test(residuals(model_ardl_final$model))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_ardl_final$model)
## W = 0.91679, p-value < 2.2e-16

# 2. Uji Autokorelasi 
bgtest(model_ardl_final$model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model_ardl_final$model
## LM test = 10.074, df = 1, p-value = 0.001504

# 3. Uji Homoskedastisitas 
bptest(model_ardl_final$model)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_ardl_final$model
## BP = 86.188, df = 19, p-value = 1.559e-10

Interpretasi: Uji Shapiro-Wilk menunjukkan residual tidak berdistribusi normal karena p-value < 0,05. Uji Breusch-Godfrey menunjukkan terdapat autokorelasi pada residual karena p-value < 0,05. Uji Breusch-Pagan menunjukkan terjadi heteroskedastisitas karena p-value < 0,05. Jadi, ketiga asumsi residual belum terpenuhi.