Regresi Peubah Lag [Individu]

Metode Peramalan Deret Waktu - Regresi Peubah lag

Tugas Individu

Muhammad Haikal Rasyadan | G1401221026 | P2

Library / Packages

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

Impor Data

data <- read_excel("/Users/user/Downloads/Documents/MPDW 💹/mpdw_regresilag_individu.xlsx")
str(data)
## tibble [72 × 4] (S3: tbl_df/tbl/data.frame)
##  $ t     : num [1:72] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Yt    : num [1:72] 39.4 39.8 40.7 41.4 42.2 ...
##  $ Y(t-1): num [1:72] NA 39.4 39.8 40.7 41.4 ...
##  $ Xt    : num [1:72] 68799030 70286312 71876439 73569145 75367369 ...
data
## # A tibble: 72 × 4
##        t    Yt `Y(t-1)`       Xt
##    <dbl> <dbl>    <dbl>    <dbl>
##  1     1  39.4     NA   68799030
##  2     2  39.8     39.4 70286312
##  3     3  40.7     39.8 71876439
##  4     4  41.4     40.7 73569145
##  5     5  42.2     41.4 75367369
##  6     6  42.9     42.2 77273504
##  7     7  43.8     42.9 79282556
##  8     8  44.5     43.8 81390843
##  9     9  45.0     44.5 83592765
## 10    10  45.9     45.0 85892868
## # ℹ 62 more rows
plot(data)

Pembagian Data

#SPLIT DATA
train<-data[1:58,]
test<-data[59:71,]
#data time series
train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(data)

Model Koyck

Pemodelan

model.koyck <- koyckDlm(x = train$Xt, y = train$Yt)
summary(model.koyck)
## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0890 -0.2127  0.1089  0.3909  3.8022 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.969e+00  2.716e+00   2.566   0.0131 *  
## Y.1         8.083e-01  8.513e-02   9.495 4.14e-13 ***
## X.t         2.886e-08  1.451e-08   1.989   0.0518 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.228 on 54 degrees of freedom
## Multiple R-Squared: 0.981,   Adjusted R-squared: 0.9802 
## Wald test:  1390 on 2 and 54 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha         beta       phi
## Geometric coefficients:  36.35482 2.885827e-08 0.8082985
AIC(model.koyck)
## [1] 190.0906
BIC(model.koyck)
## [1] 198.2628

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

\[ \hat{Y_t}=6.969+(2.886 \times 10^{-8})X_t+0.8083Y_{t-1} \]

Peramalan dan Akurasi

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

fore.koyck <- forecast(model = model.koyck, x=test$Xt, h=13)
fore.koyck
## $forecasts
##  [1] 69.01712 69.77462 70.47594 71.13414 71.75933 72.35569 72.92502 73.46953
##  [9] 73.99127 74.49147 74.97096 75.43176 75.87133
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 13)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"
mape.koyck <- MAPE(fore.koyck$forecasts, test$Yt)
mape.koyck
## [1] 0.04933724
#akurasi data training
GoF(model.koyck)
##              n       MAE           MPE       MAPE      sMAPE      MASE    MSE
## model.koyck 57 0.5564001 -0.0007164419 0.01093221 0.01076609 0.6925044 1.4286
##                  MRAE     GMRAE
## model.koyck 0.6410063 0.4102565

Nilai MAPE untuk data testing sangat bagus karena dibawah dari 10%, yaitu hanya 4.9% yang menandakan peralaman menggunakan model koyck sangat akurat.

Regression with Distributed Lag

Pemodelan (Lag=2)

model.dlm <- 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 
## -6.3133 -0.0696  0.2289  0.3415  1.6115 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.500e+01  7.891e-01  31.680  < 2e-16 ***
## x.t          2.864e-06  6.066e-07   4.722 1.81e-05 ***
## x.1         -1.528e-06  1.149e-06  -1.330   0.1893    
## x.2         -1.201e-06  5.906e-07  -2.034   0.0471 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.028 on 52 degrees of freedom
## Multiple R-squared:  0.9862, Adjusted R-squared:  0.9854 
## F-statistic:  1241 on 3 and 52 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 167.8672 177.9939
AIC(model.dlm)
## [1] 167.8672
BIC(model.dlm)
## [1] 177.9939

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}=25+(2.864 \times 10^{-6})X_t - (1.528 \times 10^{-6}) X_{t-1} - (1.201 \times 10^{-6})X_{t-2} \]

