Packages

library(dLagM)
## Warning: package 'dLagM' was built under R version 4.2.3
## Loading required package: nardl
## Warning: package 'nardl' was built under R version 4.2.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Loading required package: dynlm
## Warning: package 'dynlm' was built under R version 4.2.3
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(dynlm)
library(MLmetrics)
## Warning: package 'MLmetrics' was built under R version 4.2.3
## 
## 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

Impor Data

# Yt = TS [Earth Skin Temperature (C)]
# Xt = T2M [Temperature at 2 Meters (C)]

library(rio)
## Warning: package 'rio' was built under R version 4.2.3
datatug <- Import("https://raw.githubusercontent.com/afhkaniase/praktikum-mpdw/main/Pertemuan%203/Data%20Latihan%203.csv")

t <- datatug$Date
Yt <- datatug$TS
Xt <- datatug$T2M

data3 <- data.frame(Yt,Xt)
str(data3)
## 'data.frame':    123 obs. of  2 variables:
##  $ Yt: num  2.95 3.17 3.3 3.02 3.08 3.77 3.73 3 3.07 3.37 ...
##  $ Xt: num  4.53 7.36 5.97 2.85 2.55 4.96 3.64 0.44 1.97 2.8 ...
data3
##        Yt    Xt
## 1    2.95  4.53
## 2    3.17  7.36
## 3    3.30  5.97
## 4    3.02  2.85
## 5    3.08  2.55
## 6    3.77  4.96
## 7    3.73  3.64
## 8    3.00  0.44
## 9    3.07  1.97
## 10   3.37  2.80
## 11   3.53  2.43
## 12   3.67  3.18
## 13   3.90  3.58
## 14   4.02  5.53
## 15   4.58  6.40
## 16   4.70  5.44
## 17   4.53  6.01
## 18   4.56  6.88
## 19   4.64  6.82
## 20   5.00  7.49
## 21   5.41  7.91
## 22   6.05  9.55
## 23   6.22  7.40
## 24   6.72 10.98
## 25   6.93 11.74
## 26   7.21 12.94
## 27   7.14 12.00
## 28   7.07 11.75
## 29   7.22  9.82
## 30   7.32  8.29
## 31   6.65  7.81
## 32   6.51  9.73
## 33   7.05 12.24
## 34   7.14 12.13
## 35   7.15  9.80
## 36   7.42 10.91
## 37   8.05 10.32
## 38   8.20  8.83
## 39   8.72 10.96
## 40   9.15 12.66
## 41   9.57 14.78
## 42   8.87 11.76
## 43   8.07  8.33
## 44   7.88  8.66
## 45   8.25  9.36
## 46   8.92 10.63
## 47   9.88 12.51
## 48  10.75 13.58
## 49  11.71 15.44
## 50  12.38 16.07
## 51  13.31 16.65
## 52  13.31 16.65
## 53  13.91 16.36
## 54  14.64 16.32
## 55  14.31 14.86
## 56  14.21 15.84
## 57  14.06 16.58
## 58  13.82 16.44
## 59  13.73 15.37
## 60  14.05 16.02
## 61  14.40 16.00
## 62  15.86 17.33
## 63  17.37 20.02
## 64  17.67 19.36
## 65  17.66 18.79
## 66  18.33 19.95
## 67  19.41 21.21
## 68  20.32 22.07
## 69  20.99 21.59
## 70  21.87 22.79
## 71  21.80 22.20
## 72  21.00 21.03
## 73  19.69 19.35
## 74  18.46 19.03
## 75  18.26 19.15
## 76  18.48 20.34
## 77  18.51 19.15
## 78  19.42 20.95
## 79  19.70 21.63
## 80  19.67 21.71
## 81  19.16 19.08
## 82  18.95 18.04
## 83  19.19 19.61
## 84  19.30 19.10
## 85  19.80 19.86
## 86  20.57 21.54
## 87  20.99 22.63
## 88  20.92 22.07
## 89  20.73 21.37
## 90  20.80 20.80
## 91  20.85 20.04
## 92  20.94 20.65
## 93  21.10 21.20
## 94  20.89 20.53
## 95  20.44 18.05
## 96  20.15 17.97
## 97  19.75 18.46
## 98  19.92 19.90
## 99  20.28 20.34
## 100 20.65 21.14
## 101 21.16 22.05
## 102 21.73 22.79
## 103 21.37 21.13
## 104 21.44 21.54
## 105 21.72 21.63
## 106 22.15 22.50
## 107 22.11 22.29
## 108 21.87 21.39
## 109 21.76 21.01
## 110 21.26 19.19
## 111 21.06 19.30
## 112 21.24 21.77
## 113 21.25 21.00
## 114 21.30 21.68
## 115 21.74 22.34
## 116 21.85 22.57
## 117 21.71 21.35
## 118 21.46 20.94
## 119 21.39 21.62
## 120 20.78 19.52
## 121 20.61 19.40
## 122 20.13 17.69
## 123 20.05 20.02

