Packages

library(dLagM)
## Warning: package 'dLagM' was built under R version 4.3.3
## Loading required package: nardl
## Warning: package 'nardl' was built under R version 4.3.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.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.2
## 
## 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.3.3
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:dLagM':
## 
##     MAPE
## The following object is masked from 'package:base':
## 
##     Recall
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.2
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData

Import Data

Data ini terdiri atas angka harapan hidup sebagai peubah dependen (Y) dan supply pangan sebagai peubah independen (X) Indonesia periode 1970 - 2019.

library(readxl)
## Warning: package 'readxl' was built under R version 4.3.2
data <- read_excel("D:/Sem 5/MPDW/Tugas/Individu/3/AngkaHarapanHidup~SupplyPangan Indonesia.xlsx")
str(data)
## tibble [50 × 4] (S3: tbl_df/tbl/data.frame)
##  $ t     : num [1:50] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Yt    : num [1:50] 53 53.6 54.2 54.8 55.4 ...
##  $ Y(t-1): num [1:50] NA 53 53.6 54.2 54.8 ...
##  $ Xt    : num [1:50] 1.15e+08 1.18e+08 1.22e+08 1.25e+08 1.28e+08 ...
data
## # A tibble: 50 × 4
##        t    Yt `Y(t-1)`        Xt
##    <dbl> <dbl>    <dbl>     <dbl>
##  1     1  53.0     NA   115228390
##  2     2  53.6     53.0 118347140
##  3     3  54.2     53.6 121504136
##  4     4  54.8     54.2 124709060
##  5     5  55.4     54.8 127945190
##  6     6  56.0     55.4 131213220
##  7     7  56.5     56.0 134521020
##  8     8  57.1     56.5 137861540
##  9     9  57.6     57.1 141250960
## 10    10  58.2     57.6 144693090
## # ℹ 40 more rows

Pembagian Data

# Split Data
training <- data[1:40,]
testing <- data[41:50,]
# Data Time Series
training.ts <- ts(training)
testing.ts <- ts(testing)
data.ts <- ts(data)

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

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

koyckDlm(x , y , intercept)

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.

# Model Koyck
modelkoyck <- koyckDlm(x = training$Xt, y = training$Yt)
summary(modelkoyck)
## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92147 -0.07364  0.01528  0.11772  1.28354 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.294e+00  3.791e+00   2.188   0.0352 *  
## Y.1         8.180e-01  9.491e-02   8.619 2.81e-10 ***
## X.t         1.873e-08  1.180e-08   1.587   0.1212    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4114 on 36 degrees of freedom
## Multiple R-Squared: 0.992,   Adjusted R-squared: 0.9916 
## Wald test:  2237 on 2 and 36 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha         beta       phi
## Geometric coefficients:  45.57815 1.873248e-08 0.8180297
AIC(modelkoyck)
## [1] 46.27588
BIC(modelkoyck)
## [1] 52.93012

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\) sementara \(x_t\) tidak berpengaruh signifikan terhadap \(y_t\). Artinya, angka harapan hidup satu periode sebelumnya berpengaruh signifikan terhadap angka harapan hidup periode sekarang. Adapun model keseluruhannya adalah sebagai berikut

\[ \hat{Y_t}=8.294+0.00000001873X_t+0.818Y_{t-1} \]

Peramalan

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

forekoyck <- forecast(model = modelkoyck, x = training$Xt, h = 10)
forekoyck
## $forecasts
##  [1] 66.48218 64.89520 63.65615 62.70260 61.98319 61.45591 61.08654 60.84696
##  [9] 60.71447 60.67057
## 
## $call
## forecast.koyckDlm(model = modelkoyck, x = training$Xt, h = 10)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

Akurasi

# Akurasi Data Testing
mapekoyck <- MAPE(forekoyck$forecasts, testing$Yt)
mapekoyck
## [1] 0.1018948
# Akurasi Data Training
GoF(modelkoyck)
##             n       MAE           MPE        MAPE       sMAPE     MASE
## modelkoyck 39 0.1942202 -4.634154e-05 0.003058909 0.003053095 0.404581
##                  MSE      MRAE    GMRAE
## modelkoyck 0.1562271 0.5061085 0.250326

