Import Data

data <- read_excel("C:/SEMESTER 5/MPDW/gold.xlsx")
knitr::kable(data,align = "c")

str(data)

## tibble [198 x 7] (S3: tbl_df/tbl/data.frame)
##  $ Date    : chr [1:198] "2000-01-04" "2000-01-05" "2000-01-06" "2000-01-07" ...
##  $ Close   : num [1:198] 284 282 282 283 283 ...
##  $ High    : num [1:198] 290 285 283 284 284 ...
##  $ Low     : num [1:198] 280 281 280 282 282 ...
##  $ Open    : num [1:198] 290 284 282 282 282 ...
##  $ Volume  : num [1:198] 21621 25448 19055 11266 30603 ...
##  $ Currency: chr [1:198] "USD" "USD" "USD" "USD" ...

data1 <- data[,c(2,5)]
data.y<-ts(data$Open)
plot(data.y, main = "Time Series Plot Harga Emas (Harian)", xlab = "Periode", ylab="Open")
points(data.y)

Grafik di atas merupakan time series plot dari data Harga Emas Harian selama 10 bulan pada tahun 2000. Berdasarkan grafik tersebut, dapat dilihat bahwa data harga emas harian tersebut cenderung memiliki pola gabungan.

Split Data

Data dibagi menjadi dua bagian dengan proporsi 150 data sebagai data train dan 48 data sebagai data test.

train<-data1[1:150,]
test<-data1[151:198,]

Data Time Series

train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(data)

Hipotesis:

H0: β_i=0

H1: β_i≠0

Berdasarkan output di atas, hanya nilai p-value Xt < 0.05 yang mengindikasikan bahwa penduka parameter ini signifikan. Dari pengujian hipotesis, Xt berhasil menolak H0 sehingga Xt berpengaruh, sedangkan b0 dan Yt-1 tidak berhasil menolak H0. Menurut model koyck, nilai Close dipengaruhi oleh nilai Open namun tidak dipengaruhi oleh peubah lag-nya.

#ramalan
(fore.koyck <- dLagM::forecast(model = model.koyck, x=test$Close, h=48))

## $forecasts
##  [1] 276.6745 277.2587 280.3314 279.2641 279.3537 282.4680 282.5533 281.5056
##  [9] 279.6280 277.8200 274.9656 277.2865 278.3635 278.1575 277.0185 277.9235
## [17] 282.3568 280.5971 279.3574 277.3803 277.3150 276.2799 276.3668 276.1283
## [25] 276.0327 275.4626 275.2792 274.2540 274.9182 271.5094 272.9518 274.6807
## [33] 277.0792 276.8778 281.8665 278.3258 276.3106 275.7861 274.0958 272.3878
## [41] 272.7625 271.2279 272.2817 274.7370 272.1027 279.0301 273.8234 273.2200
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Close, h = 48)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

#mape data testing
mape.koyck <- MAPE(fore.koyck$forecasts, test$Open)

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

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

## $MAPE_testing
## [1] 0.006201109
## 
## $MAPE_training.MAPE
## [1] 0.007829506

Output di atas menunjukkan bahwa nilai MAPE testing lebih kecil dari nilai MAPE training. Selisih antara MAPE testing dan MAPE training juga tidak jauh berbeda yang mengindikasikan bahwa model ini tidak mengalami underfitted atau overfitted.

Regression with Distributed Lag (lag=2)

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

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7737 -0.3257 -0.0721  0.1853  6.1234 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.44396    2.41897   0.597  0.55149    
## x.t          0.03007    0.02192   1.372  0.17213    
## x.1          1.02370    0.02949  34.708  < 2e-16 ***
## x.2         -0.05834    0.02194  -2.659  0.00873 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8713 on 144 degrees of freedom
## Multiple R-squared:  0.9905, Adjusted R-squared:  0.9903 
## F-statistic:  4984 on 3 and 144 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 385.1715 400.1575

AIC(model.dlm)

## [1] 385.1715

BIC(model.dlm)

## [1] 400.1575

Hipotesis:

H0: β_i=0

H1: β_i≠0

Berdasarkan output di atas, hanya nilai p-value X.1 dan X.2 yang mengindikasikan bahwa penduga parameter ini signifikan. Dari pengujian hipotesis, X.1 dan X.2 berhasil menolak H0 sehingga X.1 dan X.2 berpengaruh, sedangkan b0 dan X.t tidak berhasil menolak H0. Menurut model distribusi lag, nilai Open tidak dipengaruhi oleh nilai Close pada saat itu, namun dipengaruhi oleh nilai Close pada satu dan dua periode sebelumnya.

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

