Tugas STA1341 : Regresi dengan Peubah Lag

Package

dlagM
dynlM
MLmetrics
lmtest
car
readxl

library(dLagM) #bisa otomatis timeseries rainfallnya

## Loading required package: nardl

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Loading required package: dynlm

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(dynlm) #rainfall harus timeseries
library(MLmetrics) #MAPE

## 
## Attaching package: 'MLmetrics'

## The following object is masked from 'package:dLagM':
## 
##     MAPE

## The following object is masked from 'package:base':
## 
##     Recall

library(lmtest)
library(car)

## Loading required package: carData

library(readxl)

Data

Data yang digunakan merupakan data yang bersumber dari kaggle.com berupa project dari Prediction of Worldwide Energy Resources (POWER) yang menyediakan data meteorologi dari riset NASA untuk mendukung energi terbarukan, membangun efisiensi energi, dan mendukung bidang agrikultur. Periode waktu yang digunakan dalam data tersebut mulai dari 1 Januari 2000 hingga 1 Desember 2020. Data terdiri atas dua peubah, yaitu Specific Humidity sebagai peubah respon (Y) dan Temperature sebagai peubah penjelas (X).

rainfall <- read_excel("C:/Users/LENOVO/Documents/MK SEM 5/STA1341 - MPDW/Kuliah/4/Data Hujan.xlsx", sheet = 3)
str(rainfall)

## tibble [252 x 3] (S3: tbl_df/tbl/data.frame)
##  $ t : num [1:252] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Yt: num [1:252] 8.06 8.73 8.48 13.79 17.4 ...
##  $ Xt: num [1:252] 23.9 25.8 26.7 22.5 19.1 ...

rainfall

## # A tibble: 252 x 3
##        t    Yt    Xt
##    <dbl> <dbl> <dbl>
##  1     1  8.06 23.9 
##  2     2  8.73 25.8 
##  3     3  8.48 26.7 
##  4     4 13.8  22.5 
##  5     5 17.4  19.1 
##  6     6 19.5   7.91
##  7     7 18.8   6.67
##  8     8 18.9   7.07
##  9     9 18.4  10.6 
## 10    10 16.7  15.4 
## # ... with 242 more rows

#scaterplot Y
rainfall1 <- rainfall[,-1]
rainfall1 <- rainfall1[,-2]
rainfall1.ts <- ts(rainfall1)
ts.plot(rainfall1.ts, 
        xlab="Specific Humidity", 
        ylab="Time Period", 
        main= "Time Series Plot of Specific Humidity 2000 - 2020")
points(rainfall1.ts)

Kelembaban spesifik atau rasio kelembaban adalah besaran massa uap air yang terkandung di udara per satuan massa udara kering yang diukur dalam gram per kilogram dari udara kering (Widodo et al. 2008).

#scatterplot X
rainfall2 <- rainfall[,-1:-2]
rainfall2.ts <- ts(rainfall2)
ts.plot(rainfall2.ts, 
        xlab="Temperature", 
        ylab="Time Period", 
        main= "Plot Time Series of Temperature 2000 - 2020")
points(rainfall2.ts)

Suhu adalah derajat panas atau dingin yang diukur berdasarkan skala tertentu dengan menggunakan termometer. Satuan suhu yang biasa digunakan adalah derajat celcius (0C). Sedangkan di Inggris dan beberapa Negara lainnya dinyatakan dalam derajat Fahrenheit (0F). Suhu juga bisa diartikan sebagai suatu sifat fisika dari suatu benda yang menggambarkan Energy kinetic rata-rata dari pergerakan molekul-molekul.

# Diagram pencar identifikasi model
plot(rainfall$Xt,rainfall$Yt, pch = 20, 
     col = "steelblue", 
     main = "Scatter Plot Specify Humidity vs Temperature",
     ylab = "Specify Humidity", xlab = "rainfall")

Berdasarkan scatterplot di atas, terlihat hubungan antara specific humidity dan temperature memiliki hubungan linear negatif. Arti dari hubungan negatif tersebut adalah peningkatan suhu akan menurunkan kelembaban spesifik.

# Korelasi rainfall dan Income
cor(rainfall$Xt, rainfall$Yt)

## [1] -0.9058671