Regresion 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

finiteDLMauto(formula = Yt ~ Xt,
              data = data.frame(training), q.min = 1, q.max = 10,
              model.typ = "dlm", error.type = "AIC", trace = FALSE)
##    q - k    MASE      AIC      BIC   GMRAE   MBRAE R.Adj.Sq Ljung-Box
## 10    10 0.50952 37.61211 55.82768 0.41145 -0.1325  0.98336 0.2917478

Berdasarkan output tersebut, lag optimum didapatkan ketika lag = 10. Selanjutnya dilakukan pemodelan untuk lag = 10.

Pemodelan dengan Lag = 10

modeldlm <- dlm(x = training$Xt, y = training$Yt, q = 10)
summary(modeldlm)
## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.77888 -0.19534  0.01772  0.20152  0.66379 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  2.119e+01  8.122e+00   2.609   0.0178 *
## x.t          5.366e-06  2.879e-06   1.864   0.0788 .
## x.1         -1.429e-05  7.888e-06  -1.811   0.0868 .
## x.2          1.861e-05  1.157e-05   1.609   0.1250  
## x.3         -1.553e-05  1.347e-05  -1.153   0.2639  
## x.4          1.173e-05  1.436e-05   0.817   0.4248  
## x.5         -6.058e-06  1.497e-05  -0.405   0.6906  
## x.6         -2.646e-06  1.613e-05  -0.164   0.8715  
## x.7          5.389e-06  1.812e-05   0.297   0.7696  
## x.8         -8.511e-06  1.738e-05  -0.490   0.6303  
## x.9          1.300e-05  1.141e-05   1.139   0.2696  
## x.10        -6.942e-06  3.500e-06  -1.984   0.0627 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3791 on 18 degrees of freedom
## Multiple R-squared:  0.9897, Adjusted R-squared:  0.9834 
## F-statistic: 156.8 on 11 and 18 DF,  p-value: 1.563e-15
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 37.61211 55.82768
AIC(modeldlm)
## [1] 37.61211
BIC(modeldlm)
## [1] 55.82768

Dari hasil tersebut peubah \(x\) secara simultan berpengaruh terhadap peubah \(y_t\) pada taraf nyata 5%. Namun, tidak terdapat peubah yang berpengaruh signifikan secara parsial terhadap taraf nyata 5% terhadap \(y_t\). Adapun keseluruhan model yang terbentuk adalah

\[ \hat{Y_t}=21.119+0.000005366X_t-...-0.000006942X_{t-10} \]

Peramalan

Berikut merupakan hasil peramalan \(y\) untuk 10 periode kedepan.

foredlm <- forecast(model = modeldlm, x = training$Xt, h = 10)
foredlm
## $forecasts
##  [1]  -621.27843  1219.30294 -1177.75095   822.39679  -688.72014    93.04647
##  [7]   434.71199  -258.68342   838.07115  -836.62050
## 
## $call
## forecast.dlm(model = modeldlm, x = training$Xt, h = 10)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

Akurasi

# Akurasi Data Testing
mapedlm <- MAPE(foredlm$forecasts, testing$Yt)
mapedlm
## [1] 10.06573
# Akurasi Data Training
GoF(modeldlm)
##           n       MAE           MPE       MAPE       sMAPE      MASE      MSE
## modeldlm 30 0.2295306 -2.175937e-05 0.00354587 0.003545361 0.5095214 0.086224
##            MRAE     GMRAE
## modeldlm 1.0362 0.4114548

Model Autoregressive

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

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.

Lag Optimum

modelardl <- ardlBoundOrders(data = data.frame(data), ic = "AIC", 
                                  formula = Yt ~ Xt )
min_p=c()
for(i in 1:5){
  min_p[i]=min(modelardl$Stat.table[[i]])
}
q_opt=which(min_p==min(min_p, na.rm = TRUE))
p_opt=which(modelardl$Stat.table[[q_opt]] == 
              min(modelardl$Stat.table[[q_opt]], na.rm = TRUE))