## $forecasts
##  [1] 278.8064 277.4056 277.8756 280.6010 279.7269 279.8490 282.6220 282.7445
##  [9] 281.7546 280.0156 278.2963 275.7880 277.8155 278.8335 278.6369 277.6401
## [17] 278.5323 282.4478 280.9557 279.7432 277.9645 277.8377 276.9280 276.9745
## [25] 276.7668 276.6610 276.1490 275.9464 275.0518 275.5289 272.5611 273.8021
## [33] 275.4478 277.6066 277.6076 281.9283 278.8533 276.9454 276.4027 274.8541
## [41] 273.3097 273.5686 272.2502 273.2112 275.3535 273.2498 279.1895 274.7292
## 
## $call
## forecast.dlm(model = model.dlm, x = test$Close, h = 48)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

#mape data testing
mape.dlm <- MAPE(fore.dlm$forecasts, test$Open)

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

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

## $MAPE_testing
## [1] 0.0007290361
## 
## $MAPE_training.MAPE
## [1] 0.001551427

Output di atas menunjukkan bahwa nilai MAPE testing lebih kecil dari nilai MAPE training. Selisih antara MAPE testing dan MAPE training juga tidak jauh berbeda yang mengindikasikan bahwa model ini tidak mengalami underfitted atau overfitted.

Hipotesis:

H0: β_i=0

H1: β_i≠0

Berdasarkan output di atas, hanya nilai p-value X.1 < 0.05 yang mengindikasikan bahwa penduga parameter ini signifikan. Dari pengujian hipotesis, X.1 berhasil menolak H0 sehingga X.1 berpengaruh, sedangkan b0, X.t, X.2, X.3, X.4, dan X.5 tidak berhasil menolak H0. Menurut model distribusi lag (optimum), nilai Open hanya dipengaruhi oleh nilai Close pada satu periode sebelumnya.

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

## $forecasts
##  [1] 278.8844 277.4318 277.8639 280.6224 279.8144 279.8529 282.5849 282.8144
##  [9] 281.7692 279.9624 278.2696 275.7965 277.8153 278.9084 278.7148 277.6247
## [17] 278.5197 282.4676 281.0541 279.7202 277.8899 277.8393 276.9617 277.0053
## [25] 276.7858 276.6880 276.1633 275.9550 275.0689 275.5299 272.5969 273.7875
## [33] 275.4776 277.7001 277.6566 281.9023 278.9026 276.9137 276.3030 274.9044
## [41] 273.3379 273.5645 272.3019 273.2397 275.3733 273.3478 279.1371 274.8116
## 
## $call
## forecast.dlm(model = model.dlm2, x = test$Close, h = 48)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.0007382915
## 
## $MAPE_training.MAPE
## [1] 0.001573662

Output di atas menunjukkan bahwa nilai MAPE testing lebih kecil dari nilai MAPE training. Selisih antara MAPE testing dan MAPE training juga tidak jauh berbeda yang mengindikasikan bahwa model ini tidak mengalami underfitted atau overfitted.

MODEL AUTOREGRESSIVE

#library(dLagM)
model.ardl = ardlDlm(x = train$Close, y = train$Open, 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 = 150
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7910 -0.3085 -0.0893  0.1764  6.1289 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.31021    2.41780   0.542   0.5887    
## X.t          0.02994    0.02191   1.367   0.1739    
## X.1          1.01833    0.02931  34.745   <2e-16 ***
## Y.1         -0.05233    0.02158  -2.425   0.0166 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8719 on 145 degrees of freedom
## Multiple R-squared:  0.9904, Adjusted R-squared:  0.9902 
## F-statistic:  4980 on 3 and 145 DF,  p-value: < 2.2e-16

AIC(model.ardl)

## [1] 387.9335

BIC(model.ardl)

## [1] 402.9533

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

## $forecasts
##  [1] 278.7925 277.4129 277.8714 280.5760 279.7216 279.8471 282.5990 282.7336
##  [9] 281.7591 280.0281 278.3097 275.8059 277.8028 278.8155 278.6307 277.6430
## [17] 278.5243 282.4137 280.9522 279.7528 277.9766 277.8390 276.9297 276.9713
## [25] 276.7624 276.6566 276.1470 275.9430 275.0522 275.5211 272.5734 273.7939
## [33] 275.4252 277.5801 277.5961 281.8932 278.8590 276.9632 276.4065 274.8602
## [41] 273.3178 273.5633 272.2506 273.1998 275.3275 273.2514 279.1465 274.7378
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Close, h = 48)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.0007137695
## 
## $MAPE_training.MAPE
## [1] 0.001525473

