Regresi dengan Peubah Lag

Kali ini akan dilakukan pemodelan regresi dengan peubah lag dengan mencobakan metode-metode seperti koyck, ADRL, DLM, Autoregressive dan perbandingan terhadap model-model tersebut.

Preprocessing

Packages

Impor Data

library(rio)
data<- import("https://raw.githubusercontent.com/fax17/MPDW/main/pertemuan%203/data%20prak%203.csv")
data

##     t    yt y(t-1)    xt
## 1   1  8100     NA 11769
## 2   2 10253   8100 11654
## 3   3 10784  10253 12468
## 4   4 10973  10784 11905
## 5   5 10060  10973 10631
## 6   6 10842  10060 11366
## 7   7 11367  10842 12699
## 8   8 12961  11367 11888
## 9   9 11546  12961  9306
## 10 10 11623  11546  8394
## 11 11 11648  11623  8274
## 12 12 11918  11648  7696
## 13 13 11335  11918  6256
## 14 14 10425  11335  5636
## 15 15 10848  10425  5025
## 16 16 10627  10848  4691
## 17 17 11697  10627  5522
## 18 18 13196  11697  7201
## 19 19 14224  13196  7592
## 20 20 12374  14224  7728
## 21 21 13270  12374  7164
## 22 22 12774  13270  4165
## 23 23 12584  12774  3632
## 24 24 14614  12584  3599
## 25 25 14475  14614  2429
## 26 26 14355  14475  1941
## 27 27 16780  14355  2003
## 28 28 20623  16780  1982
## 29 29 23859  20623  2217
## 30 30 26979  23859  2420
## 31 31 30191  26979  1826
## 32 32 27835  30191  1428
## 33 33 23252  27835  1362
## 34 34 21424  23252  1389
## 35 35 21704  21424  1580
## 36 36 21032  21704  1392
## 37 37 20779  21032  1431
## 38 38 21517  20779  1810
## 39 39 23133  21517  2086
## 40 40 25070  23133  2073
## 41 41 35600  25070  2211
## 42 42 37321  35600  2498
## 43 43 39201  37321  2462
## 44 44 36039  39201  3333
## 45 45 35397  36039  3356
## 46 46 37502  35397  4171
## 47 47 36585  37502  4000
## 48 48 34392  36585  3768
## 49 49 34316  34392  1514
## 50 50 38142  34316  1621
## 51 51 40072  38142  1341
## 52 52 30972  40072  1425
## 53 53 37633  30972  1458
## 54 54 33335  37633  1442
## 55 55 38293  33335  1417
## 56 56 41091  38293  1751
## 57 57 42954  41091  1903
## 58 58 45630  42954  2109
## 59 59 41815  45630  2453
## 60 60 38301  41815  2431
## 61 61 40388  38301  2439
## 62 62 40600  40388  3079
## 63 63 29859  40600  3129
## 64 64 31592  29859  3828
## 65 65 30722  31592  4244

Pembagian Data

#SPLIT DATA
train<-data[1:50,]
test<-data[51:65,]

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

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
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 
## -5435.0 -1397.2  -313.4   995.2  9698.1 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2325.82769 1381.24485   1.684    0.099 .  
## Y.1            0.95533    0.04451  21.463   <2e-16 ***
## X.t           -0.16913    0.12288  -1.376    0.175    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2308 on 46 degrees of freedom
## Multiple R-Squared: 0.9474,  Adjusted R-squared: 0.9451 
## Wald test: 414.6 on 2 and 46 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha       beta       phi
## Geometric coefficients:  52062.85 -0.1691289 0.9553265

AIC(model.koyck)

## [1] 902.8919

BIC(model.koyck)

## [1] 910.4592

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

\[ \hat{Y_t}=2325.827+-0.169X_t+0.955Y_{t-1} \]

Peramalan dan Akurasi

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

fore.koyck <- forecast(model = model.koyck, x=test$xt, h=15)
fore.koyck

