Regresi dengan Peubah Lag

Package

Terdapat 6 packages yang digunakan, di antaranya:

readxl DLagM dnylm MLmetrics lmtest car

Data

Data yang digunakan merupakan data pembelian bersih sebuah grocery store periode 1 Januari 2018 sampai dengan 24 Desember 2018 yang diperoleh dari situs web Kaggle. Data terdiri atas dua peubah, yaitu Net Purchases sebagai peubah respon (Y) dan Gross Sales sebagai peubah penjelas (X).

purchases <- read_excel("C:\\Users\\ASUS\\Downloads\\purchases.xlsx", sheet = 3)
str(purchases)

## tibble [355 x 3] (S3: tbl_df/tbl/data.frame)
##  $ t : num [1:355] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Yt: num [1:355] 1218 1995 2161 1990 3706 ...
##  $ Xt: num [1:355] 1770 2758 2972 2705 5087 ...

purchases

## # A tibble: 355 x 3
##        t    Yt    Xt
##    <dbl> <dbl> <dbl>
##  1     1 1218. 1770.
##  2     2 1995. 2758.
##  3     3 2161. 2972.
##  4     4 1990. 2705.
##  5     5 3706. 5087.
##  6     6 1434. 2097.
##  7     7 1129. 1548.
##  8     8 2428. 3269.
##  9     9 2205. 3079.
## 10    10 2430. 3258.
## # ... with 345 more rows

Time Series Plot

purchases1 <- purchases[,-1]
purchases1 <- purchases1[,-2]
purchases1.ts <- ts(purchases1)
ts.plot(purchases1.ts, 
        xlab="Time Period", 
        ylab="Net Purchase", 
        main= "Time Series Plot of Net Purchases\n Grocery Shop 2018")
points(purchases1.ts)

Pembelian bersih didefinisikan sebagai seluruh pembelian barang dagangan yang dilakukan oleh perusahaan, baik pembelian barang dagangan secara tunai maupun pembelian barang dagangan secara kredit, ditambah dengan biaya angkut pembelian tersebut serta dikurangi dengan potongan pembelian dan retur pembelian yang terjadi. Pada plot di atas, terlihat data untuk pembelian bersih berpola stasioner.

purchases2 <- purchases[,-1:-2]
purchases2.ts <- ts(purchases2)
ts.plot(purchases2.ts, 
        xlab="Time Period", 
        ylab="Gross purchases", 
        main= "Time Series Plot of Gross purchases\n Grocery Shop 2018")
points(purchases2.ts)

Penjualan kotor didefinisikan sebagai jumlah keseluruhan dari transaksi penjualan dalam jangka waktu tertentu untuk suatu perusahaan sebelum dikurangi dengan retur, diskon, komisi penjualan, dan lainnya. Pada plot di atas, terlihat data untuk penjualan kotor berpola stasioner.

Scatterplot

# Diagram pencar identifikasi model
plot(purchases$Xt, purchases$Yt, pch = 20, 
     col = "steelblue", 
     main = "Scatterplot Net Purchases \nvs Gross purchases 2018",
     ylab = "Net Purchases", xlab = "Gross purchases")

Berdasarkan scatterplot di atas, terlihat hubungan antara net purchases dan gross purchases merupakan linear positif. Artinya, peningkatan gross purchases akan meningkatkan net purchases pada grocery store tersebut.

# Korelasi
cor(purchases$Xt, purchases$Yt)

## [1] 0.9821884

Korelasi antara kedua peubah mendekati 1, yaitu 0.9821884. Artinya, terdapat hubungan linear positif yang kuat antara kedua peubah tersebut.

Splitting Data

Data training adalah data yang digunakan dalam proses konstruksi model atau digunakan untuk melatih dan membangun model, sedangkan data training digunakan untuk menguji performa model yang didapatkan pada saat tahapan training. Pembagian dataset untuk training dan testing umumnya memiliki rasio 90:10, 80:20, dan 70:30. Pada data ini, pembagian dataset untuk training dan testing adalah 80:20 (80% training 20% testing).