Output di atas menunjukkan bahwa nilai MAPE testing lebih kecil dari nilai MAPE training. Selisih antara MAPE testing dan MAPE training juga tidak jauh berbeda yang mengindikasikan bahwa model ini tidak mengalami underfitted atau overfitted.

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

## $p
##   Close
## 1    14
## 
## $q
## [1] 4
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  388.3069 388.2546 388.6754 388.2191 387.4607 387.7174 387.7701 387.0566
## p = 2  388.3073 390.2056 390.6190 390.1672 389.3971 389.6751 389.7392 388.9914
## p = 3  389.7309 389.7309 391.7030 391.3806 390.5694 390.8701 390.9067 390.1089
## p = 4  383.4454 384.6310 384.6310 386.1779 385.8579 386.1069 386.0749 385.5465
## p = 5  387.1548 388.3276 383.7000 383.7000 385.6973 386.0385 385.8957 385.1800
## p = 6  388.9442 390.2352 385.8584 385.0006 385.0006 386.6522 386.5192 385.8296
## p = 7  390.0265 390.9053 386.1040 385.1772 385.8541 385.8541 387.5211 386.6316
## p = 8  388.3050 389.4032 385.5454 384.9660 385.7663 386.9417 386.9417 388.2578
## p = 9  387.7174 388.8970 384.5953 383.7989 384.6838 386.0638 387.4533 387.4533
## p = 10 388.3087 389.2967 385.0032 384.4778 385.3830 386.7566 388.3797 389.5831
## p = 11 388.0492 388.9358 384.8746 384.2331 385.1609 386.5669 388.1477 389.2701
## p = 12 388.1037 389.1995 385.0346 384.2143 385.2084 386.4348 388.1271 389.1036
## p = 13 387.6927 388.7839 384.6458 384.0472 385.2003 386.5174 388.2312 389.1129
## p = 14 386.6808 387.4928 383.5575 383.0334 384.2071 385.7101 387.5685 388.8390
## p = 15 388.2170 389.2788 385.3083 384.7243 385.9282 387.2564 388.9155 390.0624
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  386.0751 386.1026 386.4390 386.0212 386.4144 386.3861 386.5806
## p = 2  387.9696 388.0249 388.3410 387.9254 388.3229 388.2859 388.4520
## p = 3  389.1851 389.1781 389.5092 389.1907 389.5948 389.5900 389.7851
## p = 4  384.9331 384.9331 385.1743 385.0156 385.4956 385.2440 385.5489
## p = 5  384.8176 384.7260 384.8365 384.5285 385.0361 384.8438 385.1392
## p = 6  385.3865 385.3863 385.5406 385.3995 385.8960 385.8387 386.1978
## p = 7  386.2992 386.1990 386.2074 386.0552 386.5798 386.5962 386.9149
## p = 8  387.9187 387.8333 387.7864 387.7606 388.2807 388.3071 388.5874
## p = 9  389.4457 389.3049 389.1455 389.2024 389.7368 389.8024 390.0921
## p = 10 389.5831 391.3011 391.1449 391.1998 391.7340 391.8010 392.0920
## p = 11 391.2694 391.2694 392.6837 392.6275 393.1568 393.3249 393.5692
## p = 12 391.1032 392.7612 392.7612 394.5633 395.0868 395.2692 395.4878
## p = 13 391.1127 392.7964 394.7913 394.7913 396.6854 396.7787 397.0601
## p = 14 390.8260 392.4214 394.4001 395.8267 395.8267 397.7595 398.1621
## p = 15 392.0601 393.6994 395.6774 397.3207 398.2566 398.2566 400.0559
## 
## $min.Stat
## [1] 383.0334

PEMODELAN DLM dan ARDL dengan library dynlm

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

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

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

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

#Ringkasan Model
summary(cons_lm1)