Berdasarkan nilai korelasi di atas dapat diketahui bahwa nilai antara kedua peubah mendekati -1, yaitu -0.9058671. Artinya, terdapat hubungan linear negatif yang kuat antara peubah kelembaban spesifik dan peubah temperatur tersebut.

Splitting Data

Pada proses splitting data dilakukan pemotongan untuk data training yang digunakan untuk melatih aloritma dalam mencari model yang sesuai dan data testing yang digunakan untuk menguji performa model yang didapatkan pada tahapan testing dengan pembagian sebesar 80:20 dari jumlah amatan yang ada yaitu 252 amatan.

#SPLIT rainfall 80:20
train<-rainfall[1:202,]
test<-rainfall[203:252,]

#rainfall time series
train.ts<-ts(train)
test.ts<-ts(test)
rainfall.ts<-ts(rainfall)

Pada proses splitting dilakukan pemotongan untuk data training dan data testing dengan pembagian sebesar 80:20 dari jumlah amatan yang ada yaitu 252 amatan.

MODEL KOYCK

Metode Koyck digunakan untuk penentuan estimasi model dinamis terdistribusi lag ketika panjang (lag) tidak diketahui dan asumsikan bahwa semakin jauh jarak lag pada variabel bebas dari periode sekarang maka semakin kecil pengaruh variabel lag terhadap variabel tak bebas (Lihawa et al. 2022).

#MODEL KOYCK
model.koyck = dLagM::koyckDlm(x = train$Xt, y = train$Yt)
summary(model.koyck)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.461  -9.049  -1.078   8.103  27.931 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  164.902    137.946   1.195    0.233
## Y.1           -5.319      5.228  -1.017    0.310
## X.t           -4.549      3.860  -1.178    0.240
## 
## Residual standard error: 11.78 on 198 degrees of freedom
## Multiple R-Squared: -6.282,  Adjusted R-squared: -6.355 
## Wald test: 10.33 on 2 and 198 DF,  p-value: 5.386e-05 
## 
## Diagnostic tests:
## NULL
## 
##                            alpha      beta      phi
## Geometric coefficients:  26.0948 -4.548847 -5.31934

Berdasarkan output di atas dapat disimpulkan bahwa pemodelan koyck yang dibentuk dari data rainfall antara peubah kelembaban spesifik (Y) dan Temperature (X) adalah Yt = 164.902 - 5.319Yt - 4.549Xt. Selain itu diketahui juga semua peubah memiliki nilai p-value > 0.05 yang mengindikasikan bahwa penduga parameter tidak signifikan.

AIC(model.koyck)

## [1] 1566.761

BIC(model.koyck)

## [1] 1579.975

#ramalan
(fore.koyck <- forecast(model = model.koyck, x=test$Xt, h=50))

## $forecasts
##  [1] -1.565820e+00  8.052551e+01 -3.725220e+02  2.033070e+03 -1.077428e+04
##  [6]  5.737225e+04 -3.050945e+05  1.623011e+06 -8.633209e+06  4.592310e+07
## [11] -2.442805e+08  1.299411e+09 -6.912007e+09  3.676731e+10 -1.955778e+11
## [16]  1.040345e+12 -5.533948e+12  2.943695e+13 -1.565851e+14  8.329294e+14
## [21] -4.430634e+15  2.356805e+16 -1.253664e+17  6.668667e+17 -3.547290e+18
## [26]  1.886924e+19 -1.003719e+20  5.339123e+20 -2.840061e+21  1.510725e+22
## [31] -8.036057e+22  4.274652e+23 -2.273832e+24  1.209529e+25 -6.433894e+25
## [36]  3.422407e+26 -1.820494e+27  9.683827e+27 -5.151157e+28  2.740075e+29
## [41] -1.457539e+30  7.753145e+30 -4.124161e+31  2.193781e+32 -1.166947e+33
## [46]  6.207386e+33 -3.301919e+34  1.756403e+35 -9.342904e+35  4.969808e+36
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 50)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

Plot Data Aktual vs Model Koyck

plot(test$Xt, test$Yt, type="b", col="black", main="Data Aktual vs Model Koyck")
points(test$Xt, fore.koyck$forecasts,col="red")
lines(test$Xt, fore.koyck$forecasts,col="red")