Pembagian Data

#SPLIT DATA
data3_train<-data3[1:98,]
data3_test<-data3[99:123,]
#data time series
train.ts<-ts(data3_train)
test.ts<-ts(data3_test)
data.ts<-ts(data3)

Model Koyck

Model Koyck didasarkan pada asumsi bahwa semakin jauh jarak lag peubah independen dari periode sekarang maka semakin kecil pengaruh peubah lag terhadap peubah dependen.

Koyck mengusulkan suatu metode untuk menduga model dinamis distributed lag dengan mengasumsikan bahwa semua koefisien \(\beta\) mempunyai tanda sama.

Model kyock merupakan jenis paling umum dari model infinite distributed lag dan juga dikenal sebagai geometric lag

\[ y_t=a(1-\lambda)+\beta_0X_t+\beta_1Z_t+\lambda Y_{t-1}+V_t \]

dengan \[V_t=u_t-\lambda u_{t-1}\]

Pemodelan Model Koyck

Pemodelan model Koyck dengan R dapat menggunakan dLagM::koyckDlm() . Fungsi umum dari koyckDlm adalah sebagai berikut.

Fungsi koyckDlm() akan menerapkan model lag terdistribusi dengan transformasi Koyck satu prediktor. Nilai x dan y tidak perlu sebagai objek time series (ts). intercept dapat dibuat TRUE untuk memasukkan intersep ke dalam model.

data3_model.koyck <- koyckDlm(x = data3_train$Xt, y = data3_train$Yt)
summary(data3_model.koyck)
## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.306328 -0.256056  0.002453  0.208640  1.226161 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.10931    0.13987   0.782    0.436    
## Y.1          0.94879    0.03327  28.518   <2e-16 ***
## X.t          0.04974    0.03587   1.387    0.169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4528 on 94 degrees of freedom
## Multiple R-Squared: 0.9953,  Adjusted R-squared: 0.9952 
## Wald test: 1.002e+04 on 2 and 94 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha       beta       phi
## Geometric coefficients:  2.134555 0.04973701 0.9487908
AIC(data3_model.koyck)
## [1] 126.5151
BIC(data3_model.koyck)
## [1] 136.814

Dari hasil tersebut, didapat bahwa peubah \(y_{t-1}\) memiliki nilai \(P-Value<0.05\). Hal ini menunjukkan bahwa peubah \(y_{t-1}\) berpengaruh signifikan terhadap \(y\). Adapun model keseluruhannya adalah sebagai berikut

\[ \hat{Y_t} = 0.10931 + 0.04974X_t + 0.94879Y_{t-1} \]

Peramalan dan Akurasi Model Koyck

Berikut adalah hasil peramalan y untuk 25 periode kedepan menggunakan model koyck