#SPLIT purchases 80:20
train<-purchases[1:284,]
test<-purchases[285:355,]

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

MODEL KOYCK

Metode Koyck didasari asumsi bahwa semakin jauh jarak lag pada variabel bebas dari periode sekarang maka semakin kecil pengaruh variabel lag terhadap variabel tak bebas (Aqibah et al. 2020).

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 
## -966.79 -184.35   11.23  174.94  810.01 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.374e+03  3.293e+02   4.171 4.05e-05 ***
## Y.1         3.805e-03  4.034e-02   0.094    0.925    
## X.t         2.978e-01  1.076e-01   2.769    0.006 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 318.8 on 280 degrees of freedom
## Multiple R-Squared: 0.6643,  Adjusted R-squared: 0.6619 
## Wald test: 5.448 on 2 and 280 DF,  p-value: 0.004774 
## 
## Diagnostic tests:
## NULL
## 
##                            alpha      beta         phi
## Geometric coefficients:  1378.99 0.2977916 0.003805469

Berdasarkan output, terlihat bahwa hanya Y.1 yang tidak signifikan pada model. Persamaan regresinya adalah: Yt = 1374- 0.003805Xt + 0.2978Yt-1

AIC(model.koyck)

## [1] 4070.885

BIC(model.koyck)

## [1] 4085.467

#ramalan
(fore.koyck <- forecast(model = model.koyck, x=test$Xt, h=71))

## $forecasts
##  [1] 2413.311 2413.830 2406.959 2214.765 2282.564 2422.555 2657.185 2268.119
##  [9] 2106.936 2253.465 2297.315 2292.863 2791.585 1950.182 2265.232 2346.704
## [17] 2905.826 1558.803 2468.296 2621.537 1761.533 2274.821 2250.878 2223.944
## [25] 2303.236 2670.902 2981.379 2411.811 2512.563 2259.678 2298.438 2233.863
## [33] 2505.853 2749.498 2226.402 2165.148 2234.642 2116.594 2114.793 2529.815
## [41] 2733.696 1883.232 2157.570 2137.000 2269.457 2204.971 2331.249 2570.358
## [49] 1935.126 2191.644 2185.072 2316.036 2335.629 2350.286 2591.208 2287.386
## [57] 2337.593 2238.220 2227.562 2202.707 2372.404 2466.770 1833.641 2219.572
## [65] 2307.984 2433.587 2505.767 2562.917 3297.345 2319.748 3026.043
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 71)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

#mape purchases testing
mape.koyck <- MAPE(fore.koyck$forecasts, test$Yt)

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

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

## $MAPE_testing
## [1] 0.1884594
## 
## $MAPE_training.MAPE
## [1] 0.1170681

Nilai MAPE testing tidak berbeda signifikan dengan MAPE training sehingga tidak ada indikasi bahwa model overfitted/underfitted.

Regression with Distributed Lag

Model regresi linear yang sering ditemui biasanya tidak memerhatikan pengaruh waktu, karena pada umumnya model regresi linear cenderung mengasumsikan bahwa pengaruh variabel bebas terhadap variabel takbebas terjadi dalam kurun waktu yang sama. Namun, dalam model regresi linear juga terdapat model regresi yang memerhatikan pengaruh waktu. Waktu yang diperlukan bagi variabel bebas x dalam memengaruhi variabel takbebas Y disebut lag.

Regression with Distributed Lag (lag=2)

model.dlm = dLagM::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 
## -154.60  -65.42  -22.54   25.01  486.35 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 356.738566  40.676270   8.770  < 2e-16 ***
## x.t           0.676587   0.007925  85.378  < 2e-16 ***
## x.1          -0.061524   0.008016  -7.675 2.81e-13 ***
## x.2          -0.017992   0.007895  -2.279   0.0234 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 105.2 on 278 degrees of freedom
## Multiple R-squared:  0.9636, Adjusted R-squared:  0.9632 
## F-statistic:  2455 on 3 and 278 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 3432.201 3450.411