Berdasarkan plot diatas, dapat dilihat bahwa data hasil peramalan menggunakan model koyck sangat tidak mendekati data aktual. Data aktual memiliki rentang nilai 5 sampai 20, sementara data hasil peramalan dengan model koyck memiliki rentang nilai yang sangat jauh, yaitu -9.342904e+35 sampai 9.683827e+27.

#mape rainfall testing
mape.koyck <- MAPE(fore.koyck$forecasts, test$Yt)

#akurasi rainfall training
mape_train <- dLagM::GoF(model.koyck)["MAPE"]

c("MAPE_testing" = mape.koyck, "MAPE_training" = mape_train)

## $MAPE_testing
## [1] 9.632284e+33
## 
## $MAPE_training.MAPE
## [1] 0.7065525

Berdasarkan output di atas, dapat dilihat bahwa nilai MAPE testing lebih besar dibandingkan dengan nilai MAPE training. Nilai MAPE testing dan MAPE training memiliki perbedaan yang sangat signifikan, dengan demikian artinya bahwa model koyck ini underfitted atau overfitted.

Regression with Distributed Lag

Model regresi yang memuat variabel tak bebas yang dipengaruhi oleh variabel bebas pada waktu t, serta dipengaruhi juga oleh variabel bebas pada waktu t-1, t-2, …, t-s disebut model distribusi lag, sebab pengaruh dari satu atau beberapa variabel bebas (X) terhadap variabel tak bebas (Y) menyebar ke beberapa periode waktu.

Regression with Distributed Lag (lag=2)

model.dlm = dLagM::dlm(x = train$Xt,y = train$Yt , q = 2)
summary(model.dlm)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.518 -1.549 -0.116  1.467  5.049 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 23.36258    0.41070  56.885   <2e-16 ***
## x.t         -0.63286    0.04047 -15.637   <2e-16 ***
## x.1         -0.00559    0.06127  -0.091   0.9274    
## x.2          0.08827    0.04064   2.172   0.0311 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.779 on 196 degrees of freedom
## Multiple R-squared:  0.8341, Adjusted R-squared:  0.8315 
## F-statistic: 328.4 on 3 and 196 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 803.9586 820.4502

Berdasarkan output di atas dapat disimpulkan bahwa model dlm lag 2 seluruh peubah memiliki nilai p-value < 0.05 sehingga signifikan kecuali xt-1 yang memiliki nilai p-value > 0.05 sehingga tidak signifikan.

AIC(model.dlm)

## [1] 803.9586

BIC(model.dlm)

## [1] 820.4502

#ramalan
(fore.dlm <- forecast(model = model.dlm, x=test$Xt, h=50))

## $forecasts
##  [1] 12.339499 11.725818  9.696853  9.250291  8.005926 10.841645 14.950022
##  [8] 17.693118 21.389636 20.039527 17.298656 13.843652 14.086059 14.020455
## [15] 11.329344  9.485561 10.248988 12.390663 14.610035 16.817880 20.998715
## [22] 20.930775 16.138178 12.948662 13.471187 10.868359  9.940945  8.021405
## [29]  8.715571 10.407158 13.881114 16.703413 20.874291 20.638726 19.651158
## [36] 16.318846 15.248753 13.441493  9.928077  8.951620  9.512242 11.447826
## [43] 14.909971 17.777133 20.284177 20.981625 18.122720 15.795698 13.378900
## [50] 13.135879
## 
## $call
## forecast.dlm(model = model.dlm, x = test$Xt, h = 50)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

Plot Data Aktual vs Model DLM Lag = 2

plot(test$Xt, test$Yt, type="b", col="black", main="Data Aktual vs Model DLM Lag 2" )
points(test$Xt, fore.dlm$forecasts,col="blue")
lines(test$Xt, fore.dlm$forecasts,col="blue")

Berdasarkan plot diatas, dapat dilihat bahwa data hasil peramalan menggunakan model dlm 2 hampir mendekati data aktual. Data aktual memiliki rentang nilai 5 sampai 20, sementara data hasil peramalan dengan model dlm 2 memiliki rentang nilai yang hampir sama, yaitu 8 sampai 20.

#mape rainfall testing
mape.dlm <- MAPE(fore.dlm$forecasts, test$Yt)

#akurasi rainfall training
mape_train <- GoF(model.dlm)["MAPE"]

