Tugas STA1341 MPDW 2 - Regresi dengan Peubah Lag
Pendahuluan
Tujuan
Melakukan regresi data deret waktu untuk menganalisis hubungan antara jumlah pengunjung halaman website dengan jumlah halaman harian yang dimuat. Serta membandingkan keakuratan peramalan antara model koyck, distibute lag, dan model autoregressive.
Data yang digunakan
Jenis data yang digunakan merupakan data time series yaitu pengunjung harian situs web dari tahun 2018 - 2019. Data tersebut merupakan data sekunder yang bersumber dari https://www.kaggle.com/datasets/bobnau/daily-website-visitors. Peubah yang digunakan untuk analisis adalah Page Loads, dan Unique visits dengan jumlah amatan 729.
#install.packages("dLagM")
#install.packages("dynlm")
#install.packages("MLmetrics")
library(dLagM) #bisa otomatis timeseries datanya
library(dynlm) #data harus timeseries
library(MLmetrics) #MAPE
library(lmtest)
library(car)
library(readxl)#membuka file data
data <- read_excel("E:/SEM 5/MPDW/pertemuan 4/daily-website-visitors.xlsx")
str(data)## tibble [729 x 3] (S3: tbl_df/tbl/data.frame)
## $ Row : num [1:729] 1 2 3 4 5 6 7 8 9 10 ...
## $ Page.Loads : num [1:729] 1709 3389 3391 3698 3511 ...
## $ Unique.Visits: num [1:729] 1120 2208 2282 2310 2272 ...data## # A tibble: 729 x 3
## Row Page.Loads Unique.Visits
## <dbl> <dbl> <dbl>
## 1 1 1709 1120
## 2 2 3389 2208
## 3 3 3391 2282
## 4 4 3698 2310
## 5 5 3511 2272
## 6 6 2279 1498
## 7 7 2822 1869
## 8 8 4240 2784
## 9 9 4319 2844
## 10 10 4224 2866
## # ... with 719 more rowsSelanjutnya, mengubah nama kolom unique visits dan page loads menjadi peubah Xt dan Yt:
t <- data$Row
Xt <- data$Unique.Visits
Yt <- data$Page.Loads
data <- data.frame(t, Xt, Yt)Membagi data 80% dari jumlah seluruh amatan (584 amatan) untuk data train dan 20% lainnya (146 amatan) untuk data test.
#SPLIT DATA
train<-data[1:584,]
test<-data[585:729,]#data time series
train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(data)Pembahasan
1. Explorasi Data
Plot peubah x
#Time Series Plot
data.ts<-ts(Xt)
plot(data.ts, main = "Time Series Plot of Unique Visits", xlab = "Period", ylab="Visitors")
points(data.ts)
Dari plot di atas, terlihat bahwa peubah Jumlah pengunjung halaman web memiliki pola data yang stasioner. Peubah ini akan menjadi peubah penjelas untuk analisis selanjutnya.
Plot peubah y
data.ts1<-ts(Yt)
plot(data.ts1, main = "Time Series Plot Daily Number of Page Loads", xlab = "Period", ylab="Page Loads")
points(data.ts1)
Dari plot di atas terlihat bahwa peubah Page Loads memiliki pola data yang stasioner. Peubah page loads akan dijadikan peubah respon untuk analisis selanjutnya.
Korelasi antar peubah
#Korelasi x dan y
cor(Xt, Yt)## [1] 0.9938377plot(Xt, Yt, pch = 20, col = "darkblue", main = "Scatter Plot Page Loads & Visitors")
Dari scatter plot di atas terlihat bahwa terdapat hubungan yang kuat dan positif antara peubah jumlah pengunjung halaman web dan jumlah halaman yang dibuka, dengan nilai korelasinya sebesar 0.993.
2. Model KOYCK
model.koyck <- dLagM :: koyckDlm(x=train$Xt , y=train$Yt)
summary(model.koyck)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -508.607 -163.337 -7.981 144.809 774.455
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 288.06445 64.77465 4.447 1.04e-05 ***
## Y.1 0.14120 0.03856 3.662 0.000273 ***
## X.t 1.16179 0.07414 15.670 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 227.2 on 580 degrees of freedom
## Multiple R-Squared: 0.9745, Adjusted R-squared: 0.9744
## Wald test: 6318 on 2 and 580 DF, p-value: < 2.2e-16
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: 335.4254 1.161789 0.1411966Berdasarkan ouput, terlihat semua parameter signifikan. Persamaan yang diperoleh dari model Koyck adalah \(Y_{t}= 288.064 + 1.1617X_{t} + 0.1412Y_{t-1}\). Menurut model Koyck, nilai y (jumlah halaman yang dimuat) dipengaruhi jumlah pengunjung halaman website dan lag dari jumlah halaman yang dimuat.
AIC(model.koyck)## [1] 7985.996BIC(model.koyck)## [1] 8003.468#ramalan
(fore.koyck <- forecast(model = model.koyck, x=test$Xt, h=145))## $forecasts
## [1] 3776.820 3183.255 2207.192 2344.719 3609.575 3896.214 3762.418 3665.687
## [9] 3154.783 2215.951 2600.387 3962.842 4100.612 3663.482 3081.279 3212.844
## [17] 2344.975 2699.930 3822.380 4400.271 4539.957 4393.544 3350.498 2371.382
## [25] 2534.038 3285.446 4145.543 4763.069 4536.579 3813.335 2645.855 3122.318
## [33] 4547.724 4707.162 4862.118 4873.541 3889.958 2830.942 3437.737 5233.568
## [41] 5459.250 5387.716 5115.052 4288.860 3210.243 3808.462 5245.250 5473.679
## [49] 5519.874 5386.982 4374.889 3165.463 4098.395 5228.098 5248.194 5307.959
## [57] 5050.348 4316.901 3358.264 4073.337 5427.873 5917.709 5696.425 5509.500
## [65] 4722.136 3407.350 4269.640 5563.637 5856.715 5722.666 5632.870 4681.466
## [73] 3531.727 4366.203 5719.009 5896.079 6285.883 5816.955 4549.455 3320.231
## [81] 4072.614 5649.673 5895.584 5836.201 5397.954 4547.221 3451.198 4333.920
## [89] 5760.923 5953.116 6183.566 6376.432 5088.519 3721.646 4386.048 5779.901
## [97] 6435.615 6163.397 6005.297 5146.486 3926.173 4777.405 6421.862 6738.865
## [105] 6794.080 6431.266 4930.126 3858.446 4599.382 6127.191 5933.963 5567.438
## [113] 4818.612 4392.227 3746.481 4660.251 6376.276 6968.271 7033.270 6707.853
## [121] 5645.340 4291.704 5118.302 7110.141 6858.121 6951.495 6535.980 5012.295
## [129] 3598.190 3900.416 5082.804 5003.453 4932.998 4261.993 3104.212 2093.794
## [137] 2222.985 2621.131 2048.820 1770.508 2112.278 2118.710 2003.439 2080.107
## [145] 2628.840
##
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 145)
##
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"#mape data testing
mape.koyck <- MAPE(fore.koyck$forecasts, test$Yt)
#akurasi data training
mape_train <- dLagM::GoF(model.koyck)["MAPE"]
c("MAPE_testing" = mape.koyck, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.07187
##
## $MAPE_training.MAPE
## [1] 0.04773302Output tersebut menunjukkan bahwa nilai MAPE Testing lebih besar daripada nilai MAPE Training dengan selisih yang tidak terlalu jauh yaitu 0.02413698. Hal tersebut mengindikasikan bahwa model ini tidak underfitted atau overfitted.
3. Regression with Distributed Lag
Regression with Distributed Lag (lag optimum)
#penentuan lag optimum
finiteDLMauto(formula = Yt ~ Xt,
data = data.frame(train), q.min = 1, q.max = 4 ,
model.type = "dlm", error.type = "AIC", trace = TRUE)## q - k MASE AIC BIC GMRAE MBRAE R.Adj.Sq Ljung-Box
## 4 4 0.11851 7262.535 7293.077 0.12097 0.19016 0.99216 3.069477e-09
## 3 3 0.11928 7279.633 7305.822 0.12011 -0.33040 0.99207 1.012960e-09
## 2 2 0.12012 7292.561 7314.393 0.12429 0.03619 0.99205 8.507790e-10
## 1 1 0.12372 7317.970 7335.443 0.13345 0.03470 0.99185 7.112289e-11#model dlm dengan lag optimum
model.dlm2 = dLagM::dlm(x = train$Xt,y = train$Yt , q = 4) #terdapat lag yang tidak signifikan sehingga dapat dikurangi jumlah lagnya
summary(model.dlm2)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -451.84 -88.78 -6.12 80.78 744.27
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 134.181707 23.195852 5.785 1.19e-08 ***
## x.t 1.441256 0.010000 144.119 < 2e-16 ***
## x.1 0.018630 0.016573 1.124 0.2614
## x.2 -0.031031 0.018596 -1.669 0.0957 .
## x.3 0.006693 0.016561 0.404 0.6862
## x.4 -0.017300 0.009976 -1.734 0.0834 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 125.8 on 574 degrees of freedom
## Multiple R-squared: 0.9922, Adjusted R-squared: 0.9922
## F-statistic: 1.465e+04 on 5 and 574 DF, p-value: < 2.2e-16
##
## AIC and BIC values for the model:
## AIC BIC
## 1 7262.535 7293.077Berdasarkan output, menurut model distributed lag, nilai y (jumlah halaman yang dimuat) pada taraf nyata 10% dipengaruhi jumlah pengunjung halaman website pada waktu sekarang, jumlah pengunjung halaman website 2 periode sebelumnya, dan jumlah pengunjung halaman website 4 periode sebelumnya.
AIC(model.dlm2)## [1] 7262.535BIC(model.dlm2)## [1] 7293.077#ramalan
(fore.dlm2 <- forecast(model = model.dlm2, x=test$Xt, h=145))## $forecasts
## [1] 3765.266 3001.686 1884.919 2230.241 3802.144 3952.654 3719.806 3620.636
## [9] 2989.262 1905.578 2548.466 4201.054 4142.895 3556.550 2904.842 3154.088
## [17] 2069.922 2646.154 4016.748 4533.194 4599.291 4381.989 3095.916 2046.109
## [25] 2431.540 3358.835 4308.285 4930.457 4528.850 3646.519 2313.811 3101.570
## [33] 4825.825 4782.746 4929.662 4918.206 3671.852 2516.244 3464.028 5616.741
## [41] 5598.743 5442.405 5114.359 4110.113 2910.954 3847.645 5562.367 5601.902
## [49] 5602.135 5425.341 4170.469 2834.891 4212.502 5493.079 5315.299 5379.269
## [57] 5043.173 4156.969 3094.537 4152.062 5743.599 6116.843 5745.467 5534.264
## [65] 4578.557 3071.252 4375.717 5876.830 6008.955 5785.745 5683.216 4505.411
## [73] 3229.939 4477.531 6048.338 6032.351 6472.559 5821.047 4294.172 2988.317
## [81] 4143.791 6012.203 6051.720 5917.039 5376.331 4371.571 3158.859 4450.246
## [89] 6110.452 6098.221 6336.015 6530.471 4876.664 3382.521 4460.187 6108.389
## [97] 6688.388 6231.055 6059.221 5012.627 3630.046 4909.253 6842.912 6951.352
## [105] 6947.234 6478.916 4651.798 3570.714 4700.247 6500.222 6008.527 5561.052
## [113] 4696.695 4273.243 3564.461 4811.358 6820.401 7254.227 7208.896 6778.355
## [121] 5490.770 3984.950 5253.690 7627.471 6977.024 7102.154 6582.290 4723.553
## [129] 3232.818 3870.954 5322.841 5037.074 4949.924 4134.570 2788.413 1741.673
## [137] 2090.813 2592.803 1822.865 1571.803 2056.296 2006.519 1859.766 1979.456
## [145] 2645.436
##
## $call
## forecast.dlm(model = model.dlm2, x = test$Xt, h = 145)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"#akurasi testing
mape.dlm2 <- MAPE(fore.dlm2$forecasts, test$Yt)
#akurasi data training
mape_train <- GoF(model.dlm2)["MAPE"]
c("MAPE_testing" = mape.dlm2, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.05433837
##
## $MAPE_training.MAPE
## [1] 0.02309118Output di samping menunjukkan bahwa nilai MAPE Testing lebih besar daripada nilai MAPE Training dengan selisih yang tidak terlalu jauh yaitu 0.03124719. Hal tersebut mengindikasikan bahwa model ini tidak underfitted atau overfitted.
4. Model Autoregressive / Dynamic Regression
#library(dLagM)
model.ardl = ardlDlm(x = train$Xt, y = train$Yt, p = 1 , q = 1) #p:lag x, q:lag y
#model untuk p=1, q=1: yt=b0+b1yt-1+b2xt+b3xt-1
summary(model.ardl)##
## Time series regression with "ts" data:
## Start = 2, End = 584
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -371.66 -81.73 -7.76 73.48 737.59
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.024109 17.963494 2.896 0.00392 **
## X.t 1.458653 0.007938 183.762 < 2e-16 ***
## X.1 -0.426583 0.058810 -7.254 1.31e-12 ***
## Y.1 0.282603 0.040037 7.059 4.83e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 123 on 579 degrees of freedom
## Multiple R-squared: 0.9925, Adjusted R-squared: 0.9925
## F-statistic: 2.562e+04 on 3 and 579 DF, p-value: < 2.2e-16Berdasarkan ouput, terlihat bahwa semua parameter signifikan. Menurut model Autoregressive, pada taraf nyata 5% jumlah halaman yang dimuat dipengaruhi jumlah pengunjung halaman website pada waktu sekarang, jumlah pengunjung halaman website periode sebelumnya, dan lag dari jumlah halaman yang dimuat.
AIC(model.ardl)## [1] 7271.844BIC(model.ardl)## [1] 7293.685#ramalan
(fore.ardl <- forecast(model = model.ardl, x=test$Xt, h=145))## $forecasts
## [1] 3681.051 2968.675 1868.933 2231.459 3796.486 3917.125 3692.747 3595.595
## [9] 2972.417 1890.728 2552.371 4191.968 4106.433 3529.059 2880.062 3156.155
## [17] 2042.736 2653.070 3996.564 4509.076 4572.850 4360.956 3078.460 2047.021
## [25] 2438.758 3353.555 4291.277 4902.273 4498.954 3632.254 2303.554 3124.413
## [33] 4826.082 4755.791 4917.570 4901.531 3663.900 2520.631 3484.977 5626.936
## [41] 5569.321 5433.457 5103.372 4117.529 2920.386 3877.571 5570.374 5584.312
## [49] 5596.370 5419.620 4173.712 2847.300 4249.776 5493.849 5303.954 5372.953
## [57] 5037.518 4165.219 3101.086 4182.058 5749.172 6106.633 5733.358 5539.055
## [65] 4585.889 3084.769 4418.194 5882.026 6003.859 5777.507 5688.976 4511.386
## [73] 3248.412 4515.851 6058.573 6023.588 6476.838 5813.127 4309.702 3007.295
## [81] 4187.409 6027.581 6036.737 5912.642 5372.424 4387.167 3173.112 4490.924
## [89] 6117.553 6087.374 6337.522 6534.857 4880.560 3408.117 4503.601 6130.603
## [97] 6689.220 6220.569 6071.832 5023.507 3654.501 4956.237 6862.045 6947.436
## [105] 6953.555 6486.135 4669.721 3609.370 4745.629 6525.974 5993.094 5566.731
## [113] 4695.610 4302.790 3574.354 4845.229 6827.908 7244.560 7210.567 6787.694
## [121] 5514.723 4017.073 5313.554 7660.523 6965.755 7127.124 6587.134 4752.578
## [129] 3266.719 3917.183 5347.661 5024.206 4948.840 4119.324 2792.958 1745.089
## [137] 2101.263 2579.077 1785.485 1544.024 2026.744 1970.731 1824.393 1942.388
## [145] 2616.964
##
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 145)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"#akurasi testing
mape.ardl <- MAPE(fore.ardl$forecasts, test$Yt) #data testing
#akurasi data training
mape_train <- GoF(model.ardl)["MAPE"]
c("MAPE_testing" = mape.ardl, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.05475761
##
## $MAPE_training.MAPE
## [1] 0.02273106Output tersebut menunjukkan bahwa nilai MAPE Testing lebih besar daripada nilai MAPE Training dengan selisih yang tidak terlalu jauh. Hal tersebut mengindikasikan bahwa model ini tidak underfitted atau overfitted.
#PEMODELAN DLM dan ARDL dengan library dynlm
#library(dynlm)
#sama dengan model dlm p=1
cons_lm1 <- dynlm(Yt ~ Xt+L(Xt),data = train.ts)
#sama dengan model ardl p=0 q=1
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 = 584
##
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -428.74 -93.77 -12.50 79.70 720.06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 74.165745 18.416909 4.027 6.4e-05 ***
## Xt 1.453362 0.008228 176.634 < 2e-16 ***
## L(Xt) -0.015205 0.008202 -1.854 0.0643 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 128.1 on 580 degrees of freedom
## Multiple R-squared: 0.9919, Adjusted R-squared: 0.9919
## F-statistic: 3.542e+04 on 2 and 580 DF, p-value: < 2.2e-16summary(cons_lm2)##
## Time series regression with "ts" data:
## Start = 2, End = 584
##
## Call:
## dynlm(formula = Yt ~ Xt + L(Yt), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -429.97 -93.66 -12.57 78.54 724.35
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.869216 18.588430 3.705 0.000232 ***
## Xt 1.447486 0.008126 178.128 < 2e-16 ***
## L(Yt) -0.005188 0.005596 -0.927 0.354228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 128.4 on 580 degrees of freedom
## Multiple R-squared: 0.9918, Adjusted R-squared: 0.9918
## F-statistic: 3.526e+04 on 2 and 580 DF, p-value: < 2.2e-16summary(cons_lm3)##
## Time series regression with "ts" data:
## Start = 2, End = 584
##
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -371.66 -81.73 -7.76 73.48 737.59
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.024109 17.963494 2.896 0.00392 **
## Xt 1.458653 0.007938 183.762 < 2e-16 ***
## L(Xt) -0.426583 0.058810 -7.254 1.31e-12 ***
## L(Yt) 0.282603 0.040037 7.059 4.83e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 123 on 579 degrees of freedom
## Multiple R-squared: 0.9925, Adjusted R-squared: 0.9925
## F-statistic: 2.562e+04 on 3 and 579 DF, p-value: < 2.2e-16summary(cons_lm4)##
## Time series regression with "ts" data:
## Start = 3, End = 584
##
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -454.82 -86.89 -7.81 79.21 742.44
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 112.508454 20.913286 5.380 1.08e-07 ***
## Xt 1.435912 0.009261 155.052 < 2e-16 ***
## L(Xt) 0.025350 0.013118 1.933 0.053786 .
## L(Xt, 2) -0.035758 0.009224 -3.876 0.000118 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 126.6 on 578 degrees of freedom
## Multiple R-squared: 0.9921, Adjusted R-squared: 0.992
## F-statistic: 2.417e+04 on 3 and 578 DF, p-value: < 2.2e-16#SSE
deviance(cons_lm1)## [1] 9519492deviance(cons_lm2)## [1] 9561731deviance(cons_lm3)## [1] 8765232deviance(cons_lm4)## [1] 9263424Diperoleh hasil bahwa model ketiga dengan p=1 q=1 merupakan model DLM dan ARDL yang memiliki nilai SSE terbaik (nilainya terendah) diantara keempat model. Selain itu, pada model ketiga ini diperoleh bahwa semua peubah signifikan.
5. Uji Diagnostik Model
Uji Non Autokorelasi
#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 579 -1 49.824 4.830e-12 ***
## M2 vs. ME 579 -1 52.614 1.309e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Durbin Watson
Hipotesis : H0: tidak ada autokorelasi H1: ada autokorelasi
dwtest(cons_lm1)##
## Durbin-Watson test
##
## data: cons_lm1
## DW = 1.4549, p-value = 1.675e-11
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm2)##
## Durbin-Watson test
##
## data: cons_lm2
## DW = 1.441, p-value = 5.273e-12
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm3)##
## Durbin-Watson test
##
## data: cons_lm3
## DW = 2.0655, p-value = 0.7683
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm4)##
## Durbin-Watson test
##
## data: cons_lm4
## DW = 1.487, p-value = 2.187e-10
## alternative hypothesis: true autocorrelation is greater than 0p-value > 5% : Tak Tolak H0 sehingga tidak cukup bukti untuk menyatakan bahwa terdapat autokorelasi (positif) pada model ke 3 pada taraf 5%. Sedangkan pada model 1, 2, dan 4 terdeteksi adanya autokorelasi positif padda taraf nyata 5%.
Uji Heterogenitas
Hipotesis : H0: tidak ada heteroskedastisitas H1: ada heteroskedastisitas
#heterogenitas
bptest(cons_lm1)##
## studentized Breusch-Pagan test
##
## data: cons_lm1
## BP = 11.269, df = 2, p-value = 0.003573bptest(cons_lm2)##
## studentized Breusch-Pagan test
##
## data: cons_lm2
## BP = 10.297, df = 2, p-value = 0.005809bptest(cons_lm3)##
## studentized Breusch-Pagan test
##
## data: cons_lm3
## BP = 12.599, df = 3, p-value = 0.005588bptest(cons_lm4)##
## studentized Breusch-Pagan test
##
## data: cons_lm4
## BP = 12.994, df = 3, p-value = 0.00465Untuk semua model, p-value > 5% : Tak Tolak H0 sehingga tidak cukup bukti untuk menyatakan bahwa terdapat heteroskedastisitas pada taraf 5%.
Uji Normalitas
Hipotesis : H0: Data berdistribusi normal H1: Data tidak berdistribusi normal
#shapiro wilk
shapiro.test(residuals(cons_lm1))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm1)
## W = 0.98005, p-value = 3.844e-07shapiro.test(residuals(cons_lm2))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm2)
## W = 0.97956, p-value = 2.832e-07shapiro.test(residuals(cons_lm3))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm3)
## W = 0.98086, p-value = 6.398e-07shapiro.test(residuals(cons_lm4))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm4)
## W = 0.98112, p-value = 7.734e-075. Perbandingan Keakuratan Ramalan
#PERBANDINGAN
akurasi <- matrix(c(mape.koyck,mape.dlm2, mape.ardl))
row.names(akurasi)<- c("Koyck","DLM","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi## MAPE
## Koyck 0.07187000
## DLM 0.05433837
## Autoregressive 0.05475761MAPE (Mean Absolute Percent Error) menggambarkan seberapa besar kesalahan peramalan dibandingkan dengan nilai sebenarnya dari series tersebut. Dapat diketahui dari output tersebut, bahwa pemodelan dengan DLM pada data memberikan hasil peramalan terbaik dibandingkan pemodelan lain dengan nilai MAPE terkecil yaitu 0.0543.
#PLOT
par(mfrow=c(1,1))
plot(test$Xt, test$Yt, type="b", col="black")
points(test$Xt, fore.koyck$forecasts,col="red")
points(test$Xt, fore.dlm2$forecasts,col="orange")
points(test$Xt, fore.ardl$forecasts,col="green")
legend("topleft",c("aktual", "koyck","DLM", "autoregressive"), lty=1, col=c("black","red","orange","green"), cex=0.8)
Dari plot gabungan di atas terlihat bahwa data tersebut tidak cocok dianalisis dengan pemodelan Koyck, DLM, maupun Autoregressive.
Simpulan
Berdasarkan analisis yang telah dilakukan, diketahui bahwa model regresi distributed lag memiliki nilai MAPE yang lebih kecil dibandingkan dengan model lainnya.