DLM dengan lag = 2 menghasilkan parameter regresi yang seluruhnya signifikan pada alpha = 5%. Persamaan regresinya adalah: Yt = 356.738566 + 0.676587Xt 0.061524Xt-1 - 0.017992Xt-2 Berdasarkan persamaan model lag di atas, koefisien β bernilai positif pada lag 0, sedangkan pada lag 1 dan lag 2 bernilai negatif. Koefisien β yang positif berarti semakin meningkatnya pembelian bersih dari sebuah grocery store, maka dugaan rata-rata penjualan semakin meningkat secara signifikan. Sedangkan pada 1 dan 2 periode sebelumnya berlaku sebaliknya.

AIC(model.dlm)

## [1] 3432.201

BIC(model.dlm)

## [1] 3450.411

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

## $forecasts
##  [1] 2346.3290 2403.4917 2408.2792 1972.8946 2168.7146 2483.6417 2982.5045
##  [8] 2039.7068 1743.2812 2134.1448 2211.7566 2183.4284 3314.9084 1296.1549
## [15] 2163.3230 2330.9039 3564.7207  379.1286 2702.1293 2934.9861  893.5123
## [22] 2235.8951 2122.3331 2035.4858 2222.9811 3042.8584 3664.3940 2281.5893
## [29] 2614.6673 2052.4573 2188.8142 2048.8602 2678.3543 3177.2367 1920.0749
## [36] 1879.0322 2081.3596 1801.7263 1818.8851 2769.2655 3143.2351 1142.3330
## [43] 1936.6297 1881.6069 2170.4235 1996.6893 2289.5703 2809.5270 1307.2564
## [50] 2012.5716 1980.3588 2263.8844 2280.6608 2301.9343 2845.0186 2102.0028
## [57] 2267.1068 2048.6968 2042.8076 1994.5668 2386.0844 2565.4403 1096.5219
## [64] 2104.0550 2259.6347 2502.8158 2634.5916 2741.4309 4393.4877 2010.8966
## [71] 3782.2519
## 
## $call
## forecast.dlm(model = model.dlm, x = test$Xt, h = 71)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

#mape purchases testing
mape.dlm <- MAPE(fore.dlm$forecasts, test$Yt)

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

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

## $MAPE_testing
## [1] 0.03547406
## 
## $MAPE_training.MAPE
## [1] 0.03120183

Nilai MAPE testing tidak berbeda signifikan dengan MAPE training sehingga tidak ada indikasi bahwa model overfitted/underfitted.

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 max lag maksimum

##   q - k    MASE      AIC      BIC   GMRAE    MBRAE R.Adj.Sq Ljung-Box
## 4     4 0.13595 3400.163 3425.607 0.14455  8.96184  0.96465         0
## 3     3 0.13410 3419.542 3441.372 0.14530 -0.60622  0.96354         0
## 2     2 0.13409 3432.201 3450.411 0.13874 -2.00282  0.96324         0
## 1     1 0.13529 3446.934 3461.516 0.14164 -0.05286  0.96271         0

Dari output di atas, terlihat bahwa persamaan optimum pada saat lag 4 (q = 4), dengan nilai AIC dan BIC yang terkecil dan nilai adjusted R-squared yang terbesar.

Regression with Distributed Lag (lag=4)

