Regresi dengan Peubah Lag
Pendahuluan
Package
library(dLagM)## Warning: package 'dLagM' was built under R version 4.1.3## Loading required package: nardl## Warning: package 'nardl' was built under R version 4.1.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.1.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(dynlm)
library(MLmetrics)## Warning: package 'MLmetrics' was built under R version 4.1.3##
## Attaching package: 'MLmetrics'## The following object is masked from 'package:dLagM':
##
## MAPE## The following object is masked from 'package:base':
##
## Recalllibrary(lmtest)
library(car)## Warning: package 'car' was built under R version 4.1.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.1.3library(readxl)## Warning: package 'readxl' was built under R version 4.1.3Data
Rainfall Timeseries data(https://www.kaggle.com/datasets/poojag718/rainfall-timeseries-data) lahir atas dasar projek Prediction of Worldwide Energy Resources (POWER) menyediakan data meteorologi dari penelitian NASA untuk mendukung energi terbarukan. Data yang dikumpulkan oleh POWER ini adalah pengamatan yang dikumpulkan dalam frekuensi bulanan untuk lintang dan bujur tertentu dalam periode tahun 2000-2020. Data terdiri atas variabel-variabel diantaranya yaitu, Spesific Humidity, Relative Humidity, Temperature, Participation. Dalam tugas Regresi dengan peubah lag kali ini, saya menggunakan peubah Relative Humidity sebagai peubah respon dan Temperature sebagai peubah penjelas. Relative humidity secara umum mampu mewakili pengertian kelembaban dimana rasio antara tekanan uap air aktual pada temperatur tertentu dengan tekanan uap air jenuh pada temperatur tersebut (Fathulrohman dan Saepuloh, 2018). Kelembapan dipengaruhi oleh beberapa faktor, salah satunya yaitu temperature. Semakin meningkat kelembaban udara, semakin meningkat pula temperature dan sebaliknya semakin menurun kelembaban relatif udara maka temperature akan semakin menurun (Sandi et al.2017).
#IMPORT DATA
#membuka file data
data <- read_excel("D:/SEMESTER 5/METODE PERAMALAN DERET WAKTU/Canada_gold.xlsx", sheet = "Sheet3")
data$Y=as.numeric(data$Y)
data$X=as.numeric(data$X)
str(data)## tibble [252 x 3] (S3: tbl_df/tbl/data.frame)
## $ t: num [1:252] 1 2 3 4 5 6 7 8 9 10 ...
## $ Y: num [1:252] 48.2 50.8 42.9 55.7 70.9 ...
## $ X: num [1:252] 23.9 25.8 26.7 22.5 19.1 ...#SPLIT DATA
train<-data[1:200,]
test<-data[201:252,]#data time series
train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(data)Dilakukan splitting data dengan dua kategori data train dan data test. Data train mempunyai proporsi jumlah data sebanyak 80% dari total amatan sedangkan data test 20% dari total amatan.
Eksplorasi Data Analisis
Plot Time Series
Time Series Plot Of Relative Humidity
data.ts<-ts(data$Y)
plot(data.ts, main = "Time Series Plot of Relative Humidity 2000-2020", xlab = "Period", ylab="Humidity")
points(data.ts)
Berdasarkan plot time series diatas, terlihat bahwa peubah relative humidity atau kelembapan relatif mempunyai pola data yang stasioner.Peubah ini akan menjadi peubah respon yang akan dilakukan analisis regresi.
Time Series Plot Of Temperature
data.ts<-ts(data$X)
plot(data.ts, main = "Time Series Plot of Temperature 2000-2020", xlab = "Period", ylab="Temperature")
points(data.ts)
Berdasarkan time series plot diatas, terlihat bahwa peubah temperature memiliki pola data stasioner. Peubah temperature akan menjadi peubah penjelass untuk analisis regresi.
Scatter Plot Peubah Relative Humidity (Y) dan Temperature (X)
cor(data$X,data$Y)## [1] -0.9508894plot(data$X,data$Y,pch=20,col="red",main="Scatter Plot Peubah Relative Humidity (Y) dan Temperature (X) ")
Berdasarkan scatter plot diatas terlihat hubungan antra peubah Relative Humidity (Y) dan Temperature (X) memiliki hubungan linier negatif dengan nilai korelasi -0.9508894.
Model Koyck
Model Koyck adalah model yang mengasumsikan bahwa koefesien dari variabel yang mengalami keterlambatan dapat menurun secara geometris (Gujarti,2006). Asumsi dasar dalam model ini adalah semakin jauh jarak lag pada peubah bebas dari periode sekarang maka semakin kecil pengaruh peubah lag terhadap variabel terbebas(Aqibah et al. 2020).
#MODEL KOYCK
model.koyck = dLagM::koyckDlm(x = train$X, y = train$Y)
summary(model.koyck)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.7651 -2.9602 -0.4924 2.8894 16.0041
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.95652 5.67217 15.330 < 2e-16 ***
## Y.1 0.22091 0.04817 4.586 8.05e-06 ***
## X.t -2.12551 0.15610 -13.616 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.292 on 196 degrees of freedom
## Multiple R-Squared: 0.9382, Adjusted R-squared: 0.9376
## Wald test: 1183 on 2 and 196 DF, p-value: < 2.2e-16
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: 111.6134 -2.125506 0.2209134Berdasarkan output diatas didapatkan intersep model (b0) sebesar 86.95652, nilai X.t sebesar -2.12551 dan nilai Yt sebesar 0.22091. Sehingga dapat disimupulkan bahwa pemodelan koyck yang dibentuk dari data rainfall antara peubah Relative Humidity (Y) dan Temperature (X) adalah Yt= 86.95652 -2.12551Xt + 0.22091Yt.
Berdasarkan output tersebut, semua peubah memiliki nilai p-value <0.05 yang mengindikasikan bahwa penduga parameter signifikan. Dari pengujian hipotesis baik,b0,Xt,Yt berhasil menolak H0 sehingga berpengaruh. Menurut model Koyck, temperature mempengaruhi nilai relative humidity.
AIC(model.koyck)## [1] 1149.5BIC(model.koyck)## [1] 1162.673#Ramalan
(fore.koyck <- forecast(model = model.koyck, x=test$X, h=52))## $forecasts
## [1] 89.95786 73.99036 64.19268 57.81974 48.76005 44.73942 38.62246 46.55960
## [9] 61.29986 74.94990 92.09999 92.67917 85.60165 72.05028 68.97158 66.39975
## [17] 56.79820 48.30058 47.78366 53.96096 62.91367 73.11714 90.39855 95.00269
## [25] 82.01272 68.66432 65.60920 54.73186 49.28944 40.24402 40.22248 44.78756
## [33] 57.69888 70.56230 88.98395 93.52115 93.37571 82.22718 76.10846 67.16869
## [41] 52.95087 45.58020 44.16447 49.76063 62.60215 75.79023 88.65102 95.19052
## [49] 88.38822 79.65878 68.52689 64.32480
##
## $call
## forecast.koyckDlm(model = model.koyck, x = test$X, h = 52)
##
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"Berdasarkan output diatas, terdapat 52 nilai ramalan dari metode Koyck.
#mape data testing
mape.koyck <- MAPE(fore.koyck$forecasts, test$Y)
mape.koyck## [1] 0.07165713#akurasi data training
mape_train <- dLagM::GoF(model.koyck)["MAPE"]
c("MAPE_testing" = mape.koyck, "MAPE_taining" = mape_train)## $MAPE_testing
## [1] 0.07165713
##
## $MAPE_taining.MAPE
## [1] 0.05417177Berdasarkan output daiatas, nilai MAPE Testing lebih besar dibandingkan dengan nilai MAPE Training. Nilai MAPE testing dan MAPE Training memiliki perbedaan yang tidak jauh, artinya bahwa model ini tidak underfitted atau overfitted.
Regression With Distributed Lag
Model Distribusi Lag adalah model regresi yang tidak hanya mencakup nilai sekarang tetapi nilai masa lalu (lag) dari variabel penjelas (Aqibah et.al. 2020). Kelebihan dari model ini adalah mampu membuat teori statis menjadi dinamis karena waktu diperhitungkan dan panjang beda kala nilai lag diketahui.
Lag = 2
#REGRESSION WITH DISTRIBUTED LAG=2 -> estimasi parameter menggunakan least square
model.dlm = dLagM::dlm(x = train$X,y = train$Y, q=2)
summary(model.dlm)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.4695 -2.7333 -0.4025 2.4854 14.8152
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.84932 0.99179 112.775 <2e-16 ***
## x.t -2.11050 0.09694 -21.772 <2e-16 ***
## x.1 -0.32406 0.14654 -2.211 0.0282 *
## x.2 -0.30919 0.09717 -3.182 0.0017 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.245 on 194 degrees of freedom
## Multiple R-squared: 0.9399, Adjusted R-squared: 0.939
## F-statistic: 1012 on 3 and 194 DF, p-value: < 2.2e-16
##
## AIC and BIC values for the model:
## AIC BIC
## 1 1140.39 1156.831Berdasarkan output diatas, semua nilai variable memiliki nilai p-value < 0.05 yang mempunyai arti bahwa penduga parameter signifikan. Dari pengujian hipotesis nilai b0,Xt,X1,X2 berhasil menolak H0 dengan begitu dapat disimpulkan bahwa menurut model distribusi lag, kadar relative humidity dipengaruhi oleh kadar temperature saat itu dan pada periode sebelumnya.
#Ramalan
(fore.dlm <- forecast(model = model.dlm, x=test$X, h=52))#meramalkan 52 periode ke depan## $forecasts
## [1] 91.42607 74.84058 65.52971 58.09767 48.94613 45.16234 38.54957 46.68153
## [9] 60.23222 73.88373 91.39234 91.87244 86.28462 72.81596 69.85568 66.22056
## [17] 56.95016 48.96622 48.03070 53.55761 62.24919 72.48892 89.76806 94.03663
## [25] 82.43431 70.00839 66.26768 54.59679 50.00782 40.27257 40.59740 44.29545
## [33] 57.09851 69.51846 88.24620 92.54216 93.74014 82.59526 77.10340 67.39448
## [41] 53.54931 46.37901 44.16425 49.44856 61.90369 74.81074 87.95422 94.64885
## [49] 88.47098 80.57592 69.13600 64.95096
##
## $call
## forecast.dlm(model = model.dlm, x = test$X, h = 52)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"Berdasarkan output diatas, terdapat 52 nilai ramalan dari metode distribusi lag.
#mape data testing
mape.dlm <- MAPE(fore.dlm$forecasts,test$Y)
#akurasi data training
mape_train <- GoF(model.dlm)["MAPE"]
c("MAPE_testing" = mape.dlm, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.07190444
##
## $MAPE_training.MAPE
## [1] 0.05263041Berdasarkan output diatas, nilai MAPE Testing lebih besar dibandingkan dengan nilai MAPE Training. Nilai MAPE Testing tidak memiliki perbedaan yang cukup jauh dengan MAPE training yang berarti bahwa model ini tidak mengalami underfitteed atau overfitted.
Regression with Distributed Lag Optimum
#Penentuan lag optimum
finiteDLMauto(formula = Y ~ X,
data = data.frame(train),q.min = 1,q.max = 100,
model.type = "dlm",error.type = "AIC", trace = TRUE)## q - k MASE AIC BIC GMRAE MBRAE R.Adj.Sq Ljung-Box
## 99 99 0.00000 -Inf -Inf 0.00000 0.00000 NaN NaN
## 100 100 0.00000 -Inf -Inf 0.00000 0.00000 NaN NaN
## 98 98 0.06291 404.3415 669.4638 0.05562 0.01897 0.92392 0.044030801
## 97 97 0.06425 406.8325 670.3054 0.06193 0.02323 0.96134 0.866125145
## 96 96 0.08103 451.1310 712.9257 0.07615 -0.11947 0.96032 0.515774252
## 95 95 0.12802 555.8383 815.9264 0.11194 -0.62853 0.91822 0.545159373
## 94 94 0.14366 572.4268 830.7804 0.12976 -0.15732 0.92307 0.004039401
## 93 93 0.16070 596.5121 853.1037 0.14416 -0.12138 0.91947 0.414746990
## 92 92 0.15949 604.1935 858.9960 0.14379 -0.06760 0.92691 0.802082039
## 91 91 0.16043 606.5549 859.5416 0.14617 -0.07045 0.93566 0.944241546
## 90 90 0.15981 608.9558 860.1005 0.14726 0.01935 0.94234 0.967561387
## 89 89 0.17127 623.4370 872.7138 0.16930 0.09768 0.94110 0.552512781
## 88 88 0.17900 639.8147 887.1981 0.17348 0.84015 0.93812 0.835503850
## 87 87 0.17790 641.3554 886.8203 0.17096 0.04912 0.94305 0.823357820
## 86 86 0.18121 648.9309 892.4525 0.15298 0.05370 0.94495 0.919430525
## 85 85 0.18286 650.9150 892.4691 0.17467 -0.09668 0.94835 0.843786509
## 84 84 0.19169 662.6083 902.1706 0.16003 0.28864 0.94714 0.595703416
## 83 83 0.19573 667.9474 905.4943 0.17634 -0.07280 0.94826 0.909997804
## 82 82 0.19513 672.3437 907.8519 0.18013 -0.16023 0.94962 0.681149156
## 81 81 0.19606 679.5256 912.9720 0.18271 -0.19566 0.94978 0.458949496
## 80 80 0.19737 687.5201 918.8819 0.16455 -0.38901 0.95004 0.391677517
## 79 79 0.19952 690.6276 919.8824 0.17482 -0.46075 0.95223 0.285570255
## 78 78 0.20121 694.3974 921.5231 0.17757 -0.16770 0.95392 0.295308594
## 77 77 0.20601 702.1846 927.1593 0.18653 1.12945 0.95351 0.521075690
## 76 76 0.20271 703.6714 926.4736 0.18628 -0.75222 0.95516 0.484044428
## 75 75 0.20956 717.1313 937.7397 0.17986 -0.46996 0.95246 0.502939819
## 74 74 0.20947 722.2274 940.6211 0.16471 -4.40891 0.95317 0.437131349
## 73 73 0.22385 735.8994 952.0576 0.17262 -0.73433 0.95122 0.215724760
## 72 72 0.22675 740.3245 954.2268 0.18409 0.06237 0.95203 0.449925370
## 71 71 0.22709 745.9739 957.6000 0.20445 0.18475 0.95200 0.429737906
## 70 70 0.22786 749.4266 958.7566 0.20127 1.13497 0.95264 0.440130546
## 69 69 0.22927 756.9938 964.0080 0.19592 0.86732 0.95184 0.585897975
## 68 68 0.22832 759.4563 964.1352 0.19805 -0.00122 0.95350 0.474346420
## 67 67 0.22931 761.4573 963.7818 0.19213 1.65029 0.95509 0.477129917
## 66 66 0.23128 764.1844 964.1353 0.19435 0.06490 0.95638 0.455427225
## 65 65 0.23694 772.2177 969.7764 0.21821 0.01603 0.95536 0.546926212
## 64 64 0.23494 775.9680 971.1159 0.20149 0.00383 0.95573 0.732571064
## 63 63 0.23848 781.6536 974.3724 0.20131 -0.05881 0.95580 0.640049673
## 62 62 0.24204 787.8169 978.0884 0.23385 0.07251 0.95569 0.777258444
## 61 61 0.24468 794.1908 981.9971 0.24579 0.31132 0.95557 0.948236476
## 60 60 0.25260 798.6188 983.9423 0.27973 0.19810 0.95593 0.903833832
## 59 59 0.25583 805.1414 987.9645 0.25870 0.29126 0.95543 0.808011613
## 58 58 0.25790 808.6857 988.9911 0.26820 0.35954 0.95568 0.608182200
## 57 57 0.25616 810.6163 988.3870 0.26637 0.25063 0.95649 0.596121770
## 56 56 0.25351 812.9417 988.1607 0.26359 1.13121 0.95745 0.558428174
## 55 55 0.26177 820.6739 993.3244 0.25878 -0.46064 0.95693 0.671928726
## 54 54 0.26130 822.7472 992.8128 0.26904 2.00340 0.95791 0.736197214
## 53 53 0.26496 829.4771 996.9413 0.26891 0.57629 0.95733 0.741858658
## 52 52 0.26846 836.8154 1001.6621 0.28304 0.09338 0.95640 0.414491218
## 51 51 0.26779 840.9390 1003.1521 0.25872 0.14330 0.95649 0.371545833
## 50 50 0.27610 850.5082 1010.0719 0.29193 -0.65266 0.95582 0.458810941
## 49 49 0.28214 855.3814 1012.2800 0.29603 -0.09878 0.95639 0.285971003
## 48 48 0.27908 857.3687 1011.5866 0.30028 -0.08972 0.95719 0.280583217
## 47 47 0.28183 868.1032 1019.6251 0.30744 0.44198 0.95537 0.364398652
## 46 46 0.27845 875.2136 1024.0243 0.27713 -11.13125 0.95444 0.176977141
## 45 45 0.28351 881.4553 1027.5397 0.29816 0.04371 0.95384 0.085641217
## 44 44 0.28554 885.1429 1028.4861 0.30098 0.50906 0.95438 0.054303049
## 43 43 0.28535 889.2133 1029.8006 0.30142 -0.61499 0.95481 0.084485486
## 42 42 0.28705 891.7398 1029.5566 0.30535 0.51723 0.95566 0.091500747
## 41 41 0.29329 899.7898 1034.8215 0.31942 0.06295 0.95461 0.069499734
## 40 40 0.29168 903.7454 1035.9779 0.32162 -0.04141 0.95466 0.130050399
## 39 39 0.29124 907.2480 1036.6670 0.32271 0.29586 0.95497 0.111876680
## 38 38 0.28943 909.5675 1036.1590 0.31714 -0.63254 0.95599 0.115898995
## 37 37 0.28948 913.0390 1036.7890 0.31004 0.88968 0.95642 0.127650471
## 36 36 0.29503 917.9208 1038.8156 0.32729 -0.53177 0.95664 0.146471154
## 35 35 0.29569 922.2787 1040.3047 0.31530 -0.08969 0.95677 0.121145709
## 34 34 0.29980 934.5661 1049.7097 0.30637 0.00020 0.95457 0.077990375
## 33 33 0.29746 937.0883 1049.3360 0.30485 0.28218 0.95492 0.061389144
## 32 32 0.29836 940.5244 1049.8632 0.32087 -0.30597 0.95531 0.069580284
## 31 31 0.30271 944.6566 1051.0732 0.34104 -0.07817 0.95565 0.054924923
## 30 30 0.30387 947.5760 1051.0573 0.34577 -0.08019 0.95612 0.055419575
## 29 29 0.30985 952.7160 1053.2493 0.34744 -0.21359 0.95583 0.071100752
## 28 28 0.30634 956.0946 1053.6669 0.34094 0.00605 0.95593 0.061112765
## 27 27 0.30499 958.4654 1053.0642 0.33917 -0.00958 0.95651 0.059689486
## 26 26 0.30482 961.0883 1052.7009 0.33039 -0.01809 0.95724 0.056912104
## 25 25 0.30790 964.4481 1053.0621 0.33446 0.01134 0.95793 0.053729208
## 24 24 0.30561 966.7301 1052.3331 0.34305 -0.16378 0.95844 0.053785584
## 23 23 0.30801 974.1348 1056.7147 0.33447 0.01217 0.95761 0.073955379
## 22 22 0.31302 982.3828 1061.9274 0.31734 -0.49668 0.95647 0.055568155
## 21 21 0.31059 985.5771 1062.0744 0.32892 0.02589 0.95670 0.062427302
## 20 20 0.31024 988.3185 1061.7565 0.32588 0.06442 0.95723 0.054070990
## 19 19 0.31755 996.8900 1067.2569 0.31752 -0.05356 0.95648 0.073264592
## 18 18 0.31826 999.9084 1067.1925 0.32943 0.19693 0.95698 0.082857649
## 17 17 0.31842 1002.6763 1066.8660 0.32667 0.14613 0.95732 0.079226243
## 16 16 0.31747 1006.7889 1067.8727 0.32071 0.11291 0.95720 0.089036481
## 15 15 0.31718 1009.8947 1067.8611 0.32572 0.53726 0.95753 0.080503839
## 14 14 0.31980 1016.3441 1071.1818 0.32290 0.60169 0.95695 0.135928872
## 13 13 0.32053 1020.2282 1071.9260 0.32732 0.27542 0.95737 0.176712625
## 12 12 0.32164 1023.7866 1072.3332 0.33010 0.34946 0.95770 0.146389513
## 11 11 0.32850 1038.3030 1083.6874 0.31636 0.12028 0.95555 0.111604186
## 10 10 0.33187 1045.6245 1087.8358 0.33992 0.10313 0.95465 0.096762378
## 9 9 0.33940 1060.6102 1099.6375 0.34186 0.28623 0.95195 0.069082195
## 8 8 0.33783 1067.3949 1103.2273 0.33225 0.06007 0.95144 0.021074288
## 7 7 0.33966 1076.9633 1109.5902 0.33031 -0.45994 0.95031 0.080311921
## 6 6 0.34173 1081.7896 1111.2003 0.34515 -0.11227 0.95038 0.071510415
## 5 5 0.35214 1096.8834 1123.0674 0.34737 -0.00247 0.94761 0.114927430
## 4 4 0.36921 1122.3704 1145.3172 0.36891 -1.64344 0.94144 0.064068947
## 3 3 0.38104 1133.5746 1153.2738 0.39168 0.12835 0.93931 0.022010554
## 2 2 0.38334 1140.3896 1156.8309 0.40286 -2.07541 0.93899 0.022635493
## 1 1 0.40564 1156.2136 1169.3868 0.44181 -0.08924 0.93547 0.010367535finiteDLMauto## function (formula, data, x, y, q.min = 1, q.max = 10, k.order = NULL,
## model.type = c("dlm", "poly"), error.type = c("MASE", "AIC",
## "BIC", "GMRAE", "MBRAE", "radj"), trace = FALSE)
## UseMethod("finiteDLMauto")
## <bytecode: 0x0000000024390600>
## <environment: namespace:dLagM>Berdasarkan output diatas, 98 merupakan lag optimum bagi model
#model dlm dengan lag optimum
model.dlm2 = dLagM::dlm(x = train$X,y = train$Y, q= 98)
summary(model.dlm2)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## 1 2 3 4 5 6 7
## 0.1636657 0.6632622 -0.6562336 1.5049513 1.0433059 0.3335745 0.2919524
## 8 9 10 11 12 13 14
## 0.3091465 -0.0143676 0.1649474 -0.8235862 -1.5673230 -0.7924073 -0.8505658
## 15 16 17 18 19 20 21
## -1.0809291 -0.9329250 -0.3792062 0.0129346 0.2973160 0.6197994 0.2732358
## 22 23 24 25 26 27 28
## 0.6297693 0.8989078 0.8490823 0.0597891 0.3032690 0.1905353 0.6505340
## 29 30 31 32 33 34 35
## -0.1198680 -1.2126539 -0.1223671 -0.3088657 -0.6820348 -0.9322031 0.1373285
## 36 37 38 39 40 41 42
## -0.6685456 0.4675113 0.1732994 -1.0473314 -0.6554849 -0.3220342 0.3605649
## 43 44 45 46 47 48 49
## 0.8243481 -0.2628220 0.0104624 0.0984804 0.6288886 -0.2741534 0.1738974
## 50 51 52 53 54 55 56
## 0.2335006 -0.2019355 0.3943774 -0.6481492 0.8160190 -0.1906406 -0.5297615
## 57 58 59 60 61 62 63
## 0.0598337 -1.0683141 -0.0007629 -0.9312764 -0.3766435 -0.1711415 -0.5869109
## 64 65 66 67 68 69 70
## 0.1771825 0.1738813 0.5204170 -0.3043451 0.7561255 -0.2424952 0.4424587
## 71 72 73 74 75 76 77
## 0.8877144 0.2473953 0.7950066 0.4225305 0.4234374 -0.4745669 -0.2689798
## 78 79 80 81 82 83 84
## -0.0411868 -1.0637456 0.0567544 -0.4915377 -0.6099862 0.4000099 -0.6617855
## 85 86 87 88 89 90 91
## 0.0495255 0.0871999 -0.1401397 1.6607363 0.4870670 0.6566350 0.0368729
## 92 93 94 95 96 97 98
## 0.5181272 0.2956907 0.3256770 0.1711096 0.9610794 -0.0043906 0.3919533
## 99 100 101 102
## 1.1169757 -1.1567424 0.3671212 -2.1958288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 270.925909 548.985441 0.494 0.671
## x.t -0.946178 1.301625 -0.727 0.543
## x.1 -0.875632 0.713953 -1.226 0.345
## x.2 -0.635266 0.690505 -0.920 0.455
## x.3 0.215493 0.758290 0.284 0.803
## x.4 -0.258113 0.627270 -0.411 0.721
## x.5 0.087020 0.507313 0.172 0.880
## x.6 0.259048 0.507609 0.510 0.661
## x.7 0.366376 0.542021 0.676 0.569
## x.8 -0.268295 0.578605 -0.464 0.688
## x.9 -0.056728 0.745598 -0.076 0.946
## x.10 0.481963 0.652149 0.739 0.537
## x.11 0.380638 0.662654 0.574 0.624
## x.12 -0.686711 1.169746 -0.587 0.617
## x.13 -0.364030 1.111769 -0.327 0.774
## x.14 -0.024594 1.110903 -0.022 0.984
## x.15 0.153678 1.083895 0.142 0.900
## x.16 -0.679420 1.225409 -0.554 0.635
## x.17 -0.107173 1.194382 -0.090 0.937
## x.18 -0.349935 0.976629 -0.358 0.754
## x.19 -1.208339 1.241638 -0.973 0.433
## x.20 0.450466 1.122757 0.401 0.727
## x.21 -0.547345 0.772400 -0.709 0.552
## x.22 0.062141 0.877081 0.071 0.950
## x.23 0.359851 0.749889 0.480 0.679
## x.24 -0.220678 0.638109 -0.346 0.762
## x.25 0.006736 0.697513 0.010 0.993
## x.26 -0.246953 0.726352 -0.340 0.766
## x.27 -0.479222 0.687148 -0.697 0.558
## x.28 -0.423826 0.821357 -0.516 0.657
## x.29 0.003736 1.208283 0.003 0.998
## x.30 -0.104945 1.108836 -0.095 0.933
## x.31 -0.969312 1.133952 -0.855 0.483
## x.32 0.009542 1.514016 0.006 0.996
## x.33 -0.116857 1.294566 -0.090 0.936
## x.34 -0.328441 1.265986 -0.259 0.820
## x.35 -0.781334 1.433270 -0.545 0.640
## x.36 0.034004 1.558577 0.022 0.985
## x.37 -0.258298 1.406577 -0.184 0.871
## x.38 -0.054554 1.487792 -0.037 0.974
## x.39 0.419456 1.448619 0.290 0.799
## x.40 -0.550340 1.207132 -0.456 0.693
## x.41 -0.347577 1.039171 -0.334 0.770
## x.42 -0.547124 0.976220 -0.560 0.632
## x.43 -0.753020 1.100996 -0.684 0.565
## x.44 0.665455 0.877976 0.758 0.528
## x.45 -0.579823 0.755216 -0.768 0.523
## x.46 -0.614261 0.837428 -0.734 0.540
## x.47 -0.578395 0.638947 -0.905 0.461
## x.48 -0.626972 0.677545 -0.925 0.452
## x.49 -0.313206 0.678813 -0.461 0.690
## x.50 0.111026 0.718143 0.155 0.891
## x.51 0.086298 0.710935 0.121 0.914
## x.52 -0.195950 0.762582 -0.257 0.821
## x.53 0.351154 0.788562 0.445 0.700
## x.54 -0.468520 0.687967 -0.681 0.566
## x.55 -1.132611 1.080731 -1.048 0.405
## x.56 0.345392 1.089251 0.317 0.781
## x.57 0.325984 0.938912 0.347 0.762
## x.58 -0.553244 1.035166 -0.534 0.646
## x.59 -0.817355 1.491165 -0.548 0.639
## x.60 0.435653 1.081626 0.403 0.726
## x.61 -0.213082 0.808374 -0.264 0.817
## x.62 -0.837259 0.970671 -0.863 0.479
## x.63 0.033780 0.898856 0.038 0.973
## x.64 -0.172696 0.742690 -0.233 0.838
## x.65 -0.166053 0.857687 -0.194 0.864
## x.66 0.496125 0.775708 0.640 0.588
## x.67 0.241845 0.623435 0.388 0.735
## x.68 -0.110256 0.623215 -0.177 0.876
## x.69 0.167774 0.564983 0.297 0.795
## x.70 -0.394580 0.606191 -0.651 0.582
## x.71 0.234638 0.829998 0.283 0.804
## x.72 -0.108386 0.844275 -0.128 0.910
## x.73 -0.619559 0.633001 -0.979 0.431
## x.74 -0.187393 0.623721 -0.300 0.792
## x.75 0.877836 0.776354 1.131 0.376
## x.76 -0.647862 0.571941 -1.133 0.375
## x.77 0.031502 0.552818 0.057 0.960
## x.78 0.647138 0.541666 1.195 0.355
## x.79 0.009380 0.691814 0.014 0.990
## x.80 -0.149421 0.840479 -0.178 0.875
## x.81 0.101372 0.696459 0.146 0.898
## x.82 -0.555780 0.899478 -0.618 0.600
## x.83 0.352163 1.058226 0.333 0.771
## x.84 -0.344799 0.781702 -0.441 0.702
## x.85 -0.176775 0.680119 -0.260 0.819
## x.86 0.122881 0.704683 0.174 0.878
## x.87 0.008966 0.680734 0.013 0.991
## x.88 -0.264745 0.822862 -0.322 0.778
## x.89 0.164676 0.753854 0.218 0.847
## x.90 0.886672 1.282756 0.691 0.561
## x.91 -0.211366 1.494580 -0.141 0.900
## x.92 -0.075226 0.918824 -0.082 0.942
## x.93 0.196843 0.984065 0.200 0.860
## x.94 0.145025 1.433562 0.101 0.929
## x.95 -0.183836 0.917299 -0.200 0.860
## x.96 1.036672 0.676876 1.532 0.265
## x.97 0.578385 0.739739 0.782 0.516
## x.98 -0.114442 0.705273 -0.162 0.886
##
## Residual standard error: 4.659 on 2 degrees of freedom
## Multiple R-squared: 0.9985, Adjusted R-squared: 0.9239
## F-statistic: 13.39 on 99 and 2 DF, p-value: 0.07191
##
## AIC and BIC values for the model:
## AIC BIC
## 1 404.3415 669.4638Berdasarkan output diatas, didapatkan nilai intersep model b0 sebesar 270.925909, nilai Xt sebesar -0.946178 dst. Sehingga dapat disimpulkan bahwa pemodelan distribusi lag yang dibentuk dari data rainfall dengan peubah relative humidity dan temperature adalah sebagai berikut: Yt = (270.925909 -0.946178Xt,1 -0.875632Xt-1,1 …-0.114442xt-98,1 )
#ramalan
(fore.dlm2 <- forecast(model = model.dlm2, x=test$X, h=52))## $forecasts
## [1] 81.83019 66.61952 81.50079 55.82707 37.94272 41.84302 41.84890
## [8] 36.69681 45.98685 76.54465 78.24101 88.64021 101.24686 69.34788
## [15] 68.29671 55.86380 48.81404 46.37438 37.17842 34.71320 53.64905
## [22] 74.02093 65.62039 83.73571 89.09567 49.88711 61.30218 55.06551
## [29] 41.09800 30.29716 31.33674 44.02340 46.99753 56.77869 72.55642
## [36] 101.27135 90.15921 60.35777 71.92563 70.20807 46.35260 45.85970
## [43] 51.92533 38.58187 38.61857 75.97877 72.58423 73.63384 87.72526
## [50] 73.74997 64.72626 59.26891
##
## $call
## forecast.dlm(model = model.dlm2, x = test$X, h = 52)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"Berdasarkan output daiatas, terdapat 52 nilai ramalan dari metode Regression with Distributed Lag Optimum.
#akurasi testing
mape.dlm2 <- MAPE(fore.dlm2$forecast, test$Y)
#akurasi data training
mape_train <- GoF(model.dlm2)["MAPE"]
c("MAPE_testing" = mape.dlm2,"MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.1527344
##
## $MAPE_training.MAPE
## [1] 0.00796936Berdasarkan output diatasm nilai MAPE Testing jauh lebih besar dibandingkan dengan nilai MAPE training, oleh karena itu patut diduga bahwa model ini mengalamai underfitted atau overfitted.
Model Autoregressive/Dynamic Regression
Model Autoregressive adalah model regresi time series yang menghubungkan nilai pengamatan aktual dengan nilai pengamatan sebelumnya (Gujarati, 2006). Asumsi dasar dalam model ini adalah meregresikan pengamatan aktual dengan nilai pengamatan sebelumnya untuk peramalan nilai ke depan.
#MODEL AUTOREGRESSIVE
model.ardl = ardlDlm(x = train$X, y = train$Y, p =1, q =1)
summary(model.ardl)##
## Time series regression with "ts" data:
## Start = 2, End = 200
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.3078 -2.7272 -0.4197 2.6712 14.4046
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.84095 6.34303 13.848 < 2e-16 ***
## X.t -1.94158 0.07727 -25.126 < 2e-16 ***
## X.1 -0.18178 0.17054 -1.066 0.287769
## Y.1 0.20742 0.05855 3.543 0.000496 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.242 on 195 degrees of freedom
## Multiple R-squared: 0.94, Adjusted R-squared: 0.9391
## F-statistic: 1018 on 3 and 195 DF, p-value: < 2.2e-16AIC(model.ardl)## [1] 1145.801BIC(model.ardl)## [1] 1162.268#ramalan
(fore.ardl <- forecast(model = model.ardl, x=test$X,h=52))## $forecasts
## [1] 90.03334 75.06064 64.87664 58.38366 49.68724 45.38448 39.54301 46.36887
## [9] 60.44216 73.96631 90.57196 92.29343 85.79405 72.87917 69.09741 66.57770
## [17] 57.64155 49.19067 48.13502 53.77947 62.41978 72.37487 88.87026 94.29537
## [25] 82.69292 69.51636 65.78996 55.68835 49.94402 41.32813 40.67588 44.88405
## [33] 57.02060 69.70083 87.41899 92.85377 92.97260 82.74463 76.33286 67.75881
## [41] 54.14788 46.41384 44.64203 49.69029 61.84380 74.81275 87.47469 94.33017
## [49] 88.53512 80.02641 69.23640 64.61908
##
## $call
## forecast.ardlDlm(model = model.ardl, x = test$X, h = 52)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"Berdasarkan output diatas, terdapat 52 nilai ramalan dari model autoregresive.
#Akurasi testing
mape.ardl <- MAPE(fore.ardl$forecasts, test$Y)
#Akurasi data training
mape_train <- GoF(model.ardl)["MAPE"]
c("MAPE_testing" = mape.ardl, "MAPE_training" = mape_train) ## $MAPE_testing
## [1] 0.07211403
##
## $MAPE_training.MAPE
## [1] 0.05330417Berdasarkan output daiatas, nilai MAPE Testing lebih besar dibandingkan dengan nilai MAPE Training. Nilai MAPE testing dan MAPE Training memiliki perbedaan yang tidak jauh, artinya bahwa model ini tidak underfitted atau overfitted.
#penentuan lag optimum
ardlBoundOrders(data = data.frame(data), formula = Y ~ X)## $p
## X
## 1 15
##
## $q
## [1] 3
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7 q = 8
## p = 1 1452.279 1437.827 1419.795 1407.055 1401.118 1395.820 1392.080 1377.102
## p = 2 1445.023 1439.827 1421.574 1408.829 1403.085 1397.818 1394.056 1379.084
## p = 3 1432.718 1432.718 1422.325 1408.841 1403.629 1399.050 1395.260 1380.400
## p = 4 1407.269 1408.976 1408.976 1410.293 1405.297 1400.344 1396.583 1382.082
## p = 5 1404.612 1405.698 1404.275 1404.275 1405.976 1401.287 1397.490 1383.960
## p = 6 1400.834 1401.973 1400.504 1402.424 1402.424 1403.255 1399.432 1385.850
## p = 7 1392.183 1393.408 1392.345 1394.188 1395.774 1395.774 1394.827 1379.424
## p = 8 1383.811 1385.357 1384.999 1386.892 1388.750 1388.984 1388.984 1381.356
## p = 9 1368.894 1370.274 1369.932 1371.328 1373.222 1373.131 1374.904 1374.904
## p = 10 1364.227 1365.894 1364.851 1366.301 1367.791 1368.403 1370.400 1367.795
## p = 11 1352.082 1353.699 1351.252 1352.630 1353.999 1355.285 1357.217 1353.192
## p = 12 1347.558 1349.094 1346.762 1348.312 1349.545 1350.947 1352.805 1348.921
## p = 13 1341.704 1343.409 1341.704 1343.410 1344.958 1346.145 1348.094 1344.663
## p = 14 1334.147 1336.092 1334.482 1336.049 1337.496 1338.661 1340.615 1337.790
## p = 15 1330.108 1332.003 1329.027 1330.635 1331.994 1333.113 1335.111 1331.943
## q = 9 q = 10 q = 11 q = 12 q = 13 q = 14 q = 15
## p = 1 1371.651 1367.293 1354.665 1345.949 1340.861 1336.112 1329.035
## p = 2 1373.643 1369.255 1356.617 1347.667 1342.348 1337.886 1330.550
## p = 3 1374.863 1370.490 1357.394 1348.415 1343.009 1338.570 1329.994
## p = 4 1376.601 1372.324 1359.248 1350.153 1344.773 1340.377 1331.981
## p = 5 1378.548 1374.064 1361.163 1352.128 1346.768 1342.366 1333.980
## p = 6 1380.069 1375.629 1363.163 1354.122 1348.768 1344.340 1335.954
## p = 7 1375.079 1369.918 1356.574 1348.349 1343.311 1338.688 1330.426
## p = 8 1376.887 1371.548 1357.472 1349.007 1343.334 1338.363 1330.259
## p = 9 1376.808 1371.688 1358.405 1350.496 1345.005 1339.593 1331.449
## p = 10 1367.795 1368.218 1355.521 1348.396 1343.962 1338.000 1330.526
## p = 11 1354.255 1354.255 1355.435 1349.157 1344.185 1339.012 1331.405
## p = 12 1350.379 1351.238 1351.238 1351.065 1346.128 1341.012 1333.271
## p = 13 1345.882 1347.463 1346.484 1346.484 1348.068 1342.905 1335.271
## p = 14 1339.499 1340.915 1340.573 1340.825 1340.825 1342.711 1335.249
## p = 15 1333.297 1335.113 1334.394 1334.640 1336.639 1336.639 1336.784
##
## $min.Stat
## [1] 1329.027#PEMODELAN DLM dan ARDL dengan library dynlm
#sama dengan model dlm p=1
cons_lm1 <- dynlm(Y ~ X+L(X),data = train.ts)
#sama dengan model ardl p=0 q=3
cons_lm2 <- dynlm(Y ~ X+L(Y),data = train.ts)
#sama dengan ardl p=1 q=3
cons_lm3 <- dynlm(Y ~ X+L(X)+L(Y),data = train.ts)
#sama dengan dlm p=2
cons_lm4 <- dynlm(Y ~ X+L(X)+L(X,2),data = train.ts)#Ringkasan Model
summary(cons_lm1)##
## Time series regression with "ts" data:
## Start = 2, End = 200
##
## Call:
## dynlm(formula = Y ~ X + L(X), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.4820 -2.8133 -0.3813 2.5190 13.4948
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 110.11187 0.86952 126.636 <2e-16 ***
## X -1.91508 0.07914 -24.197 <2e-16 ***
## L(X) -0.72059 0.07939 -9.077 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.365 on 196 degrees of freedom
## Multiple R-squared: 0.9361, Adjusted R-squared: 0.9355
## F-statistic: 1436 on 2 and 196 DF, p-value: < 2.2e-16summary(cons_lm2)##
## Time series regression with "ts" data:
## Start = 2, End = 200
##
## Call:
## dynlm(formula = Y ~ X + L(Y), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1889 -2.6591 -0.5487 2.7397 15.0636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 81.74053 2.73606 29.875 <2e-16 ***
## X -1.97847 0.06912 -28.626 <2e-16 ***
## L(Y) 0.26308 0.02650 9.929 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.243 on 196 degrees of freedom
## Multiple R-squared: 0.9396, Adjusted R-squared: 0.939
## F-statistic: 1525 on 2 and 196 DF, p-value: < 2.2e-16summary(cons_lm3)##
## Time series regression with "ts" data:
## Start = 2, End = 200
##
## Call:
## dynlm(formula = Y ~ X + L(X) + L(Y), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.3078 -2.7272 -0.4197 2.6712 14.4046
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.84095 6.34303 13.848 < 2e-16 ***
## X -1.94158 0.07727 -25.126 < 2e-16 ***
## L(X) -0.18178 0.17054 -1.066 0.287769
## L(Y) 0.20742 0.05855 3.543 0.000496 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.242 on 195 degrees of freedom
## Multiple R-squared: 0.94, Adjusted R-squared: 0.9391
## F-statistic: 1018 on 3 and 195 DF, p-value: < 2.2e-16summary(cons_lm4)##
## Time series regression with "ts" data:
## Start = 3, End = 200
##
## Call:
## dynlm(formula = Y ~ X + L(X) + L(X, 2), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.4695 -2.7333 -0.4025 2.4854 14.8152
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.84932 0.99179 112.775 <2e-16 ***
## X -2.11050 0.09694 -21.772 <2e-16 ***
## L(X) -0.32406 0.14654 -2.211 0.0282 *
## L(X, 2) -0.30919 0.09717 -3.182 0.0017 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.245 on 194 degrees of freedom
## Multiple R-squared: 0.9399, Adjusted R-squared: 0.939
## F-statistic: 1012 on 3 and 194 DF, p-value: < 2.2e-16#SSE
deviance(cons_lm1)## [1] 3734.411deviance(cons_lm2)## [1] 3529.038deviance(cons_lm3)## [1] 3508.594deviance(cons_lm4)## [1] 3496.187Uji Diagnostik Model
Uji Non Autokorelasi
#Uji Model
if(require("lmtest"))encomptest(cons_lm1, cons_lm2)## Encompassing test
##
## Model 1: Y ~ X + L(X)
## Model 2: Y ~ X + L(Y)
## Model E: Y ~ X + L(X) + L(Y)
## Res.Df Df F Pr(>F)
## M1 vs. ME 195 -1 12.5504 0.0004956 ***
## M2 vs. ME 195 -1 1.1362 0.2877693
## ---
## 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.6193, p-value = 0.002647
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm2)##
## Durbin-Watson test
##
## data: cons_lm2
## DW = 2.0435, p-value = 0.5793
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm3)##
## Durbin-Watson test
##
## data: cons_lm3
## DW = 1.9418, p-value = 0.3021
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm4)##
## Durbin-Watson test
##
## data: cons_lm4
## DW = 1.6734, p-value = 0.008721
## alternative hypothesis: true autocorrelation is greater than 0Uji Heterogenitas
#heterogenitas
bptest(cons_lm1)##
## studentized Breusch-Pagan test
##
## data: cons_lm1
## BP = 4.2101, df = 2, p-value = 0.1218bptest(cons_lm2)##
## studentized Breusch-Pagan test
##
## data: cons_lm2
## BP = 6.7886, df = 2, p-value = 0.03356bptest(cons_lm3)##
## studentized Breusch-Pagan test
##
## data: cons_lm3
## BP = 7.7844, df = 3, p-value = 0.05068bptest(cons_lm4)##
## studentized Breusch-Pagan test
##
## data: cons_lm4
## BP = 8.0893, df = 3, p-value = 0.0442Uji Normalitas
#shapiro wilk
shapiro.test(residuals(cons_lm1))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm1)
## W = 0.98933, p-value = 0.145shapiro.test(residuals(cons_lm2))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm2)
## W = 0.98127, p-value = 0.009309shapiro.test(residuals(cons_lm3))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm3)
## W = 0.98133, p-value = 0.009466shapiro.test(residuals(cons_lm4))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm4)
## W = 0.97646, p-value = 0.00206#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.07165713
## DLM 1 0.07190444
## DLM 2 0.15273441
## Autoregressive 0.07211403plot(ts(test$X), ts(test$Y), type="b", col="black")
points(test$X, fore.koyck$forecasts,col="red")
#Plot
par(mfrow=c(1,1))
plot(test$X, test$Y, type="b", col="black")
points(test$X, fore.koyck$forecasts,col="red")
points(test$X, fore.dlm$forecasts,col="blue")
points(test$X, fore.dlm2$forecasts,col="orange")
points(test$X, 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)