c("MAPE_testing" = mape.dlm, "MAPE_training" = mape_train)

## $MAPE_testing
## [1] 0.1283605
## 
## $MAPE_training.MAPE
## [1] 0.11451

Berdasarkan output di atas, dapat dilihat bahwa nilai MAPE testing lebih besar dibandingkan dengan nilai MAPE training. Nilai MAPE testing dan MAPE training memiliki perbedaan yang tidak terlalu signifikan, dengan demikian artinya bahwa model dlm lag 2 ini tidak underfitted atau overfitted.

Penentuan Lag Optimum

finiteDLMauto(formula = Yt ~ Xt,
              data = data.frame(train), q.min = 1, q.max = 4 ,
              model.type = "dlm", error.type = "AIC", trace = TRUE) ##q max lag maksimum

##   q - k    MASE      AIC      BIC   GMRAE    MBRAE R.Adj.Sq    Ljung-Box
## 4     4 0.62011 732.1474 755.1653 0.75257  0.88787  0.87918 1.188254e-05
## 3     3 0.71460 779.1102 798.8700 0.91653  0.75968  0.84832 1.721757e-06
## 2     2 0.74814 803.9586 820.4502 0.93461  0.42274  0.83152 1.621071e-06
## 1     1 0.75652 809.7075 822.9207 0.99160 -0.05597  0.82984 8.562576e-06

Berdasarkaan output di atas, dapat dilihat bahwa persamaan optimum pada saat lag 4 (q = 4), dengan nilai AIC dan BIC yang terkecil dan nilai adjusted R-squared yang terbesar.

#model dlm dengan lag optimum
model.dlm2 = dLagM::dlm(x = train$Xt,y = train$Yt , q = 4) #terdapat lag yang tidak signifikan sehingga dapat dikurangi jumlah lagnya 
summary(model.dlm2)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3801 -1.1309 -0.0719  0.9887  4.1034 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.13239    0.59598  32.102  < 2e-16 ***
## x.t         -0.45365    0.03991 -11.366  < 2e-16 ***
## x.1          0.02001    0.05212   0.384   0.7014    
## x.2         -0.09178    0.05240  -1.752   0.0814 .  
## x.3         -0.04360    0.05216  -0.836   0.4043    
## x.4          0.27714    0.04001   6.926 6.35e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.507 on 192 degrees of freedom
## Multiple R-squared:  0.8822, Adjusted R-squared:  0.8792 
## F-statistic: 287.7 on 5 and 192 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 732.1474 755.1653

Berdasarkan output di atas dapat disimpulkan bahwa model dlm lag 4 seluruh peubah memiliki nilai p-value < 0.05 sehingga signifikan kecuali peubah xt-1 dan xt-3 yang memiliki nilai p-value > 0.05 sehingga tidak signifikan.

AIC(model.dlm2)

## [1] 732.1474

BIC(model.dlm2)

## [1] 755.1653

#ramalan
(fore.dlm2 <- forecast(model = model.dlm2, x=test$Xt, h=50))

## $forecasts
##  [1] 11.545353 10.103336  8.522088  9.911999  9.215866 11.552135 14.966969
##  [8] 17.620158 21.972138 20.658132 18.280699 14.586070 12.460817 11.810901
## [15] 10.590670 10.794439 10.727014 11.843018 14.509108 17.324288 20.733466
## [22] 20.450326 17.281622 14.446675 11.974529  8.955682 10.009880  9.288111
## [29]  9.457192 11.343179 14.161376 17.446519 21.027472 21.061167 20.137914
## [36] 16.810561 14.063527 11.913884  9.136101  9.395420  9.272979 11.087748
## [43] 15.027788 17.938812 20.559043 21.169298 18.502992 16.052295 12.578754
## [50] 11.329969
## 
## $call
## forecast.dlm(model = model.dlm2, x = test$Xt, h = 50)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

Plot Data Aktual vs Model DLM Lag = 4

plot(test$Xt, test$Yt, type="b", col="black", main="Data Aktual vs Model DLM Lag 4")
points(test$Xt, fore.dlm2$forecasts,col="orange")
lines(test$Xt, fore.dlm2$forecasts,col="orange")