Peramalan dan Akurasi

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

fore.dlm <- forecast(model = model.dlm, x=test$Xt, h=13)
fore.dlm
## $forecasts
##  [1] 69.55704 69.86040 70.29941 70.97058 71.67528 71.97360 71.93973 71.89733
##  [9] 71.89912 71.85581 71.76491 71.78535 71.43885
## 
## $call
## forecast.dlm(model = model.dlm, x = test$Xt, h = 13)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"
mape.dlm <- MAPE(fore.dlm$forecasts, test$Yt)
mape.dlm
## [1] 0.02867477
#akurasi data training
GoF(model.dlm)
##            n       MAE           MPE        MAPE       sMAPE      MASE
## model.dlm 56 0.4853715 -0.0004364361 0.009065696 0.008891523 0.6055158
##                 MSE      MRAE     GMRAE
## model.dlm 0.9813557 0.7887411 0.4974135

Nilai MAPE untuk data testing sangat bagus karena dibawah dari 10%, yaitu hanya 2.9% yang menandakan peralaman menggunakan model Regression with Distributed Lag (Lag = 2) sangat akurat.

Lag Optimum

#penentuan lag optimum 
finiteDLMauto(formula = Yt ~ Xt,
              data = data.frame(train),
              model.type = "dlm", error.type = "AIC", trace = TRUE) 
##    q - k    MASE      AIC      BIC   GMRAE    MBRAE R.Adj.Sq Ljung-Box
## 10    10 0.63389 164.8525 189.1781 0.53215 -0.16520  0.97069 0.6672765
## 9      9 0.61994 165.6179 188.3197 0.49850  2.83680  0.97317 0.6031268
## 8      8 0.60823 166.1362 187.1684 0.46556 -3.42695  0.97551 0.5512246
## 7      7 0.60334 166.6016 185.9198 0.38338  1.36779  0.97758 0.5058382
## 6      6 0.60908 167.0271 184.5883 0.49697 -0.60304  0.97944 0.4705092
## 5      5 0.60587 167.2660 183.0284 0.50491 -2.37957  0.98119 0.4571670
## 4      4 0.59513 167.4141 181.3370 0.46761  0.25879  0.98278 0.4608448
## 3      3 0.59537 167.6049 179.6489 0.49358  0.03621  0.98419 0.4872440
## 2      2 0.60552 167.8672 177.9939 0.49741  0.35651  0.98543 0.5563575
## 1      1 0.63740 172.3209 180.4931 0.47902  0.27938  0.98554 0.8147771

Berdasarkan output tersebut, lag optimum didapatkan ketika lag=10 dikarenakan memiliki nilai AIC terendah yaitu 164.8525. Selanjutnya dilakukan pemodelan untuk lag=10

#model dlm dengan lag optimum
model.dlm2 <- dlm(x = train$Xt,y = train$Yt , q = 10)
summary(model.dlm2)
## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1615 -0.1545  0.1386  0.4360  1.2062 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.571e+01  2.144e+00  11.991 3.93e-14 ***
## x.t          2.260e-06  9.063e-07   2.494   0.0174 *  
## x.1         -1.615e-06  1.337e-06  -1.208   0.2349    
## x.2         -6.093e-07  1.340e-06  -0.455   0.6522    
## x.3          9.311e-08  1.349e-06   0.069   0.9454    
## x.4          7.505e-08  1.353e-06   0.055   0.9561    
## x.5          1.336e-07  1.353e-06   0.099   0.9219    
## x.6          7.405e-08  1.354e-06   0.055   0.9567    
## x.7          1.603e-08  1.354e-06   0.012   0.9906    
## x.8          1.559e-07  1.354e-06   0.115   0.9090    
## x.9         -1.717e-08  1.347e-06  -0.013   0.9899    
## x.10        -4.435e-07  9.355e-07  -0.474   0.6383    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.187 on 36 degrees of freedom
## Multiple R-squared:  0.9776, Adjusted R-squared:  0.9707 
## F-statistic: 142.5 on 11 and 36 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 164.8525 189.1781
AIC(model.dlm2)
## [1] 164.8525
BIC(model.dlm2)
## [1] 189.1781