#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 
## -167.24  -64.78  -18.63   30.09  470.79 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 449.915387  48.270405   9.321  < 2e-16 ***
## x.t           0.678045   0.007807  86.849  < 2e-16 ***
## x.1          -0.059194   0.007918  -7.476 1.04e-12 ***
## x.2          -0.013531   0.007936  -1.705  0.08933 .  
## x.3          -0.010704   0.007915  -1.352  0.17738    
## x.4          -0.024835   0.007777  -3.193  0.00157 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 103.4 on 274 degrees of freedom
## Multiple R-squared:  0.9653, Adjusted R-squared:  0.9646 
## F-statistic:  1524 on 5 and 274 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 3400.163 3425.607

DLM dengan lag = 4 menghasilkan parameter regresi yang tidak seluruhnya signifikan pada alpha = 5%. Persamaan regresinya adalah: Yt = 449.915387 + 0.678045Xt - 0.059194Xt-1 - 0.013531Xt-2 - 0.010704Xt-3 - 0.024835Xt-4 Berdasarkan persamaan model lag di atas, koefisien β bernilai positif pada lag 0 dan bernilai negatif pada lag 1 sampai 4.

AIC(model.dlm2)

## [1] 3400.163

BIC(model.dlm2)

## [1] 3425.607

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

## $forecasts
##  [1] 2356.7056 2419.2082 2396.4040 1941.8672 2165.3879 2478.6326 2988.2419
##  [8] 2061.0230 1753.5337 2117.9581 2189.0353 2201.4891 3344.0521 1311.1766
## [15] 2177.3316 2317.6195 3548.2877  420.5549 2716.4215 2910.1744  881.2561
## [22] 2301.3147 2097.1156 2036.1165 2276.6036 3056.3497 3686.5863 2308.2943
## [29] 2622.2143 2009.2727 2140.1247 2040.2693 2671.7055 3192.6294 1937.8930
## [36] 1891.7066 2054.6371 1772.6863 1835.7892 2788.8539 3165.4784 1178.1255
## [43] 1955.3455 1847.7526 2154.9348 2042.6080 2315.2135 2833.3928 1323.0543
## [50] 2029.1065 1970.2235 2261.1762 2322.7020 2324.8356 2865.4200 2111.4052
## [57] 2275.8681 2042.9332 2027.8510 2001.4846 2392.7805 2582.2246 1115.2316
## [64] 2117.0793 2249.1039 2514.3514 2687.6318 2761.6758 4407.0917 2013.2744
## [71] 3783.4093
## 
## $call
## forecast.dlm(model = model.dlm2, x = test$Xt, h = 71)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.0375484
## 
## $MAPE_training.MAPE
## [1] 0.03171937

Nilai MAPE testing tidak berbeda signifikan dengan MAPE training sehingga tidak ada indikasi bahwa model overfitted/underfitted.

Model Autoregressive

Apabila peubah tak bebas dipengaruhi oleh peubah bebas pada waktu sekarang, serta dipengaruhi oleh peubah tak bebas itu sendiri pada satu waktu yang lalu, maka model tersebut disebut autoregressive (Gujarati 2004).

# 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
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 = 284
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -270.04  -43.40  -15.87   23.03  518.01 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 206.166505  31.034648   6.643  1.6e-10 ***
## X.t           0.689037   0.006767 101.824  < 2e-16 ***
## X.1          -0.401003   0.031427 -12.760  < 2e-16 ***
## Y.1           0.503264   0.045865  10.973  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 88.65 on 279 degrees of freedom
## Multiple R-squared:  0.9741, Adjusted R-squared:  0.9739 
## F-statistic:  3503 on 3 and 279 DF,  p-value: < 2.2e-16

Berdasarkan model di atas, seluruh parameter regresi signifikan pada alpha = 5%. Persamaan regresinya adalah: Yt = 206.166505 + 0.689037Xt - 0.401003Xt-1 + 0.503264Yt-1 Berdasarkan persamaan model lag di atas, koefisien β pada lag 0 bernilai positif, sedangkan pada lag 1 bernilai negatif. Koefisien γ bernilai positif. Pembelian bersih pada periode saat ini terpengaruh positif oleh penjualan pada periode saat ini, sedangkan pada satu periode sebelumnya terpengaruh negatif. Selain itu, terdapat pula pengaruh positif oleh pembelian bersih pada satu periode sebelumnya.