## $forecasts
##  [1] 38537.09 38900.32 39241.75 39570.63 39889.04 40136.75 40347.67 40514.34
##  [9] 40615.38 40715.62 40810.04 40791.99 40766.30 40623.53 40416.78
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$xt, h = 15)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

mape.koyck <- MAPE(fore.koyck$forecasts, test$yt)
#akurasi data training
GoF(model.koyck)

##              n      MAE          MPE       MAPE     sMAPE     MASE     MSE
## model.koyck 49 1627.013 -0.009471086 0.08360347 0.0834525 1.042303 5001373
##                 MRAE    GMRAE
## model.koyck 2.356914 1.348825

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.

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 
##  -9183  -5708  -2640   3927  16309 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 29725.4004  1924.8019  15.443   <2e-16 ***
## x.t            -0.3300     1.3425  -0.246   0.8070    
## x.1             0.8534     2.0787   0.411   0.6834    
## x.2            -2.2956     1.3058  -1.758   0.0857 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7576 on 44 degrees of freedom
## Multiple R-squared:  0.4455, Adjusted R-squared:  0.4077 
## F-statistic: 11.78 on 3 and 44 DF,  p-value: 8.552e-06
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 999.5866 1008.943

AIC(model.dlm)

## [1] 999.5866

BIC(model.dlm)

## [1] 1008.943

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

\[ \hat{Y_t}=29725.4004-0.3300X_t+0.8534X_{t-1}-2.2956X_{t-2} \]

Peramalan dan Akurasi

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

fore.dlm <- forecast(model = model.dlm, x=test$xt, h=15)
fore.dlm

## $forecasts
##  [1] 27190.78 26678.46 27382.04 27222.65 27141.49 27046.67 27338.96 26633.97
##  [9] 26347.34 26175.28 25364.17 25210.32 25721.67 24064.49 24409.00
## 
## $call
## forecast.dlm(model = model.dlm, x = test$xt, h = 15)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

mape.dlm <- MAPE(fore.dlm$forecasts, test$yt)
#akurasi data training
GoF(model.dlm)

##            n     MAE         MPE      MAPE     sMAPE     MASE      MSE     MRAE
## model.dlm 48 6153.44 -0.09953821 0.3208107 0.3006614 3.887463 52614911 14.34493
##              GMRAE
## model.dlm 5.595952

Lag Optimum

#penentuan lag optimum 
finiteDLMauto(formula = yt ~ xt,
              data = data.frame(train), q.min = 1, q.max = 23,
              model.type = "dlm", error.type = "AIC", trace = FALSE)

##    q - k    MASE      AIC      BIC   GMRAE    MBRAE R.Adj.Sq    Ljung-Box
## 23    23 0.07507 415.2207 448.9124 0.10746 -0.09488  0.99182 0.0003644872

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