Berdasarkan plot diatas, dapat dilihat bahwa data hasil peramalan menggunakan model dlm 4 hampir mendekati data aktual. Data aktual memiliki rentang nilai 5 sampai 20, sementara data hasil peramalan dengan model dlm 4 memiliki rentang nilai yang hampir sama, yaitu 8 sampai 21.

#akurasi testing
mape.dlm2 <- MAPE(fore.dlm2$forecasts, test$Yt)

#akurasi rainfall training
mape_train <- GoF(model.dlm2)["MAPE"]

c("MAPE_testing" = mape.dlm2, "MAPE_training" = mape_train)

## $MAPE_testing
## [1] 0.1100986
## 
## $MAPE_training.MAPE
## [1] 0.09446757

Berdasarkan output di atas, dapat dilihat bahwa nilai MAPE testing lebih besar dibandingkan dengan nilai MAPE training. Nilai MAPE testing dan MAPE training memiliki perbedaan yang tidak terlalu signifikan, dengan demikian artinya bahwa model dlm lag 4 ini tidak underfitted atau overfitted.

Model Autoregressive

Model Autoregressive(AR) merupakan salah satu model yang mengasumsikan bahwa data pada periode sekarang dipengaruhi oleh data pada periode sebelumnya (Fitriani et al. 2013).

#MODEL AUTOREGRESSIVE 
#library(dLagM)
model.ardl = ardlDlm(x = train$Xt, y = train$Yt, p = 1 , q = 1) #p:lag x, q:lag y
#model untuk p=1, q=1: yt=b0+b1yt-1+b2xt+b3xt-1
#model untuk p=2, q=3: yt=b0+b1yt-1+b2yt-2+b3xt+b4xt-1+b5xt-2

summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 2, End = 202
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0923 -1.1458 -0.0271  1.2168  4.5718 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.65084    2.15379   4.017 8.40e-05 ***
## X.t         -0.49371    0.03955 -12.482  < 2e-16 ***
## X.1          0.31560    0.04097   7.704 6.31e-13 ***
## Y.1          0.60178    0.08450   7.122 1.95e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.601 on 197 degrees of freedom
## Multiple R-squared:  0.866,  Adjusted R-squared:  0.864 
## F-statistic: 424.5 on 3 and 197 DF,  p-value: < 2.2e-16

Berdasarkan output di atas dapat disimpulkan bahwa model autoregressive yang dibentuk dari data rainfall antara peubah kelembaban spesifik (Y) dan Temperature untuk semua peubah memiliki nilai p-value < 0.05 yang mengindikasikan bahwa penduga parameter signifikan.

AIC(model.ardl)

## [1] 765.6621

BIC(model.ardl)

## [1] 782.1787

#ramalan
(fore.ardl <- forecast(model = model.ardl, x=test$Xt, h=50))

## $forecasts
##  [1] 13.806307 12.704444 10.888904 10.463480  9.292767 11.522126 14.501097
##  [8] 16.779705 19.890803 18.918774 17.136713 14.349674 14.432715 14.055912
## [15] 12.011786 10.641843 11.080206 12.603394 14.348382 16.182442 19.555287
## [22] 19.536390 16.154717 13.765510 13.858947 11.561161 10.987269  9.271439
## [29]  9.862591 10.986298 13.748747 15.969144 19.437740 19.320353 18.913679
## [36] 16.257285 15.460123 13.760701 11.020955 10.207599 10.395549 11.849599
## [43] 14.542895 16.844855 19.003714 19.684924 17.630137 15.939525 13.857430
## [50] 13.566168
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 50)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

Plot Data Aktual vs Model Autoregressive

plot(test$Xt, test$Yt, type="b", col="black", main="Data Aktual vs Model Autoregressive")
points(test$Xt, fore.ardl$forecasts,col="green")
lines(test$Xt, fore.ardl$forecasts,col="green")

Berdasarkan plot diatas, dapat dilihat bahwa data hasil peramalan menggunakan model autoregressive hampir mendekati data aktual. Data aktual memiliki rentang nilai 5 sampai 20, sementara data hasil peramalan dengan model autoregressive memiliki rentang nilai yang hampir sama, yaitu 9 sampai 20.

#akurasi testing
mape.ardl <- MAPE(fore.ardl$forecasts, test$Yt) #rainfall testing

#akurasi rainfall training
mape_train <- GoF(model.ardl)["MAPE"]