AIC(model.ardl)

## [1] 3347.406

BIC(model.ardl)

## [1] 3365.633

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

## $forecasts
##  [1] 2305.6467 2368.0119 2379.1241 1949.3291 2150.3681 2482.5779 3003.2636
##  [8] 2047.7792 1722.5080 2114.3248 2213.5447 2194.4938 3345.1183 1301.3403
## [15] 2144.7417 2326.3763 3602.6842  370.8170 2677.3541 2953.0998  898.9316
## [22] 2219.2327 2128.1783 2055.1152 2238.1958 3073.4393 3714.2496 2299.9298
## [29] 2594.8580 2018.6748 2161.6623 2030.3777 2681.3648 3203.7530 1927.4562
## [36] 1853.6616 2057.6651 1792.6834 1815.5180 2789.1346 3188.3375 1147.2048
## [43] 1908.5057 1867.8518 2183.1599 2012.9972 2307.6266 2837.6750 1311.1723
## [50] 1998.6926 1978.5497 2281.6647 2302.0034 2320.2690 2867.1452 2112.9085
## [57] 2262.5319 2038.2943 2035.7293 1990.8636 2394.5683 2585.9522 1090.2832
## [64] 2089.1692 2270.1080 2533.9346 2663.9335 2764.4030 4437.2131 2021.9421
## [71] 3769.4599
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 71)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.03483691
## 
## $MAPE_training.MAPE
## [1] 0.02430779

Nilai MAPE testing tidak berbeda signifikan dengan MAPE training sehingga tidak ada indikasi bahwa model overfitted/underfitted.

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

## $p
##   Xt
## 1 14
## 
## $q
## [1] 15
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  4098.660 4081.777 4072.008 4061.799 4052.069 4040.662 4023.267 4010.904
## p = 2  4086.103 4050.826 4040.046 4028.921 4019.278 4007.101 3988.972 3975.998
## p = 3  4040.977 4040.977 4032.349 4016.543 4007.457 3995.865 3978.489 3964.307
## p = 4  4040.757 4017.325 4017.325 4017.582 4008.885 3997.290 3979.776 3965.645
## p = 5  4047.922 4019.494 4008.674 4008.674 4008.318 3995.804 3978.744 3964.633
## p = 6  4030.349 4007.103 3997.448 3998.999 3998.999 3918.054 3905.885 3892.009
## p = 7  3911.520 3902.405 3902.735 3904.371 3906.338 3906.338 3906.260 3891.204
## p = 8  4014.942 3981.781 3970.511 3971.747 3971.935 3893.294 3893.294 3884.100
## p = 9  4004.811 3971.005 3958.503 3959.666 3959.117 3873.793 3874.274 3874.274
## p = 10 3974.810 3960.181 3949.745 3951.344 3951.934 3874.668 3874.187 3867.888
## p = 11 3967.282 3942.957 3939.357 3941.156 3941.751 3866.696 3865.391 3857.072
## p = 12 3975.255 3944.702 3933.428 3935.017 3935.166 3856.023 3855.158 3848.347
## p = 13 3962.498 3935.272 3924.347 3925.979 3926.578 3847.538 3846.591 3839.491
## p = 14 3885.394 3868.225 3863.822 3865.785 3867.537 3817.202 3815.222 3807.276
## p = 15 3925.523 3889.448 3874.409 3875.424 3875.442 3788.681 3790.604 3786.177
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  3997.404 3988.840 3979.202 3969.832 3960.350 3933.911 3917.647
## p = 2  3964.188 3955.650 3946.842 3937.840 3928.346 3900.151 3881.207
## p = 3  3952.961 3944.584 3935.637 3926.829 3917.817 3889.151 3868.220
## p = 4  3954.547 3946.175 3937.289 3928.496 3919.371 3890.717 3869.584
## p = 5  3953.586 3945.232 3936.455 3927.642 3918.381 3890.063 3868.732
## p = 6  3880.606 3872.361 3862.702 3853.270 3844.862 3814.652 3791.724
## p = 7  3879.698 3871.435 3861.882 3852.320 3843.927 3814.134 3792.059
## p = 8  3874.374 3866.080 3856.244 3846.398 3837.574 3808.927 3786.862
## p = 9  3876.214 3867.888 3858.055 3848.235 3839.404 3810.587 3788.611
## p = 10 3867.888 3868.037 3858.750 3848.839 3839.920 3811.510 3789.740
## p = 11 3858.802 3858.802 3860.570 3850.728 3841.794 3813.447 3791.603
## p = 12 3850.285 3850.835 3850.835 3851.940 3843.252 3815.113 3793.315
## p = 13 3841.462 3841.566 3843.521 3843.521 3841.275 3811.312 3788.628
## p = 14 3808.366 3810.123 3811.455 3812.205 3812.205 3806.158 3784.549
## p = 15 3788.129 3788.759 3790.752 3792.733 3784.797 3784.797 3786.467
## 
## $min.Stat
## [1] 3784.549