data.frame("q_optimum" = q_opt, "p_optimum" = p_opt, 
           "AIC"=modelardl$min.Stat)
##   q_optimum p_optimum       AIC
## 1         1         4 -75.75334

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

Selanjutnya dapat dilakukan pemodelan dengan nilai \(p\) dan \(q\) optimum.

Pemodelan dengan Lag Optimum

modelardl <- ardlDlm(formula = Yt ~ Xt, data = training, p = 4, q = 1)
summary(modelardl)
## 
## Time series regression with "ts" data:
## Start = 5, End = 40
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24409 -0.16700 -0.01354  0.22831  0.55611 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.074e+00  4.886e+00   1.652 0.109269    
## Xt.t         6.805e-06  2.000e-06   3.402 0.001969 ** 
## Xt.1        -1.883e-05  5.960e-06  -3.159 0.003687 ** 
## Xt.2         2.038e-05  8.080e-06   2.523 0.017377 *  
## Xt.3        -9.657e-06  6.203e-06  -1.557 0.130373    
## Xt.4         1.337e-06  2.138e-06   0.625 0.536753    
## Yt.1         6.631e-01  1.493e-01   4.442 0.000119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3585 on 29 degrees of freedom
## Multiple R-squared:  0.9931, Adjusted R-squared:  0.9917 
## F-statistic: 700.2 on 6 and 29 DF,  p-value: < 2.2e-16
AIC(modelardl)
## [1] 36.52248
BIC(modelardl)
## [1] 49.19063

Dari hasil tersebut, didapat bahwa peubah \(x_t\), \(x_{t-1}\), \(x_{t-2}\), dan \(y_{t-1}\) memiliki nilai \(P-Value<0.05\). Hal ini menunjukkan bahwa peubah \(x_t\), \(x_{t-1}\), \(x_{t-2}\), dan \(y_{t-1}\) berpengaruh signifikan terhadap \(y_t\) sementara \(x_{t-3}\) dan \(x_{t-4}\) tidak berpengaruh signifikan terhadap \(y_t\). Artinya, supply pangan pada periode sekarang, satu periode sebelumnya, dua periode sebelumnya, dan angka harapan hidup satu periode sebelumnya berpengaruh signifikan terhadap angka harapan hidup periode sekarang. Adapun model keseluruhannya adalah sebagai berikut

\[ \hat{Y_t}=8.074+0.000006805X_t-...+0.6631Y_{t-1} \]

Peramalan

foreardl <- forecast(model = modelardl, x = training$Xt, h = 10)
foreardl
## $forecasts
##  [1] -807.058041 1037.532878 -365.345969  -51.190781  -15.039224    9.260458
##  [7]   25.595721   36.536743   44.162002   49.352246
## 
## $call
## forecast.ardlDlm(model = modelardl, x = training$Xt, h = 10)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

Akurasi

# Akurasi Data Testing
mapeardl <- MAPE(foreardl$forecasts, testing$Yt)
mapeardl
## [1] 3.872911
# Akurasi Data Training
GoF(modelardl)
##            n       MAE           MPE        MAPE       sMAPE      MASE
## modelardl 36 0.2320023 -2.698237e-05 0.003615361 0.003612415 0.4954591
##                 MSE      MRAE     GMRAE
## modelardl 0.1035389 0.7500387 0.3739476

Pemodelan DLM & ARDL dengan Library dynlm

Pemodelan regresi dengan peubah lag tidak hanya dapat dilakukan dengan fungsi pada packages dLagM , tetapi terdapat packages dynlm yang dapat digunakan. Fungsi dynlm secara umum adalah sebagai berikut.

dynlm(formula, data, subset, weights, na.action, method = "qr",
  model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE,
  contrasts = NULL, offset, start = NULL, end = NULL, ...)

Untuk menentukan formula model yang akan digunakan, tersedia fungsi tambahan yang memungkinkan spesifikasi dinamika (melalui d() dan L()) atau pola linier/siklus dengan mudah (melalui trend(), season(), dan harmon()). Semua fungsi formula baru mengharuskan argumennya berupa objek deret waktu (yaitu, "ts" atau "zoo").