c("MAPE_testing" = mape.ardl, "MAPE_training" = mape_train)

## $MAPE_testing
## [1] 0.1493687
## 
## $MAPE_training.MAPE
## [1] 0.1009329

Berdasarkan output di atas, dapat dilihat bahwa nilai MAPE testing lebih besar dibandingkan dengan nilai MAPE training. Nilai MAPE testing dan MAPE training memiliki perbedaan yang tidak terlalu signifikan, dengan demikian artinya bahwa model autoregressive ini tidak underfitted atau overfitted.

#penentuan lag optimum
ardlBoundOrders(data = data.frame(rainfall) , formula = Yt ~ Xt ) #yang digunakan harusnya data train, tetapi karena keterbatasan data jika menggunakan data train menyebabkan error sehingga dicontohkan menggunakan keseluruhan data

## $p
##   Xt
## 1 11
## 
## $q
## [1] 15
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  821.7318 768.2995 735.0955 717.9291 715.6051 714.0595 712.2159 707.8496
## p = 2  795.9056 768.3062 737.0803 718.3067 716.0591 714.5613 712.2180 707.5992
## p = 3  746.6906 746.6906 735.0557 719.3957 717.6911 716.0451 714.1267 709.5782
## p = 4  724.8624 724.9997 724.9997 718.7086 717.8331 715.6758 713.5137 710.4713
## p = 5  707.8228 705.5583 706.9728 706.9728 708.8708 705.1079 698.5882 695.6146
## p = 6  710.4605 706.8435 706.2001 705.1655 705.1655 702.8424 692.6235 685.6752
## p = 7  703.9281 701.3817 701.3583 700.1345 702.0910 702.0910 693.3151 683.7601
## p = 8  646.1175 647.8571 649.6956 650.8766 652.8482 654.4888 654.4888 655.6525
## p = 9  626.8011 628.7967 630.5736 632.3638 634.2750 636.2725 636.6988 636.6988
## p = 10 626.4653 628.4030 630.2707 632.0449 633.9426 635.9378 635.6254 636.2009
## p = 11 610.1032 612.1027 614.0505 615.5168 617.5089 619.4721 619.3900 620.0434
## p = 12 587.4542 589.4439 591.3852 592.8169 594.8097 596.7302 597.0488 597.1561
## p = 13 610.5051 612.4892 614.4857 615.9614 617.9519 619.7716 620.0956 620.7371
## p = 14 604.5349 606.5345 608.5339 610.3879 612.3607 614.0382 614.5126 615.7036
## p = 15 597.6890 599.3933 600.7758 602.4750 603.8489 605.6739 606.8070 608.2086
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  698.8784 677.3433 653.3122 639.9131 629.8367 621.1121 612.2147
## p = 2  699.3227 677.7006 654.0172 640.2293 629.9587 622.4780 613.2531
## p = 3  701.2656 679.6731 656.0079 642.2265 631.8890 624.3846 615.2201
## p = 4  703.0506 681.6707 657.8810 644.0180 633.7950 626.3426 617.2092
## p = 5  692.5768 674.1361 652.7990 640.0081 628.6459 621.5469 613.4236
## p = 6  682.8035 669.6446 648.5838 634.3534 622.7451 615.2012 607.6944
## p = 7  678.4204 663.3791 645.8572 630.6783 618.6321 611.2316 602.1955
## p = 8  654.8926 648.3995 633.6966 621.9726 612.7075 606.9978 598.5437
## p = 9  638.0755 635.4399 626.0605 612.2615 604.7276 598.6913 590.1093
## p = 10 636.2009 635.9907 624.2739 612.8736 605.5086 599.3472 589.8893
## p = 11 621.7933 621.7933 602.1940 597.7306 592.5939 588.8697 578.7533
## p = 12 599.1404 601.1395 601.1395 599.7177 594.5111 590.4681 580.4976
## p = 13 622.7015 622.5079 596.9391 596.9391 596.4557 592.0525 581.9458
## p = 14 617.7021 618.1085 593.7870 593.8647 593.8647 593.8676 583.9411
## p = 15 610.2056 610.3491 587.3360 585.1852 586.9736 586.9736 585.5723
## 
## $min.Stat
## [1] 578.7533