#PEMODELAN DLM dan ARDL dengan 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 = 284
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -175.90  -62.86  -24.30   23.88  476.27 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 308.474562  35.353477   8.725 2.42e-16 ***
## Xt            0.675052   0.007937  85.048  < 2e-16 ***
## L(Xt)        -0.063910   0.007910  -8.080 1.95e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 105.9 on 280 degrees of freedom
## Multiple R-squared:  0.963,  Adjusted R-squared:  0.9627 
## F-statistic:  3642 on 2 and 280 DF,  p-value: < 2.2e-16

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 284
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -156.03  -65.27  -28.28   21.80  503.27 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 272.445367  38.434390   7.089 1.10e-11 ***
## Xt            0.669981   0.008291  80.811  < 2e-16 ***
## L(Yt)        -0.068828   0.012141  -5.669 3.57e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 111.4 on 280 degrees of freedom
## Multiple R-squared:  0.959,  Adjusted R-squared:  0.9588 
## F-statistic:  3278 on 2 and 280 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 284
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -270.04  -43.40  -15.87   23.03  518.01 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 206.166505  31.034648   6.643  1.6e-10 ***
## Xt            0.689037   0.006767 101.824  < 2e-16 ***
## L(Xt)        -0.401003   0.031427 -12.760  < 2e-16 ***
## L(Yt)         0.503264   0.045865  10.973  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 88.65 on 279 degrees of freedom
## Multiple R-squared:  0.9741, Adjusted R-squared:  0.9739 
## F-statistic:  3503 on 3 and 279 DF,  p-value: < 2.2e-16

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 284
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -154.60  -65.42  -22.54   25.01  486.35 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 356.738566  40.676270   8.770  < 2e-16 ***
## Xt            0.676587   0.007925  85.378  < 2e-16 ***
## L(Xt)        -0.061524   0.008016  -7.675 2.81e-13 ***
## L(Xt, 2)     -0.017992   0.007895  -2.279   0.0234 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 105.2 on 278 degrees of freedom
## Multiple R-squared:  0.9636, Adjusted R-squared:  0.9632 
## F-statistic:  2455 on 3 and 278 DF,  p-value: < 2.2e-16

#SSE
deviance(cons_lm1)

## [1] 3138578

deviance(cons_lm2)

## [1] 3471889

deviance(cons_lm3)

## [1] 2192433

deviance(cons_lm4)

## [1] 3077091

Uji 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    279 -1 120.40 < 2.2e-16 ***
## M2 vs. ME    279 -1 162.82 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan hasil pengujian, terlihat bahwa peubah lag bagi Yt dan Xt berpengaruh signifikan.