#model dlm dengan lag optimum
model.dlm2 <- dlm(x = train$xt,y = train$yt , q = 23)
summary(model.dlm2)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##       1       2       3       4       5       6       7       8       9      10 
##  -34.61  114.35 -158.97  178.68 -182.50  151.95 -112.56  113.40 -167.39  188.55 
##      11      12      13      14      15      16      17      18      19      20 
##  -61.33   17.06 -130.69   85.45  -78.12  225.01 -244.04  159.81  -70.74 -118.69 
##      21      22      23      24      25      26      27 
##  123.84 -268.25  570.96  -93.73 -404.53  308.64 -111.54 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.002e+04  2.250e+03  17.790  0.00314 **
## x.t          2.922e-01  4.561e-01   0.641  0.58729   
## x.1         -2.672e-02  6.407e-01  -0.042  0.97053   
## x.2         -2.000e+00  1.058e+00  -1.890  0.19926   
## x.3          1.136e+00  6.739e-01   1.685  0.23401   
## x.4          1.838e+00  7.048e-01   2.608  0.12092   
## x.5         -2.111e+00  8.441e-01  -2.502  0.12948   
## x.6          1.024e+00  9.571e-01   1.070  0.39673   
## x.7         -2.836e-01  9.072e-01  -0.313  0.78416   
## x.8          2.668e+00  9.049e-01   2.948  0.09835 . 
## x.9         -3.212e+00  8.428e-01  -3.811  0.06248 . 
## x.10         8.287e-01  7.989e-01   1.037  0.40858   
## x.11        -5.856e-02  8.425e-01  -0.070  0.95091   
## x.12         1.910e+00  8.596e-01   2.222  0.15635   
## x.13        -1.781e+00  8.215e-01  -2.168  0.16243   
## x.14         2.718e-01  8.578e-01   0.317  0.78138   
## x.15        -1.876e+00  8.934e-01  -2.100  0.17051   
## x.16         4.388e-01  7.147e-01   0.614  0.60177   
## x.17        -1.778e+00  6.681e-01  -2.661  0.11696   
## x.18         1.376e+00  6.254e-01   2.200  0.15877   
## x.19        -2.727e+00  6.113e-01  -4.461  0.04675 * 
## x.20         1.131e+00  5.422e-01   2.086  0.17225   
## x.21        -1.843e+00  5.018e-01  -3.672  0.06681 . 
## x.22         1.639e+00  4.454e-01   3.680  0.06656 . 
## x.23         4.711e-01  2.999e-01   1.571  0.25685   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 741.5 on 2 degrees of freedom
## Multiple R-squared:  0.9994, Adjusted R-squared:  0.9918 
## F-statistic: 132.4 on 24 and 2 DF,  p-value: 0.007524
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 415.2207 448.9124

AIC(model.dlm2)

## [1] 415.2207

BIC(model.dlm2)

## [1] 448.9124

Dari hasil tersebut terdapat beberapa peubah yang berpengaruh signifikan terhadap taraf nyata 5% yaitu \(x_{t-19}\) . Adapun keseluruhan model yang terbentuk adalah

\[ \hat{Y_t}=0.0004002+0.2922X_t+...+0.4711X_{t-23} \]

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

#peramalan dan akurasi
fore.dlm2 <- forecast(model = model.dlm2, x=test$xt, h=15)
mape.dlm2<- MAPE(fore.dlm2$forecasts, test$yt)
#akurasi data training
GoF(model.dlm2)

##             n      MAE          MPE        MAPE       sMAPE       MASE      MSE
## model.dlm2 27 165.7548 -4.14146e-05 0.006230376 0.006232622 0.07506748 40730.83
##                 MRAE     GMRAE
## model.dlm2 0.3354409 0.1074563

Model tersebut merupakan model yang sangat baik dengan nilai MAPE yang kurang dari 10%.

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.

model.ardl <- ardlDlm(x = train$xt, y = train$yt, p = 1 , q = 1)
summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 2, End = 50
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5398.3  -946.8  -252.5   977.4  9578.0 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2246.25515 1340.51432   1.676    0.101    
## X.t            0.49630    0.38241   1.298    0.201    
## X.1           -0.62784    0.37823  -1.660    0.104    
## Y.1            0.95652    0.04343  22.024   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2259 on 45 degrees of freedom
## Multiple R-squared:  0.9507, Adjusted R-squared:  0.9474 
## F-statistic: 289.1 on 3 and 45 DF,  p-value: < 2.2e-16

AIC(model.ardl)

## [1] 901.701

BIC(model.ardl)

## [1] 911.1601

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

\[ \hat{Y}=2246.25515 +0.4963X_t-0.62784X_{t-1}+0.95652_{t-1} \]

Peramalan dan Akurasi

fore.ardl <- forecast(model = model.ardl, x=test$xt, h=15)
fore.ardl

## $forecasts
##  [1] 38377.59 38820.41 39207.62 39549.33 39873.83 40365.67 40701.87 41030.25
##  [9] 41385.75 41498.90 41624.91 42058.05 42095.36 42446.56 42550.10
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$xt, h = 15)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

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

mape.ardl <- MAPE(fore.ardl$forecasts, test$yt)
mape.ardl

## [1] 0.1356542

#akurasi data training
GoF(model.ardl)