Dari hasil tersebut intercept dan \(x_t\) berpengaruh signifikan terhadap taraf nyata 10%. Adapun keseluruhan model yang terbentuk adalah

\[ \hat{Y_t}=25.71+(2.260 \times 10^{-6})X_t-...-(4.435 \times 10^{-7})X_{t-10} \]

Adapun hasil peramalan 13 periode kedepan menggunakan model tersebut adalah sebagai berikut

#peramalan dan akurasi
fore.dlm2 <- forecast(model = model.dlm2, x=test$Xt, h=13)
mape.dlm2<- MAPE(fore.dlm2$forecasts, test$Yt)
mape.dlm2
## [1] 0.01756039
#akurasi data training
GoF(model.dlm2)
##             n       MAE           MPE        MAPE       sMAPE      MASE
## model.dlm2 48 0.5169128 -0.0004769541 0.009843321 0.009650366 0.6338864
##                 MSE      MRAE     GMRAE
## model.dlm2 1.056418 0.7392689 0.5321502

Nilai MAPE untuk data testing sangat bagus karena dibawah dari 10%, yaitu hanya 1.8% yang menandakan peralaman menggunakan model Regression with Optimal Distributed Lag (Lag = 10) sangat akurat.

Model Autoregressive

Pemodelan

model.ardl <- ardlDlm(formula = Yt ~ Xt, 
                         data = train,p = 1 , q = 1)
summary(model.ardl)
## 
## Time series regression with "ts" data:
## Start = 2, End = 58
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2310 -0.0082  0.1588  0.2871  1.8763 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.951e+01  3.444e+00   5.665 6.14e-07 ***
## Xt.t         3.060e-06  6.319e-07   4.843 1.15e-05 ***
## Xt.1        -2.955e-06  6.162e-07  -4.797 1.35e-05 ***
## Yt.1         2.275e-01  1.396e-01   1.630    0.109    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.035 on 53 degrees of freedom
## Multiple R-squared:  0.9867, Adjusted R-squared:  0.986 
## F-statistic:  1313 on 3 and 53 DF,  p-value: < 2.2e-16
AIC(model.ardl)
## [1] 171.5329
BIC(model.ardl)
## [1] 181.7481

Hasil di atas menunjukkan bahwa selain peubah \(y_{t-1}\), hasil uji t menunjukkan nilai-p pada peubah \(\le0.05\) Hal ini menunjukkan bahwa selain peubah \(y_{t-1}\) berpengaruh signifikan terhadap \(y_t\), sementara itu \(y_{t-1}\) tidak berpengaruh signifikan terhadap \(y_t\). Model keseluruhannya adalah sebagai berikut:

\[ \hat{Y}=19.51+(3.060 \times 10^{-6})X_t-(2.955 \times 10^{-6})X_{t-1}+0.2275Y_{t-1} \] ### Peramalan dan Akurasi

fore.ardl <- forecast(model = model.ardl, x=test$Xt, h=13)
fore.ardl
## $forecasts
##  [1] 69.37467 69.86683 70.37105 71.05251 71.73487 72.00070 71.99731 71.99849
##  [9] 72.01868 71.98661 71.91566 71.96582 71.59455
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 13)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

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

mape.ardl <- MAPE(fore.ardl$forecasts, test$Yt)
mape.ardl
## [1] 0.02973151
#akurasi data training
GoF(model.ardl)
##             n       MAE           MPE        MAPE       sMAPE      MASE
## model.ardl 57 0.4593649 -0.0004603633 0.008649055 0.008478512 0.5717329
##                 MSE      MRAE     GMRAE
## model.ardl 0.996041 0.7023964 0.3868153

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

Lag Optimum