fore.koyck <- forecast(model = data3_model.koyck, x=data3_test$Xt, h=25) 
fore.koyck
## $forecasts
##  [1] 20.02087 20.15637 20.33019 20.53191 20.64074 20.76439 20.88618 21.04501
##  [9] 21.18526 21.27357 21.33845 21.30949 21.28748 21.38945 21.44790 21.53718
## [17] 21.65471 21.77766 21.83364 21.86636 21.93123 21.88832 21.84165 21.71231
## [25] 21.70548
## 
## $call
## forecast.koyckDlm(model = data3_model.koyck, x = data3_test$Xt, 
##     h = 25)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"
#akurasi data testing
mape.koyck.test <- MAPE(fore.koyck$forecasts, data3_test$Yt)
mape.koyck.test
## [1] 0.02981854
#akurasi data training
mape.koyck.train <- GoF(data3_model.koyck)
mape.koyck.train
##                    n       MAE          MPE       MAPE     sMAPE      MASE
## data3_model.koyck 97 0.3315851 -0.003443496 0.03339684 0.0331023 0.8322136
##                         MSE     MRAE     GMRAE
## data3_model.koyck 0.1986777 26649597 0.9770295

Berdasarkan akurasi di atas, terlihat bahwa nilai MAPE keduanya tidak jauh berbeda. Artinya, model regresi dengan distribusi lag ini tidak  overfitted atau underfitted.

Regression with Distributed Lag

Pemodelan model Regression with Distributed Lag dengan R dapat menggunakan dLagM::dlm() . Fungsi umum dari dlm adalah sebagai berikut.

dlm(formula , data , x , y , q , remove )

Fungsi dlm() akan menerapkan model lag terdistribusi dengan satu atau lebih prediktor. Nilai x dan y tidak perlu sebagai objek time series (ts). \(q\) adalah integer yang mewakili panjang lag yang terbatas.

Lag Optimum Distribution Lag Model

#penentuan lag optimum 
finiteDLMauto(formula = Yt ~ Xt,
              data = data.frame(data3_train),
              model.type = "dlm", error.type = "AIC")
##    q - k    MASE      AIC      BIC   GMRAE   MBRAE R.Adj.Sq    Ljung-Box
## 10    10 2.25323 295.2537 327.4591 3.20525 0.52278  0.96212 2.220446e-16

Diperoleh lag optimum untuk peubah Xt = T2M [Temperature at 2 Meters (C)] adalah 10 hari sebelumnya. Selanjutnya dilakukan pemodelan kembali dengan \(q=10\)

Pemodelan (Lag=10) Distribution Lag Model

model.dlm <- dlm(x = data3_train$Xt,y = data3_train$Yt , q = 10)
summary(model.dlm)
## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.54205 -0.77300  0.09346  0.73706  2.58880 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.31919    0.37736  -6.146 3.39e-08 ***
## x.t          0.49601    0.10758   4.610 1.59e-05 ***
## x.1          0.10451    0.15465   0.676   0.5012    
## x.2          0.10045    0.15370   0.654   0.5154    
## x.3          0.02629    0.14969   0.176   0.8611    
## x.4          0.06991    0.14891   0.469   0.6401    
## x.5          0.11779    0.14813   0.795   0.4290    
## x.6          0.03569    0.14565   0.245   0.8071    
## x.7         -0.01543    0.14543  -0.106   0.9158    
## x.8         -0.04417    0.14569  -0.303   0.7626    
## x.9          0.02276    0.14328   0.159   0.8742    
## x.10         0.16909    0.09796   1.726   0.0884 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.202 on 76 degrees of freedom
## Multiple R-squared:  0.9669, Adjusted R-squared:  0.9621 
## F-statistic: 201.9 on 11 and 76 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 295.2537 327.4591
AIC(model.dlm)
## [1] 295.2537
BIC(model.dlm)
## [1] 327.4591

Dari hasil diatas, didapat bahwa \(P-value\) dari intercept dan \(x_{t}<0.05\). Hal ini menunjukkan bahwa intercept dan \(x_{t}\) berpengaruh signifikan terhadap \(y\). Adapun model keseluruhan yang terbentuk adalah sebagai berikut

\[ \hat{Y_t} = - 2.31919 + 0.49601X_t + 0.10451X_{t-1} + 0.10045X_{t-2} + 0.02629X_{t-3} + 0.06991X_{t-4} + 0.11779X_{t-5} + 0.03569X_{t-6} - 0.01543X_{t-7} - 0.04417X_{t-8} + 0.02276X_{t-9} + 0.16909X_{t-10} \]