Pemodelan

#sama dengan model dlm q=10
cons_lm1 <- dynlm(Yt ~ Xt+L(Xt)+L(Xt,2)+L(Xt,3)+L(Xt,4)+L(Xt,5)+L(Xt,6)+L(Xt,7)+L(Xt,8)+L(Xt,9)+L(Xt,10),data = training.ts)
cons_lm1
## 
## Time series regression with "ts" data:
## Start = 11, End = 40
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2) + L(Xt, 3) + L(Xt, 
##     4) + L(Xt, 5) + L(Xt, 6) + L(Xt, 7) + L(Xt, 8) + L(Xt, 9) + 
##     L(Xt, 10), data = training.ts)
## 
## Coefficients:
## (Intercept)           Xt        L(Xt)     L(Xt, 2)     L(Xt, 3)     L(Xt, 4)  
##   2.119e+01    5.366e-06   -1.429e-05    1.861e-05   -1.553e-05    1.173e-05  
##    L(Xt, 5)     L(Xt, 6)     L(Xt, 7)     L(Xt, 8)     L(Xt, 9)    L(Xt, 10)  
##  -6.058e-06   -2.646e-06    5.389e-06   -8.511e-06    1.300e-05   -6.942e-06
summary(cons_lm1)
## 
## Time series regression with "ts" data:
## Start = 11, End = 40
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2) + L(Xt, 3) + L(Xt, 
##     4) + L(Xt, 5) + L(Xt, 6) + L(Xt, 7) + L(Xt, 8) + L(Xt, 9) + 
##     L(Xt, 10), data = training.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.77888 -0.19534  0.01772  0.20152  0.66379 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  2.119e+01  8.122e+00   2.609   0.0178 *
## Xt           5.366e-06  2.879e-06   1.864   0.0788 .
## L(Xt)       -1.429e-05  7.888e-06  -1.811   0.0868 .
## L(Xt, 2)     1.861e-05  1.157e-05   1.609   0.1250  
## L(Xt, 3)    -1.553e-05  1.347e-05  -1.153   0.2639  
## L(Xt, 4)     1.173e-05  1.436e-05   0.817   0.4248  
## L(Xt, 5)    -6.058e-06  1.497e-05  -0.405   0.6906  
## L(Xt, 6)    -2.646e-06  1.613e-05  -0.164   0.8715  
## L(Xt, 7)     5.389e-06  1.812e-05   0.297   0.7696  
## L(Xt, 8)    -8.511e-06  1.738e-05  -0.490   0.6303  
## L(Xt, 9)     1.300e-05  1.141e-05   1.139   0.2696  
## L(Xt, 10)   -6.942e-06  3.500e-06  -1.984   0.0627 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3791 on 18 degrees of freedom
## Multiple R-squared:  0.9897, Adjusted R-squared:  0.9834 
## F-statistic: 156.8 on 11 and 18 DF,  p-value: 1.563e-15
#sama dengan model ardl p=4 q=1
cons_lm2 <- dynlm(Yt ~ Xt+L(Xt)+L(Xt,2)+L(Xt,3)+L(Xt,4)+L(Yt),data = training.ts)
cons_lm2
## 
## Time series regression with "ts" data:
## Start = 5, End = 40
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2) + L(Xt, 3) + L(Xt, 
##     4) + L(Yt), data = training.ts)
## 
## Coefficients:
## (Intercept)           Xt        L(Xt)     L(Xt, 2)     L(Xt, 3)     L(Xt, 4)  
##   8.074e+00    6.805e-06   -1.883e-05    2.039e-05   -9.657e-06    1.337e-06  
##       L(Yt)  
##   6.631e-01
summary(cons_lm2)
## 
## Time series regression with "ts" data:
## Start = 5, End = 40
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2) + L(Xt, 3) + L(Xt, 
##     4) + L(Yt), data = training.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24409 -0.16700 -0.01354  0.22831  0.55611 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.074e+00  4.886e+00   1.652 0.109269    
## Xt           6.805e-06  2.000e-06   3.402 0.001969 ** 
## L(Xt)       -1.883e-05  5.960e-06  -3.159 0.003687 ** 
## L(Xt, 2)     2.038e-05  8.080e-06   2.523 0.017377 *  
## L(Xt, 3)    -9.657e-06  6.203e-06  -1.557 0.130373    
## L(Xt, 4)     1.337e-06  2.138e-06   0.625 0.536753    
## L(Yt)        6.631e-01  1.493e-01   4.442 0.000119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3585 on 29 degrees of freedom
## Multiple R-squared:  0.9931, Adjusted R-squared:  0.9917 
## F-statistic: 700.2 on 6 and 29 DF,  p-value: < 2.2e-16

SSE

deviance(cons_lm1)
## [1] 2.58672
deviance(cons_lm2)
## [1] 3.727399

Uji Diagnostik

Autokorelasi

dwtest(cons_lm1)
## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 1.6142, p-value = 0.03518
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(cons_lm2)
## 
##  Durbin-Watson test
## 
## data:  cons_lm2
## DW = 2.5269, p-value = 0.8099
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil uji diagnostik model yaitu autokorelasi diperoleh \(P-Value<0.05\) untuk kedua model. Artinya, tidak terjadi autokorelasi pada model dlm dan ardl.

Heterogenitas

bptest(cons_lm1)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 18.569, df = 11, p-value = 0.06929
bptest(cons_lm2)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 11.614, df = 6, p-value = 0.07114

Berdasarkan hasil uji diagnostik model yaitu heterogenitas diperoleh \(P-Value<0.05\) untuk kedua model. Artinya, tidak terjadi heterogenitas pada model dlm dan ardl.

Normalitas

shapiro.test(residuals(cons_lm1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm1)
## W = 0.98891, p-value = 0.9845
ks.test(residuals(cons_lm1), "pnorm", mean=mean(residuals(cons_lm1)), sd=sd(residuals(cons_lm1)))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  residuals(cons_lm1)
## D = 0.074553, p-value = 0.9918
## alternative hypothesis: two-sided
shapiro.test(residuals(cons_lm2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.90139, p-value = 0.003724
ks.test(residuals(cons_lm2), "pnorm", mean=mean(residuals(cons_lm2)), sd=sd(residuals(cons_lm2)))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  residuals(cons_lm2)
## D = 0.11885, p-value = 0.646
## alternative hypothesis: two-sided

Berdasarkan hasil uji diagnostik model yaitu normalitas dengan menggunakan uji shapiro-wilk diperoleh \(P-Value<0.05\) untuk model dlm tetapi \(P-Value>0.05\) untuk model ardl. Artinya, normalitas terpenuhi pada model dlm tetapi tidak pada model ardl. Namun, berdasarkan hasil uji diagnostik model yaitu normalitas dengan menggunakan uji kolmogorov-smirnov diperoleh \(P-Value<0.05\) untuk kedua model. Artinya, normalitas terpenuhi pada model dlm dan ardl.

Perbandingan Model

akurasi <- matrix(c(mapekoyck, mapedlm, mapeardl))
row.names(akurasi)<- c("Koyck","DLM","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi
##                      MAPE
## Koyck           0.1018948
## DLM            10.0657312
## Autoregressive  3.8729112

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

Plot

par(mfrow=c(1,1))
plot(testing$Xt, testing$Yt, type="b", col="black", ylim=c(0,200))
points(testing$Xt, forekoyck$forecasts,col="red")
lines(testing$Xt, forekoyck$forecasts,col="red")
points(testing$Xt, foredlm$forecasts,col="blue")
lines(testing$Xt, foredlm$forecasts,col="blue")
points(testing$Xt, foreardl$forecasts,col="green")
lines(testing$Xt, foreardl$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 plot yang paling mendekati data aktualnya adalah Model koyck, sehingga dapat disimpulkan model terbaik dalam hal ini adalah model regresi koyck.