#penentuan lag optimum
model.ardl.opt <- ardlBoundOrders(data = data.frame(data), ic = "AIC", 
                                  formula = Yt ~ Xt )
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("q_optimum" = q_opt, "p_optimum" = p_opt, 
           "AIC"=model.ardl.opt$min.Stat)
##   q_optimum p_optimum      AIC
## 1         1        15 89.29585

Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=15\) dan \(q=1\), yaitu sebesar 89.29585. Artinya, model autoregressive optimum didapat ketika \(p=15\) dan \(q=1\).

Selanjutnya dapat dilakukan pemodelan dengan nilai \(p\) dan \(q\) optimum seperti inisialisasi di langkah sebelumnya.

model.ardl.opt <- ardlDlm(formula = Yt ~ Xt, 
                         data = train,p = 15 , q = 1)
summary(model.ardl.opt)
## 
## Time series regression with "ts" data:
## Start = 16, End = 58
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3536 -0.2897  0.1245  0.6255  1.4918 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.732e+01  6.380e+00   4.282  0.00024 ***
## Xt.t         2.039e-06  1.142e-06   1.785  0.08637 .  
## Xt.1        -1.636e-06  1.497e-06  -1.093  0.28479    
## Xt.2        -6.330e-07  1.553e-06  -0.408  0.68707    
## Xt.3         1.475e-07  1.512e-06   0.098  0.92307    
## Xt.4         1.675e-07  1.515e-06   0.111  0.91279    
## Xt.5         1.802e-07  1.514e-06   0.119  0.90624    
## Xt.6         8.405e-08  1.515e-06   0.055  0.95620    
## Xt.7         7.042e-09  1.514e-06   0.005  0.99633    
## Xt.8         1.189e-07  1.515e-06   0.078  0.93806    
## Xt.9        -5.147e-09  1.517e-06  -0.003  0.99732    
## Xt.10       -4.524e-08  1.518e-06  -0.030  0.97646    
## Xt.11        1.300e-07  1.518e-06   0.086  0.93241    
## Xt.12       -1.715e-08  1.519e-06  -0.011  0.99108    
## Xt.13        1.715e-07  1.520e-06   0.113  0.91106    
## Xt.14       -8.516e-08  1.510e-06  -0.056  0.95548    
## Xt.15       -5.162e-07  1.083e-06  -0.477  0.63772    
## Yt.1        -2.712e-02  2.151e-01  -0.126  0.90067    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.327 on 25 degrees of freedom
## Multiple R-squared:  0.9715, Adjusted R-squared:  0.9522 
## F-statistic: 50.21 on 17 and 25 DF,  p-value: 3.944e-15
AIC(model.ardl.opt)
## [1] 161.0219
BIC(model.ardl.opt)
## [1] 194.4847

Peramalan dan Akurasi

fore.ardl.opt <- forecast(model = model.ardl.opt, x=test$Xt, h=13)
fore.ardl.opt
## $forecasts
##  [1] 68.30485 68.40792 68.63849 69.01141 69.33876 69.37896 69.32799 69.39512
##  [9] 69.56061 69.72182 69.88350 70.17022 70.10401
## 
## $call
## forecast.ardlDlm(model = model.ardl.opt, x = test$Xt, h = 13)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

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

mape.ardl.opt <- MAPE(fore.ardl.opt$forecasts, test$Yt)
mape.ardl.opt
## [1] 0.004361242
#akurasi data training
GoF(model.ardl.opt)
##                 n       MAE           MPE       MAPE      sMAPE      MASE
## model.ardl.opt 43 0.6035012 -0.0004840767 0.01105772 0.01090388 0.8766208
##                     MSE     MRAE     GMRAE
## model.ardl.opt 1.023367 1.008545 0.7258419

Nilai MAPE untuk data testing sangat bagus karena dibawah dari 10%, yaitu hanya 0.4% yang menandakan peralaman menggunakan model optimal Autoregressive dengan \(p=15\) dan \(q=1\) sangat akurat.

Pemodelan DLM & ARDL dengan Library dynlm