Berdasarkan output di atas, dapat diketahui bahwa lag optimum berada pada saat p = 11 dan q = 15. Namun, pada metode lanjutan tidak digunakan lag optimum tersebut dikarenakan akan membuat koefisien pada model tersebut banyak yang tidak signifikan sehingga kurang baik untuk digunakan.

#PEMODELAN DLM dan ARDL dengan library dynlm
#library(dynlm)

#sama dengan model dlm p=1
cons_lm1 <- dynlm(Yt ~ Xt+L(Xt),data = train.ts)

#sama dengan model ardl p=0 q=1
cons_lm2 <- dynlm(Yt ~ Xt+L(Yt),data =  train.ts)

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

#sama dengan dlm p=2
cons_lm4 <- dynlm(Yt ~ Xt+L(Xt)+L(Xt,2),data =  train.ts)

#Ringkasan Model
summary(cons_lm1)

## 
## Time series regression with "ts" data:
## Start = 2, End = 202
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.534 -1.368 -0.300  1.545  5.682 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 23.82355    0.35306  67.477  < 2e-16 ***
## Xt          -0.68623    0.03230 -21.249  < 2e-16 ***
## L(Xt)        0.10795    0.03219   3.354 0.000955 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.791 on 198 degrees of freedom
## Multiple R-squared:  0.8315, Adjusted R-squared:  0.8298 
## F-statistic: 488.7 on 2 and 198 DF,  p-value: < 2.2e-16

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 202
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7283 -1.3686 -0.1845  1.3186  5.1010 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 20.83278    1.66386   12.52   <2e-16 ***
## Xt          -0.51660    0.04488  -11.51   <2e-16 ***
## L(Yt)        0.13845    0.06754    2.05   0.0417 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.822 on 198 degrees of freedom
## Multiple R-squared:  0.8257, Adjusted R-squared:  0.8239 
## F-statistic: 468.9 on 2 and 198 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 202
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0923 -1.1458 -0.0271  1.2168  4.5718 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.65084    2.15379   4.017 8.40e-05 ***
## Xt          -0.49371    0.03955 -12.482  < 2e-16 ***
## L(Xt)        0.31560    0.04097   7.704 6.31e-13 ***
## L(Yt)        0.60178    0.08450   7.122 1.95e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.601 on 197 degrees of freedom
## Multiple R-squared:  0.866,  Adjusted R-squared:  0.864 
## F-statistic: 424.5 on 3 and 197 DF,  p-value: < 2.2e-16

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 202
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2), data = train.ts)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.518 -1.549 -0.116  1.467  5.049 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 23.36258    0.41070  56.885   <2e-16 ***
## Xt          -0.63286    0.04047 -15.637   <2e-16 ***
## L(Xt)       -0.00559    0.06127  -0.091   0.9274    
## L(Xt, 2)     0.08827    0.04064   2.172   0.0311 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.779 on 196 degrees of freedom
## Multiple R-squared:  0.8341, Adjusted R-squared:  0.8315 
## F-statistic: 328.4 on 3 and 196 DF,  p-value: < 2.2e-16

#SSE
deviance(cons_lm1)

## [1] 635.2523

deviance(cons_lm2)

## [1] 657.3868

deviance(cons_lm3)

## [1] 505.1932

deviance(cons_lm4)

## [1] 620.3187

Uji Diagnostik Model

Uji Model

#uji model
if(require("lmtest")) encomptest(cons_lm1, cons_lm2)

## Encompassing test
## 
## Model 1: Yt ~ Xt + L(Xt)
## Model 2: Yt ~ Xt + L(Yt)
## Model E: Yt ~ Xt + L(Xt) + L(Yt)
##           Res.Df Df      F    Pr(>F)    
## M1 vs. ME    197 -1 50.717 1.954e-11 ***
## M2 vs. ME    197 -1 59.348 6.311e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan hasil uji di atas, dapat dilihat bahwa peubah lag bagi kelembaban spesifik (Yt) dan temperatur (Xt) berpengaruh signifikan.

Uji Non-Autokorelasi

H0: tidak ada autokorelasi || H1: terdapat autokorelasi

#Diagnostik
#durbin watson
dwtest(cons_lm1)

## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 1.3734, p-value = 2.659e-06
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm2)