Peramalan dan Akurasi Distribution Lag Model

Berikut merupakan hasil peramalan \(y\) untuk 25 periode kedepan

fore.dlm <- forecast(model = model.dlm, x=data3_test$Xt, h=25)
fore.dlm
## $forecasts
##  [1] 19.49596 19.62555 20.05032 20.94346 20.71624 20.82975 20.47847 21.03166
##  [9] 21.12570 20.81076 20.56632 19.76291 19.89415 20.94923 20.35982 20.69972
## [17] 20.92367 21.46341 21.19640 20.86884 20.99687 19.66869 19.51093 18.65288
## [25] 19.37090
## 
## $call
## forecast.dlm(model = model.dlm, x = data3_test$Xt, h = 25)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"
#akurasi data testing
mape.dlm.test <- MAPE(fore.dlm$forecasts, data3_test$Yt)
mape.dlm.test
## [1] 0.04159384
#akurasi data training
mape.dlm.train <- GoF(model.dlm)
mape.dlm.train
##            n       MAE         MPE      MAPE     sMAPE     MASE      MSE
## model.dlm 88 0.9269329 0.009663972 0.1142461 0.1263623 2.253232 1.248352
##               MRAE    GMRAE
## model.dlm 96906248 3.205253

Berdasarkan akurasi di atas, terlihat bahwa nilai MAPE antara data training dengan data testing didapatkan jauh berbeda, artinya, model regresi dengan distribution lag ini  overfitted .

Model Autoregressive Distributed Lag (ARDL)

Peubah dependen dipengaruhi oleh peubah independen pada waktu sekarang, serta dipengaruhi juga oleh peubah dependen itu sendiri pada satu waktu yang lalu maka model tersebut disebut autoregressive (Gujarati 2004).

Pemodelan Autoregressive Distributed Lag

Pemodelan Autoregressive dilakukan menggunakan fungsi dLagM::ardlDlm() . Fungsi tersebut akan menerapkan autoregressive berordo \((p,q)\) dengan satu prediktor. Fungsi umum dari ardlDlm() adalah sebagai berikut.

ardlDlm(formula = NULL , data = NULL , x = NULL , y = NULL , p = 1 , q = 1 , 
         remove = NULL )

Dengan \(p\) adalah integer yang mewakili panjang lag yang terbatas dan \(q\) adalah integer yang merepresentasikan ordo dari proses autoregressive.

Penentuan Lag Optimum untuk Autoregressive Distributed Lag

#penentuan lag optimum
model.ardl.opt <- ardlBoundOrders(data = data.frame(data3), ic = "AIC", 
                                  formula = Yt ~ Xt )
min_p=c()
for(i in 1:15){
  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("q_optimum" = q_opt, "p_optimum" = p_opt, 
           "AIC"=model.ardl.opt$min.Stat)
##   q_optimum p_optimum      AIC
## 1         2         1 30.82842

Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=1\) dan \(q=2\), yaitu sebesar 30.82842. Artinya, model autoregressive optimum didapat ketika \(p=1\) dan \(q=2\). Selanjutnya nilai ini akan dimasukkan ke dalam proses pembentukan model ardl.

model.ardl <- ardlDlm(x = data3_train$Xt, y = data3_train$Yt, p = 1 , q = 2)
summary(model.ardl)
## 
## Time series regression with "ts" data:
## Start = 3, End = 98
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66789 -0.20311  0.03346  0.17545  0.94769 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.03384    0.08123   0.417    0.678    
## X.t          0.21466    0.02292   9.367 5.37e-15 ***
## X.1         -0.18353    0.02535  -7.239 1.40e-10 ***
## Y.1          1.42815    0.07840  18.215  < 2e-16 ***
## Y.2         -0.46053    0.06785  -6.787 1.14e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.295 on 91 degrees of freedom
## Multiple R-squared:  0.998,  Adjusted R-squared:  0.998 
## F-statistic: 1.16e+04 on 4 and 91 DF,  p-value: < 2.2e-16
AIC(model.ardl)
## [1] 44.89524
BIC(model.ardl)
## [1] 60.28133