##             n      MAE          MPE       MAPE      sMAPE      MASE     MSE
## model.ardl 49 1448.906 -0.009097063 0.06899382 0.06903931 0.9282035 4686060
##                MRAE    GMRAE
## model.ardl 2.069443 0.902694

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: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         1        15 941.023

Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=15\) dan \(q=1\), yaitu sebesar 941.023. 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.

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").

#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 = 50
## 
## Call:
## dynlm(formula = yt ~ xt + L(xt), data = train.ts)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -10400  -5280  -2414   3694  14688 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 29076.7763  1898.9677  15.312   <2e-16 ***
## xt              0.3635     1.2980   0.280    0.781    
## L(xt)          -2.0699     1.2646  -1.637    0.108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7668 on 46 degrees of freedom
## Multiple R-squared:  0.419,  Adjusted R-squared:  0.3938 
## F-statistic: 16.59 on 2 and 46 DF,  p-value: 3.762e-06

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 50
## 
## Call:
## dynlm(formula = yt ~ xt + L(yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5332.7 -1206.0  -289.6  1119.2  9787.4 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1761.39454 1333.03361   1.321    0.193    
## xt            -0.10924    0.11689  -0.935    0.355    
## L(yt)          0.96900    0.04358  22.234   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2302 on 46 degrees of freedom
## Multiple R-squared:  0.9477, Adjusted R-squared:  0.9454 
## F-statistic: 416.4 on 2 and 46 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 50
## 
## Call:
## dynlm(formula = yt ~ xt + L(xt) + L(yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5398.3  -946.8  -252.5   977.4  9578.0 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2246.25515 1340.51432   1.676    0.101    
## xt             0.49630    0.38241   1.298    0.201    
## L(xt)         -0.62784    0.37823  -1.660    0.104    
## L(yt)          0.95652    0.04343  22.024   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2259 on 45 degrees of freedom
## Multiple R-squared:  0.9507, Adjusted R-squared:  0.9474 
## F-statistic: 289.1 on 3 and 45 DF,  p-value: < 2.2e-16

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 50
## 
## Call:
## dynlm(formula = yt ~ xt + L(xt) + L(xt, 2), data = train.ts)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -9183  -5708  -2640   3927  16309 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 29725.4004  1924.8019  15.443   <2e-16 ***
## xt             -0.3300     1.3425  -0.246   0.8070    
## L(xt)           0.8534     2.0787   0.411   0.6834    
## L(xt, 2)       -2.2956     1.3058  -1.758   0.0857 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7576 on 44 degrees of freedom
## Multiple R-squared:  0.4455, Adjusted R-squared:  0.4077 
## F-statistic: 11.78 on 3 and 44 DF,  p-value: 8.552e-06

SSE

deviance(cons_lm1)

## [1] 2704774387

deviance(cons_lm2)

## [1] 243676749

deviance(cons_lm3)

## [1] 229616956

deviance(cons_lm4)

## [1] 2525515743

Diagnostik 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     45 -1 485.0778 <2e-16 ***
## M2 vs. ME     45 -1   2.7554 0.1039    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Autokorelasi

#durbin watson
dwtest(cons_lm1)

## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 0.15362, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm2)

## 
##  Durbin-Watson test
## 
## data:  cons_lm2
## DW = 1.3377, p-value = 0.003431
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm3)

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

dwtest(cons_lm4)

## 
##  Durbin-Watson test
## 
## data:  cons_lm4
## DW = 0.14559, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Heterogenitas

bptest(cons_lm1)

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

bptest(cons_lm2)

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

bptest(cons_lm3)

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

bptest(cons_lm4)

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

Kenormalan

shapiro.test(residuals(cons_lm1))

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

shapiro.test(residuals(cons_lm2))

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

shapiro.test(residuals(cons_lm3))

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

shapiro.test(residuals(cons_lm4))

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

Perbandingan Model

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          0.1223624
## DLM 1          0.2899022
## DLM 2          0.1797005
## Autoregressive 0.1356542

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

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)

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

Regresi dengan Peubah Lag