#sama dengan model dlm q=1
cons_lm1 <- dynlm(Yt ~ Xt+L(Xt),data = train.ts)
#sama dengan model ardl p=1 q=0
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 = 58
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9237 -0.0935  0.2053  0.3340  2.4151 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.500e+01  7.477e-01   33.43   <2e-16 ***
## Xt           3.954e-06  3.193e-07   12.38   <2e-16 ***
## L(Xt)       -3.817e-06  3.215e-07  -11.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.051 on 54 degrees of freedom
## Multiple R-squared:  0.9861, Adjusted R-squared:  0.9855 
## F-statistic:  1909 on 2 and 54 DF,  p-value: < 2.2e-16
summary(cons_lm2)
## 
## Time series regression with "ts" data:
## Start = 2, End = 58
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0800 -0.2217  0.1149  0.4056  3.7653 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.165e+00  2.715e+00   2.639   0.0108 *  
## Xt          2.999e-08  1.451e-08   2.067   0.0435 *  
## L(Yt)       8.018e-01  8.510e-02   9.421  5.4e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.228 on 54 degrees of freedom
## Multiple R-squared:  0.981,  Adjusted R-squared:  0.9802 
## F-statistic:  1391 on 2 and 54 DF,  p-value: < 2.2e-16
summary(cons_lm3)
## 
## Time series regression with "ts" data:
## Start = 2, End = 58
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2310 -0.0082  0.1588  0.2871  1.8763 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.951e+01  3.444e+00   5.665 6.14e-07 ***
## Xt           3.060e-06  6.319e-07   4.843 1.15e-05 ***
## L(Xt)       -2.955e-06  6.162e-07  -4.797 1.35e-05 ***
## L(Yt)        2.275e-01  1.396e-01   1.630    0.109    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.035 on 53 degrees of freedom
## Multiple R-squared:  0.9867, Adjusted R-squared:  0.986 
## F-statistic:  1313 on 3 and 53 DF,  p-value: < 2.2e-16
summary(cons_lm4)
## 
## Time series regression with "ts" data:
## Start = 3, End = 58
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3133 -0.0696  0.2289  0.3415  1.6115 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.500e+01  7.891e-01  31.680  < 2e-16 ***
## Xt           2.864e-06  6.066e-07   4.722 1.81e-05 ***
## L(Xt)       -1.528e-06  1.149e-06  -1.330   0.1893    
## L(Xt, 2)    -1.201e-06  5.906e-07  -2.034   0.0471 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.028 on 52 degrees of freedom
## Multiple R-squared:  0.9862, Adjusted R-squared:  0.9854 
## F-statistic:  1241 on 3 and 52 DF,  p-value: < 2.2e-16

SSE