Hasil di atas menunjukkan bahwa terdapat 4 peubah yang berpengaruh signifikan terhadap nilai T2M [Temperature at 2 Meters (C)] pada selang kepercayaan 95% yaitu peubah \(x_t\), \(x_{t-1}\), \(y_{t-1}\), \(y_{t-2}\). Artinya, menurut model ARDL dengan \(p=1\) dan \(q=2\), nilai T2M [Temperature at 2 Meters (C)] saat ini dipengaruhi oleh TS [Earth Skin Temperature (C)] pada saat ini, serta 8 hari sebelumnya. Model ini sangat baik dengan nilai R-Square sebesar 99.8%. Model keseluruhannya adalah sebagai berikut:

\[ \hat{Y} = 0.03384 + 0.21466X_t - 0.18353X_{t-1} + 1.42815Y_{t-1} - 0.46053Y_{t-2} \]

Peramalan dan Akurasi Model Autoregressive Distributed Lag

fore.ardl <- forecast(model = model.ardl, x=data3_test$Xt, h=25)
fore.ardl
## $forecasts
##  [1] 20.10104 20.37227 20.72478 21.09513 20.96957 21.01235 21.07536 21.31587
##  [9] 21.42559 21.31687 21.19468 20.74931 20.52716 20.92502 20.97693 21.15512
## [17] 21.40257 21.60215 21.46913 21.32313 21.39710 20.99440 20.74486 20.22891
## [25] 20.42095
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = data3_test$Xt, h = 25)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

Data di atas merupakan hasil peramalan untuk 25 periode ke depan menggunakan Model Autoregressive dengan \(p=1\) dan \(q=2\).

#akurasi data testing
mape.ardl.test <- MAPE(fore.ardl$forecasts, data3_test$Yt)
mape.ardl.test
## [1] 0.0171622
#akurasi data training
mape.ardl.train <- GoF(model.ardl)
mape.ardl.train
##             n       MAE          MPE       MAPE      sMAPE      MASE        MSE
## model.ardl 96 0.2300173 -0.000136184 0.02438607 0.02433235 0.5732331 0.08247844
##                MRAE     GMRAE
## model.ardl 57837119 0.7241981

Berdasarkan akurasi di atas, terlihat bahwa nilai MAPE keduanya jauh berbeda. Artinya, model regresi dengan distribusi lag ini  underfitted.

Perbandingan Model

akurasi <- matrix(c(mape.koyck.test, mape.dlm.test, mape.ardl.test))
row.names(akurasi)<- c("Koyck","DLM","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi
##                      MAPE
## Koyck          0.02981854
## DLM            0.04159384
## Autoregressive 0.01716220

Berdasarkan nilai MAPE, model paling optimum didapat pada Model Koyck karena memiliki nilai MAPE yang terkecil dibandingkan Model DLM dan Model ARDL.

Plot

par(mfrow=c(1,1))
plot(data3_test$Xt, data3_test$Yt, type="b", col="black")
points(data3_test$Xt, fore.koyck$forecasts,col="red")
lines(data3_test$Xt, fore.koyck$forecasts,col="red")
points(data3_test$Xt, fore.dlm$forecasts,col="blue")
lines(data3_test$Xt, fore.dlm$forecasts,col="blue")
points(data3_test$Xt, fore.ardl$forecasts,col="green")
lines(data3_test$Xt, fore.ardl$forecasts,col="green")
legend("topleft",c("Aktual", "Koyck","DLM", "Autoregressive"), lty=1, col=c("black","red","blue","green"), cex=0.8)

Berdasarkan plot tersebut, terlihat bahwa model Autoregressive Distributed Lag merupakan metode yang paling sesuai untuk melakukan peramalan pada data suhu di US karena memiliki pola data aktual. Sehingga dapat disimpulkan model terbaik dalam hal ini adalah model regresi Autoregressive Distributed Lag.