Uji Non-Autokorelasi (Durbin-Watson Test)

Hipotesis

H0 : Tidak ada autokorelasi

H1 : Ada autokorelasi

dwtest(cons_lm1)

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

dwtest(cons_lm2)

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

dwtest(cons_lm3)

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

dwtest(cons_lm4)

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

Di antara keempat model, H0 tidak ditolak pada model 3. Oleh karena itu, tidak terjadi autokorelasi pada model autoregressive dengan p = 0 dan q = 1.

Uji Heteroskedastisitas (Breusch-Pagan Test)

Hipotesis

H0 : Tidak terdapat gejala heteroskedastisitas

H1 : Terdapat gejala heteroskedastisitas

bptest(cons_lm1)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 40.262, df = 2, p-value = 1.808e-09

bptest(cons_lm2)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 26.779, df = 2, p-value = 1.531e-06

bptest(cons_lm3)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm3
## BP = 29.465, df = 3, p-value = 1.788e-06

bptest(cons_lm4)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm4
## BP = 38.334, df = 3, p-value = 2.402e-08

Pada keempat model di atas, p-value bernilai < 0.05 sehingga H0 ditolak. Artinya, terdapat gejala heteroskedastisitas atau ragam menyebar secara tidak homogen.

Uji Normalitas (Shapiro-Wilk Normality Test)

Hipotesis

H0 : Sisaan menyebar normal

H1 : Sisaan tidak menyebar normal

shapiro.test(residuals(cons_lm1))

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

shapiro.test(residuals(cons_lm2))

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

shapiro.test(residuals(cons_lm3))

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

shapiro.test(residuals(cons_lm4))

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

Pada keempat model di atas, p-value bernilai < 0.05 sehingga H0 ditolak. Artinya, sisaan menyebar tidak normal.

Perbandingan Keakuratan Ramalan

Metrik Akurasi

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.18845941
## DLM 1          0.03547406
## DLM 2          0.03754840
## Autoregressive 0.03483691

Nilai MAPE terkecil dimiliki oleh model autoregressive dengan p = 1 dan q = 1. Nilai MAPE pada model autoregressive tersebut tidak memiliki perbedaan signifikan dengan model DLM 1 (p = 2) dan DLM 2 (p = 4). Sementara itu, model KOYCK memiliki MAPE yang berbeda signifikan dengan ketiga model lainnya. Berdasarkan perbandingan nilai MAPE, metode terbaik untuk memprediksi data testing adalah model autoregressive. Sebaliknya, metode KOYCK memprediksi dengan tidak tepat atau tidak mendekati data aktual.

Plot Perbandingan

par(mfrow=c(2,2))

plot(test$Xt, test$Yt, type="b", col="black", main = "Model KOYCK")
points(test$Xt, fore.koyck$forecasts,col="red")
lines(test$Xt, fore.koyck$forecasts,col="red")

plot(test$Xt, test$Yt, type="b", col="black", main = "DLM Lag 2")
points(test$Xt, fore.dlm$forecasts,col="blue")
lines(test$Xt, fore.dlm$forecasts,col="blue")

plot(test$Xt, test$Yt, type="b", col="black", main = "DLM Lag 4")
points(test$Xt, fore.dlm2$forecasts,col="orange")
lines(test$Xt, fore.dlm2$forecasts,col="orange")

plot(test$Xt, test$Yt, type="b", col="black", main = "Autoregressive Linear Model")
points(test$Xt, fore.ardl$forecasts,col="green")
lines(test$Xt, fore.ardl$forecasts,col="green")

Pada keempat plot di atas, hasil peramalan dari model DLM lag 2, DLM lag 4, dan autoregressive mendekati data aktual. Sementara itu, garis merah pada model KOYCK menyimpang dari data aktual. Hal ini sejalan dengan pernyataan sebelumnya.

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)