#library 'dynlm'
deviance(cons_lm1)
## [1] 59.62033
deviance(cons_lm2)
## [1] 81.42099
deviance(cons_lm3)
## [1] 56.77434
deviance(cons_lm4)
## [1] 54.95592
#Manual
SSE.ardl <- sum(residuals(model.ardl)^2)
## Time Series:
## Start = 2 
## End = 58 
## Frequency = 1 
##            2            3            4            5            6            7 
## -0.462941620 -0.113574442 -0.077340528  0.031158125  0.069217607  0.262573802 
##            8            9           10           11           12           13 
##  0.251171421  0.130728882  0.287133290  0.138919407  0.135148674  0.202790699 
##           14           15           16           17           18           19 
##  0.040020163  0.107788695 -6.230956447  0.930171191  1.331314600  0.374349809 
##           20           21           22           23           24           25 
##  0.222626910  0.152774979  0.151668481  0.260027307  0.249155809  0.294832398 
##           26           27           28           29           30           31 
##  0.291222042  0.243864112  0.245318517  0.064506541 -0.006123204 -0.008235985 
##           32           33           34           35           36           37 
## -0.188035627 -0.126474369  1.876305587 -2.293542143  0.060517186  0.359841050 
##           38           39           40           41           42           43 
##  0.351998060  0.158804117  0.230961923  0.398359587  0.366552471  0.566494667 
##           44           45           46           47           48           49 
##  0.573948609  0.305286435  0.253715026 -0.122849917 -0.166725083 -0.247763190 
##           50           51           52           53           54           55 
## -0.041370407  0.099734429  0.215381525  0.290941828  0.265753656 -1.523098246 
##           56           57           58 
##  0.248820282 -0.640161057 -0.842707633
SSE.ardl.opt <- sum(residuals(model.ardl.opt)^2)
## Time Series:
## Start = 16 
## End = 58 
## Frequency = 1 
##          16          17          18          19          20          21 
## -5.35356046 -0.23855939  1.49175993  0.91194699  0.70064417  0.65480317 
##          22          23          24          25          26          27 
##  0.67612059  0.74800973  0.71788860  0.69812973  0.59623611  0.48749851 
##          28          29          30          31          32          33 
##  0.44378867  0.31208570  0.29078072  0.13179744 -0.28145912 -0.36148792 
##          34          35          36          37          38          39 
##  0.92526812 -1.65977995 -0.72801601  0.01528953 -0.06500614 -0.34873686 
##          40          41          42          43          44          45 
## -0.46678084 -0.41601380 -0.34821210 -0.22512750 -0.13256764 -0.12124360 
##          46          47          48          49          50          51 
## -0.16342492 -0.22468298 -0.29784189 -0.47350278  0.20629517  0.44082731 
##          52          53          54          55          56          57 
##  0.58350413  0.67326935  0.70163735 -1.06927221  0.39412874  0.12451173 
##          58 
##  0.04905463
SSE.dlm <- sum(residuals(model.dlm)^2)
##           1           2           3           4           5           6 
## -0.12154884 -0.02993216  0.09005836  0.14484050  0.33860404  0.36213773 
##           7           8           9          10          11          12 
##  0.23765517  0.35955677  0.22976852  0.18344702  0.23490235  0.07700522 
##          13          14          15          16          17          18 
##  0.10358928 -6.31326977 -0.45553293  1.33189657  0.61512793  0.30433884 
##          19          20          21          22          23          24 
##  0.14874744  0.11932376  0.22804249  0.24536503  0.28457790  0.29811507 
##          25          26          27          28          29          30 
##  0.25647096  0.23929274  0.06273256 -0.05960709 -0.09977444 -0.28502334 
##          31          32          33          34          35          36 
## -0.25799075  1.61151170 -1.47255264 -0.94163962  0.27324245  0.39743757 
##          37          38          39          40          41          42 
##  0.23864221  0.25616207  0.42154991  0.45367158  0.66132232  0.71846504 
##          43          44          45          46          47          48 
##  0.46953975  0.35007553 -0.04575842 -0.18304981 -0.30420496 -0.11688369 
##          49          50          51          52          53          54 
##  0.10304894  0.29071730  0.40771253  0.40805925 -1.38426828  0.03364825 
##          55          56 
## -0.56450630 -0.95485959
SSE.dlm2 <- sum(residuals(model.dlm2)^2)
##           1           2           3           4           5           6 
##  0.70535556  0.63397985  0.65346504  0.46385062  0.47279902 -6.16148113 
##           7           8           9          10          11          12 
## -0.67268316  1.06330057  0.51996918  0.34194987  0.29959517  0.34458493 
##          13          14          15          16          17          18 
##  0.47392765  0.53603421  0.62896762  0.49638591  0.23565403  0.11173565 
##          19          20          21          22          23          24 
## -0.09359347 -0.17192382 -0.16038183 -0.32157945 -0.27576591  1.20618091 
##          25          26          27          28          29          30 
## -1.52782173 -0.60891735  0.10096682  0.04136317 -0.15253318 -0.17215558 
##          31          32          33          34          35          36 
## -0.04064040  0.10453054  0.26242046  0.14149387  0.35742774  0.39211153 
##          37          38          39          40          41          42 
##  0.09787274  0.04137402 -0.02447891  0.13569290  0.25424678  0.35563610 
##          43          44          45          46          47          48 
##  0.42665224  0.41043554 -1.39362377  0.09594612 -0.25809181 -0.37023485
SSE.dlm
## [1] 54.95592
SSE.dlm2
## [1] 50.70804
SSE.ardl
## [1] 56.77434
SSE.ardl.opt
## [1] 44.00479