## 
## Time series regression with "ts" data:
## Start = 2, End = 150
## 
## Call:
## dynlm(formula = Open ~ Close + L(Close), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7577 -0.2787 -0.0828  0.1172  6.1023 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.23678    2.41632   0.098    0.922    
## Close        0.02842    0.02226   1.277    0.204    
## L(Close)     0.97122    0.02231  43.534   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8863 on 146 degrees of freedom
## Multiple R-squared:   0.99,  Adjusted R-squared:  0.9899 
## F-statistic:  7225 on 2 and 146 DF,  p-value: < 2.2e-16

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 150
## 
## Call:
## dynlm(formula = Open ~ Close + L(Open), data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.3187  -1.2429  -0.1319   1.1273  19.9969 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.17515    7.35632   0.432    0.667    
## Close        0.54441    0.04914  11.078  < 2e-16 ***
## L(Open)      0.44471    0.04918   9.042 8.69e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.653 on 146 degrees of freedom
## Multiple R-squared:  0.9104, Adjusted R-squared:  0.9091 
## F-statistic: 741.3 on 2 and 146 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 150
## 
## Call:
## dynlm(formula = Open ~ Close + L(Close) + L(Open), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7910 -0.3085 -0.0893  0.1764  6.1289 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.31021    2.41780   0.542   0.5887    
## Close        0.02994    0.02191   1.367   0.1739    
## L(Close)     1.01833    0.02931  34.745   <2e-16 ***
## L(Open)     -0.05233    0.02158  -2.425   0.0166 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8719 on 145 degrees of freedom
## Multiple R-squared:  0.9904, Adjusted R-squared:  0.9902 
## F-statistic:  4980 on 3 and 145 DF,  p-value: < 2.2e-16

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 150
## 
## Call:
## dynlm(formula = Open ~ Close + L(Close) + L(Close, 2), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7737 -0.3257 -0.0721  0.1853  6.1234 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.44396    2.41897   0.597  0.55149    
## Close        0.03007    0.02192   1.372  0.17213    
## L(Close)     1.02370    0.02949  34.708  < 2e-16 ***
## L(Close, 2) -0.05834    0.02194  -2.659  0.00873 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8713 on 144 degrees of freedom
## Multiple R-squared:  0.9905, Adjusted R-squared:  0.9903 
## F-statistic:  4984 on 3 and 144 DF,  p-value: < 2.2e-16

#SSE
deviance(cons_lm1)

## [1] 114.6956

deviance(cons_lm2)

## [1] 1027.93

deviance(cons_lm3)

## [1] 110.2263

deviance(cons_lm4)

## [1] 109.3198

Uji Non Autokorelasi

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

## Encompassing test
## 
## Model 1: Open ~ Close + L(Close)
## Model 2: Open ~ Close + L(Open)
## Model E: Open ~ Close + L(Close) + L(Open)
##           Res.Df Df         F  Pr(>F)    
## M1 vs. ME    145 -1    5.8792 0.01655 *  
## M2 vs. ME    145 -1 1207.2166 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Diagnostik
#durbin watson
dwtest(cons_lm1)

## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 1.9216, p-value = 0.2886
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm2)

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

dwtest(cons_lm3)

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

dwtest(cons_lm4)

## 
##  Durbin-Watson test
## 
## data:  cons_lm4
## DW = 1.948, p-value = 0.347
## alternative hypothesis: true autocorrelation is greater than 0

Uji Heterogenitas

#heterogenitas
bptest(cons_lm1)

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

bptest(cons_lm2)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 32.672, df = 2, p-value = 8.043e-08

bptest(cons_lm3)

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

bptest(cons_lm4)

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

Uji Normalitas

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

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm1)
## W = 0.60252, p-value < 2.2e-16

shapiro.test(residuals(cons_lm2))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.76871, p-value = 4.709e-14

shapiro.test(residuals(cons_lm3))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm3)
## W = 0.6518, p-value < 2.2e-16

shapiro.test(residuals(cons_lm4))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm4)
## W = 0.66408, p-value < 2.2e-16

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          0.0062011089
## DLM 1          0.0007290361
## DLM 2          0.0007382915
## Autoregressive 0.0007137695

PLOT

par(mfrow=c(1,1))
plot(test$Close, test$Open, type="b", col="black", ylim=c(270.7,318.3))
points(test$Close, fore.koyck$forecasts,col="red")
lines(test$Close, fore.koyck$forecasts,col="red")
points(test$Close, fore.dlm$forecasts,col="blue")
lines(test$Close, fore.dlm$forecasts,col="blue")
points(test$Close, fore.dlm2$forecasts,col="orange")
lines(test$Close, fore.dlm2$forecasts,col="orange")
points(test$Close, fore.ardl$forecasts,col="green")
lines(test$Close, 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)