## 
##  Durbin-Watson test
## 
## data:  cons_lm2
## DW = 1.0007, p-value = 3.071e-13
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm3)

## 
##  Durbin-Watson test
## 
## data:  cons_lm3
## DW = 1.5947, p-value = 0.001367
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm4)

## 
##  Durbin-Watson test
## 
## data:  cons_lm4
## DW = 1.3232, p-value = 5.664e-07
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan keempat model, dapat dilihat nilai p-value < 0.05 sehingga terjadi autokorelasi pada semua model autoregressive.

Uji Heterogenitas

H0: tidak ada gejala heteroskedastisitas || H1: terdapat gejala heteroskedastisitas

#heterogenitas
bptest(cons_lm1)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 3.7582, df = 2, p-value = 0.1527

bptest(cons_lm2)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 7.2931, df = 2, p-value = 0.02608

bptest(cons_lm3)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm3
## BP = 11.151, df = 3, p-value = 0.01093

bptest(cons_lm4)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm4
## BP = 5.2038, df = 3, p-value = 0.1575

Berdasarkan keempat model di atas, diketahui pada saat p = 1 dan p = 2 memiliki p-value > 0.05 sehingga tak tolak H0. Dengan demikian dapat dikatakan ragam menyebar secara homogen.

Uji Normalitas

H0: sisaan menyebar normal || H1: sisaan tidak menyebar normal

#shapiro wilk
shapiro.test(residuals(cons_lm1))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm1)
## W = 0.97847, p-value = 0.003511

shapiro.test(residuals(cons_lm2))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.9796, p-value = 0.005071

shapiro.test(residuals(cons_lm3))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm3)
## W = 0.99056, p-value = 0.2127

shapiro.test(residuals(cons_lm4))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm4)
## W = 0.97909, p-value = 0.004435

Berdasarkan keempat model di atas, diketahui pada saat p = 1 dan q = 1 memiliki p-value > 0.05 sehingga tak tolak H0. Dengan demikian dapat dikatakan bahwa sisaan menyebar normal.

Perbandingan Keakuratan Ramalan

#PERBANDINGAN
akurasi <- matrix(c(mape.koyck, mape.dlm, mape.dlm2, mape.ardl))
row.names(akurasi)<- c("Koyck","DLM 1","DLM 2","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi

##                        MAPE
## Koyck          9.632284e+33
## DLM 1          1.283605e-01
## DLM 2          1.100986e-01
## Autoregressive 1.493687e-01

Berdasarkan perbandingan nilai MAPE, dapat diketahui bahwa metode terbaik untuk memprediksi data testing adalah model dlm 2. Sebaliknya, metode KOYCK memprediksi dengan tidak tepat atau tidak mendekati data aktual. Pernyataan tersebut sesuai dengan plot perbandingan yang dihasilkan bahwa yang mendekati plot data aktual adalah model dlm 2.

#PLOT
par(mfrow=c(1,1))
plot(test$Xt, test$Yt, type="b", col="black")
points(test$Xt, fore.koyck$forecasts,col="red")
lines(test$Xt, fore.koyck$forecasts,col="red")
points(test$Xt, fore.dlm$forecasts,col="blue")
lines(test$Xt, fore.dlm$forecasts,col="blue")
points(test$Xt, fore.dlm2$forecasts,col="orange")
lines(test$Xt, fore.dlm2$forecasts,col="orange")
points(test$Xt, fore.ardl$forecasts,col="green")
lines(test$Xt, fore.ardl$forecasts,col="green")
legend("topleft",c("aktual", "koyck","DLM 1","DLM 2", "autoregressive"), lty=1, col=c("black","red","blue","orange","green"), cex=0.8)

Daftar Pustaka

Sumber : https://www.kaggle.com/datasets/poojag718/rainfall-timeseries-data

Fitriani, Herdiani ET, Saleh AF. 2013. Pemodelan autoregessive (AR) pada data hilang dan aplikasinya pada data kurs mata uang rupiah. Jurnal Matematika, Statistika, dan Komputasi. 9(2) : 69 - 85.

Lihawa SH, Resmawan, Isa DR, Nashar LO. 2022. Distributed lag model pengaruh jumlah uang beredar terhadap nilai tukar rupiah menggunakan metode koyck dan almon. Jambura Journal Probability and Statistic. 3(1): 40 - 45.