Nilai SSE pada cara manual lebih kecil daripada nilai SSE yang didapatkan ketika menggunakan library dynlm, sehingga lebih direkomendasikan menggunakan model dengan cara manual

Uji Diagnostik

#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     53 -1  2.6568     0.109    
## M2 vs. ME     53 -1 23.0082 1.352e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Autokorelasi

  • \(H_0\): Tidak terdapat Autokorelasi
  • \(H_1\): Terdapat Autokorelasi
#durbin watson
dwtest(cons_lm1)
## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 2.0357, p-value = 0.4494
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(cons_lm2)
## 
##  Durbin-Watson test
## 
## data:  cons_lm2
## DW = 2.391, p-value = 0.8942
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(cons_lm3)
## 
##  Durbin-Watson test
## 
## data:  cons_lm3
## DW = 2.2466, p-value = 0.7527
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(cons_lm4)
## 
##  Durbin-Watson test
## 
## data:  cons_lm4
## DW = 1.83, p-value = 0.1991
## alternative hypothesis: true autocorrelation is greater than 0

p-value pada semua model lebih besar daripada nilai \(\alpha = 0.05\), sehingga hipotesis 0 diterima atau tidak terdapat autokorelasi pada data.

Heterogenitas

  • \(H_0\): Tidak terjadi heteroskedastisitas
  • \(H_1\): Terjadi heteroskedastisitas
bptest(cons_lm1)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 0.46432, df = 2, p-value = 0.7928
bptest(cons_lm2)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 1.4182, df = 2, p-value = 0.4921
bptest(cons_lm3)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm3
## BP = 1.6523, df = 3, p-value = 0.6476
bptest(cons_lm4)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm4
## BP = 2.2032, df = 3, p-value = 0.5313

p-value pada semua model lebih besar daripada nilai \(\alpha = 0.05\), sehingga hipotesis 0 diterima atau tidak terjadi heteroskedastisitas pada data.

Kenormalan

  • \(H_0\): Data berdistribusi normal
  • \(H_1\): Data tidak berdistribusi normal
shapiro.test(residuals(cons_lm1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm1)
## W = 0.58626, p-value = 2.027e-11
shapiro.test(residuals(cons_lm2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.57752, p-value = 1.496e-11
shapiro.test(residuals(cons_lm3))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm3)
## W = 0.52621, p-value = 2.725e-12
shapiro.test(residuals(cons_lm4))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm4)
## W = 0.54273, p-value = 6.019e-12

p-value pada semua model lebih kecil daripada nilai \(\alpha = 0.05\), sehingga hipotesis 0 ditolak atau data tidak berdistribusi secara normal.

Perbandingan Model

akurasi <- matrix(c(mape.koyck, mape.dlm, mape.dlm2, mape.ardl, mape.ardl.opt))
row.names(akurasi)<- c("Koyck","DLM 1","DLM 2","ARDL 1", "ARDL 2")
colnames(akurasi) <- c("MAPE")
akurasi
##               MAPE
## Koyck  0.049337236
## DLM 1  0.028674773
## DLM 2  0.017560387
## ARDL 1 0.029731506
## ARDL 2 0.004361242

Berdasarkan nilai MAPE, model paling optimum didapat pada Model ARDL 2 karena memiliki nilai MAPE yang terkecil.

Plot

par(mfrow=c(1,1))
plot(test$Xt, test$Yt, type="b", col="black",ylim=c(68, 76))
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")
points(test$Xt, fore.ardl.opt$forecasts,col="pink")
lines(test$Xt, fore.ardl.opt$forecasts,col="pink")
legend("topleft",c("aktual", "koyck","DLM 1","DLM 2", "ARDL1", "ARDL 2"), lty=1, col=c("black","red","blue","orange","green","pink"), cex=0.8)

Berdasarkan plot tersebut, terlihat bahwa plot yang paling mendekati data aktualnya adalah Model ARDL 2, sehingga dapat disimpulkan model terbaik dalam hal ini adalah model Optimal Autoregressive dengan \(p=15\) dan \(q=1\).