Tugas MPDW Pertemuan 3

Pada Tugas Kali Ini, Saya akan membuat persamaan regresi dengan Peubah Lag dalam rangka mencari hubungan antara Jumlah Produksi Ikan di Indonesia dengan Tingkat Konsumsi Ikan per-Kapitanya dari jangka waktu 1961 sampai 2021. Produksi Ikan disini berarti perhitungan seluruh Ikan yang ada baik yang ditangkap ataupun yang dibudidayakan secara konvensional. Data ini diambil dari Website World In Data yang mengambil langsung dari FAO ( Food Agriculture Organization).

Import Library

library(dLagM)

## Warning: package 'dLagM' was built under R version 4.3.3

## Loading required package: nardl

## Warning: package 'nardl' was built under R version 4.3.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.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.2

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(dynlm)
library(MLmetrics)

## Warning: package 'MLmetrics' was built under R version 4.3.3

## 
## Attaching package: 'MLmetrics'

## The following object is masked from 'package:dLagM':
## 
##     MAPE

## The following object is masked from 'package:base':
## 
##     Recall

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.3.3

library(car)

## Warning: package 'car' was built under R version 4.3.2

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.3.2

library(readxl)

## Warning: package 'readxl' was built under R version 4.3.2

Import Data

dataikan <- read_excel("C:/Users/Admin/Downloads/Adriano Excel Putra/data mpdw ikan.xlsx", sheet = 1)

head(dataikan)

## # A tibble: 6 × 4
##   Years Consumption `Yt-1` Production
##   <dbl>       <dbl>  <dbl>      <dbl>
## 1  1961       10.3   NA        906800
## 2  1962       10.1   10.3      943200
## 3  1963        9.68  10.1      936200
## 4  1964       10.1    9.68    1000000
## 5  1965       10.5   10.1     1066800
## 6  1966       11.5   10.5     1201600

Visualisasi

model <- lm (Consumption ~ Production, data = dataikan)
summary(model)

## 
## Call:
## lm(formula = Consumption ~ Production, data = dataikan)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.73904 -0.46515  0.04785  0.43478  1.56746 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.534e+00  1.500e-01   43.58   <2e-16 ***
## Production  2.819e-06  2.543e-08  110.84   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7418 on 59 degrees of freedom
## Multiple R-squared:  0.9952, Adjusted R-squared:  0.9951 
## F-statistic: 1.229e+04 on 1 and 59 DF,  p-value: < 2.2e-16

plotproduksi <- plot(dataikan$Years, dataikan$Production)

plotkonsumsi<- plot(dataikan$Years, dataikan$Consumption)

Dapat dilihat secara visualisasi bahwa terdapat korelasi yang sangat kuat antara Produksi dan Konsumsi Ikan dilihat dari persebaran datanya yang cenderung Mirip. Untuk itu data ini akan dilanjutkan ke pemodelan. Untuk melihat apakah ada pengaruh dari data-data ditahun tahun sebelumnya terhadap nilai produksi dan utamanya konsumsi ikan.

Splitting Data

Data akan di split menggunakan pembagian 80 : 20. 80 untuk Data Training dan 20 Untuk Data Testing. Sehingga saya memakai 48 data untuk training dan 13 data untuk testing.

#SPLIT DATA
train<-dataikan[1:48,]
test<-dataikan[49:61,]
head(train)

## # A tibble: 6 × 4
##   Years Consumption `Yt-1` Production
##   <dbl>       <dbl>  <dbl>      <dbl>
## 1  1961       10.3   NA        906800
## 2  1962       10.1   10.3      943200
## 3  1963        9.68  10.1      936200
## 4  1964       10.1    9.68    1000000
## 5  1965       10.5   10.1     1066800
## 6  1966       11.5   10.5     1201600

Converting data ke Time Series

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

##      Years Consumption      Yt-1 Production
## [1,]  1961   10.284873        NA     906800
## [2,]  1962   10.143362 10.284873     943200
## [3,]  1963    9.675091 10.143362     936200
## [4,]  1964   10.062971  9.675091    1000000
## [5,]  1965   10.453741 10.062971    1066800
## [6,]  1966   11.488996 10.453741    1201600

Model Koyck

Model Koyck didasarkan pada asumsi bahwa semakin jauh jarak lag peubah independen dari periode sekarang maka semakin kecil pengaruh peubah lag terhadap peubah dependen. Koyck mengusulkan suatu metode untuk menduga model dinamis distributed lag dengan mengasumsikan bahwa semua koefisien \(\beta\) mempunyai tanda sama.

Pemodelan

model.koyck <- koyckDlm(x= train$Production, y=train$Consumption)
summary(model.koyck)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73447 -0.25741  0.03308  0.23165  1.00246 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.145e+00  8.430e-01   2.545  0.01452 *  
## Y.1         7.020e-01  1.134e-01   6.190 1.78e-07 ***
## X.t         8.342e-07  2.722e-07   3.064  0.00372 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3653 on 44 degrees of freedom
## Multiple R-Squared: 0.9935,  Adjusted R-squared: 0.9933 
## Wald test:  3384 on 2 and 44 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha         beta       phi
## Geometric coefficients:  7.199906 8.341517e-07 0.7020385

AIC(model.koyck)

## [1] 43.61587

BIC(model.koyck)

## [1] 51.01646

Dari hasil tersebut, didapat bahwa peubah \(x_t\) dan \(y_{t-1}\) memiliki nilai \(P-Value<0.05\). Hal ini menunjukkan bahwa peubah \(x_t\) dan \(y_{t-1}\) berpengaruh signifikan terhadap \(y\). Selain itu dapat dilihat bahwa interceptnya juga berpengaruh secara signfikan terhadap model . Sehingga Model keseluruhannya adalah sebagai berikut

\[ \hat{Y_t}=2145+8.342*10^-7X_t+0.7020Y_{t-1} \]

Selain dari Model, didapat pula AIC model ini cukup bagus di Angka 43.61587 yang nantinya akan dibandingkan dengan model-model lainnya. Namun Model Koyck sudah mengset standar yang cukup untuk ukuran kebaikan modelnya.

Forecasting

fore.koyck <- forecast(model = model.koyck, x=test$Production, h=13)
fore.koyck

## $forecasts
##  [1] 25.22736 26.27145 27.65535 29.00482 30.93192 32.79665 34.37230 35.82103
##  [9] 37.50843 39.02413 40.74384 41.95113 42.79869
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Production, h = 13)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

mape.koyck <- MAPE(fore.koyck$forecasts, test$Consumption)
mape.koyck

## [1] 0.08192264

#akurasi data training
GoF(model.koyck)

##              n       MAE           MPE      MAPE      sMAPE      MASE       MSE
## model.koyck 47 0.2780909 -0.0007959849 0.0198101 0.01982738 0.6114765 0.1249186
##                 MRAE     GMRAE
## model.koyck 92.06818 0.6179788

Didapatkan MAPE Forecasting untuk 13 data kedepan adalah 0.0819. Nilai MAPE yang kecil ini mengindikasikan bahwa forecasting/ peramalan data sangat bagus dalam hal akursainya. Tentunya Model ini belum dapat dikatakan model yang sempurna dan perlu di lihat asumsi-asumsi yang ada mengenai model ini.

Dapat dilihat pula pada Gof( Model.Koyck) Untuk akurasi data traininginya, bernilai di angka 0.0198, tentunya MAPE ini mengindikasikan model menghasilkan output yang sudah sangat baik.

Regression With Distribution LAG

Pada regresi ini, Model akan diikuti dengan LAG. Lag sendiri adalah Waktu yang diperlukan oleh X dalam mempengengaruhi Y. Artinya dalam konteks data yang saya gunakan, waktu / periode yang dibutuhkan bagi data produksi ikan dalam mempengaruhi jumlah konsumsi ikan.

Pada tahap awal saya akan mencoba memakai lag 2, yaitu produksi ikan 2 periode yang lalu mempengaruhi konsumsi ikan.

Lag= 2 (Model)

model.dlm2 <- dlm(x = train$Production,y = train$Consumption , q = 2)
summary(model.dlm2)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.77295 -0.35433 -0.00729  0.26913  1.08036 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.270e+00  1.348e-01  53.919  < 2e-16 ***
## x.t          3.388e-06  1.097e-06   3.088  0.00357 ** 
## x.1          4.632e-07  1.555e-06   0.298  0.76729    
## x.2         -1.427e-06  1.142e-06  -1.250  0.21837    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.461 on 42 degrees of freedom
## Multiple R-squared:   0.99,  Adjusted R-squared:  0.9892 
## F-statistic:  1381 on 3 and 42 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 65.10808 74.25128

AIC(model.dlm2)

## [1] 65.10808

BIC(model.dlm2)

## [1] 74.25128

Dapat dilihat jikalau p-value dari Intercept dan X.t < alpha 5%. Hal ini mengindikasikan bahwa intercept dan Xt berpengaruh secara signifikan terhadap model. Model keseleruhan yang terbentuk mendapat tambahan 2 lag yang tadi dimasukkan ke dalam model sehingga modelnya menjadi :

\[ \hat{Y_t} = 7270 + 3.388 \times 10^{-6} X_t + 4.632 \times 10^{-7} X_{t-1} - 1.427 \times 10^{-6} X_{t-2} \]

Model ini memiliki nilai AIC yang cukup tinggi di angka 65. Dimana model ini bisa saja dibandingkan dengan lag - lag yang berbeda untuk mendapat Nilai AIC yang lebih rendah. Namun jika dibandingkan dengan Model Koyck sebelumnya tentu Model dengan DLM ( Distributed Lag Model) dengan Lag 2 ini masih kurang efektif , mengingat prinsip parsimoni yang digunakaan dalam pemilihan model.

Lag Optimum (Model)

Ketika melakukan pencarian Lag Optimum by Trial and Error, Data yang dipakai ini mampu menoleransi hingga lag 22, dimana akan Mendapatkan AIC paling rendah, Namun jika dilihat dalam pemodelan , tentu saja model dengan Lag 22 sangat tidak relevan karena tidak terdapat peubah yang signfikan terhadap model :

Percobaan Lag Maksimum (22)

model.dlm22 <- dlm(x = train$Production,y = train$Consumption , q = 22)
summary(model.dlm22)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##         1         2         3         4         5         6         7         8 
## -0.179263  0.001692  0.085419  0.168092  0.057261 -0.041116 -0.157061 -0.101235 
##         9        10        11        12        13        14        15        16 
## -0.008131  0.092005  0.097520  0.120058  0.003877 -0.029751 -0.067667 -0.023993 
##        17        18        19        20        21        22        23        24 
## -0.091583  0.076012  0.046104  0.056418 -0.134952 -0.075604  0.065531  0.047915 
##        25        26 
##  0.016699 -0.024247 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  3.888e+00  1.503e+00   2.587    0.123
## x.t          3.902e-06  1.581e-06   2.467    0.132
## x.1          5.349e-07  1.616e-06   0.331    0.772
## x.2         -3.528e-08  1.933e-06  -0.018    0.987
## x.3         -6.452e-07  1.846e-06  -0.349    0.760
## x.4         -7.180e-07  1.971e-06  -0.364    0.751
## x.5          5.765e-07  2.146e-06   0.269    0.813
## x.6          2.683e-06  2.335e-06   1.149    0.369
## x.7         -1.530e-06  2.338e-06  -0.654    0.580
## x.8         -1.391e-06  2.484e-06  -0.560    0.632
## x.9         -3.363e-07  2.567e-06  -0.131    0.908
## x.10        -3.260e-06  2.489e-06  -1.310    0.320
## x.11        -2.081e-06  2.636e-06  -0.790    0.513
## x.12         4.722e-07  2.754e-06   0.171    0.880
## x.13         1.123e-06  2.602e-06   0.432    0.708
## x.14        -1.372e-06  2.727e-06  -0.503    0.665
## x.15        -2.798e-06  2.894e-06  -0.967    0.436
## x.16         2.105e-06  2.871e-06   0.733    0.540
## x.17         2.628e-06  2.511e-06   1.046    0.405
## x.18         2.860e-06  2.263e-06   1.264    0.334
## x.19        -2.361e-07  2.745e-06  -0.086    0.939
## x.20        -1.041e-07  2.318e-06  -0.045    0.968
## x.21         1.744e-06  2.476e-06   0.705    0.554
## x.22         1.531e-06  2.912e-06   0.526    0.652
## 
## Residual standard error: 0.3138 on 2 degrees of freedom
## Multiple R-squared:  0.9993, Adjusted R-squared:  0.9917 
## F-statistic: 130.9 on 23 and 2 DF,  p-value: 0.007605
## 
## AIC and BIC values for the model:
##         AIC    BIC
## 1 -3.171416 28.281

AIC(model.dlm22)

## [1] -3.171416

BIC(model.dlm22)

## [1] 28.281

Model ini tentu menjadi Model dengan AIC paling kecil, Namun tentunya relevansi jjuga perlu diperhatikan dalam Model. Dapat dilihat walaupun secara Uji statsitik p-valuenya dibawah 0.05 namun secara partial tidak ada satupun koefisien yang berpengaruh terhadap model. Tentunya model ini tidak dapat dikatakan model yang paling baik . Maka dari itu dalam pencarian Lag Optimum, saya menurunkan LAG maximummnya menjadi 15 , karena di lag 15 ke atas sudah tidak adalagi koefisien yang berpengaruh terhadap model.

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

##    q - k    MASE      AIC      BIC   GMRAE   MBRAE R.Adj.Sq   Ljung-Box
## 13    13 0.39579 30.52319 55.40876 0.45403 0.34481   0.9945 0.001323713

Dari informasi ini, didapat Lag optimum di 13 sehingga saya akan memodelkan Lag 13

model.dlmopt <- dlm(x = train$Production,y = train$Consumption , q = 13)
summary(model.dlmopt)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34814 -0.20061 -0.02353  0.14307  0.52561 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.990e+00  1.516e-01  46.114   <2e-16 ***
## x.t          2.664e-06  9.699e-07   2.746   0.0124 *  
## x.1          4.264e-07  1.358e-06   0.314   0.7568    
## x.2         -3.983e-07  1.402e-06  -0.284   0.7792    
## x.3         -1.108e-06  1.467e-06  -0.755   0.4591    
## x.4          2.553e-07  1.504e-06   0.170   0.8669    
## x.5          1.058e-06  1.527e-06   0.693   0.4964    
## x.6          2.471e-06  1.811e-06   1.364   0.1876    
## x.7          4.608e-08  2.011e-06   0.023   0.9819    
## x.8         -2.292e-06  1.858e-06  -1.233   0.2319    
## x.9         -1.044e-06  1.879e-06  -0.556   0.5845    
## x.10        -2.138e-06  1.908e-06  -1.121   0.2758    
## x.11         1.558e-07  1.903e-06   0.082   0.9356    
## x.12         8.787e-07  1.886e-06   0.466   0.6462    
## x.13         1.683e-06  1.352e-06   1.245   0.2277    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3134 on 20 degrees of freedom
## Multiple R-squared:  0.9968, Adjusted R-squared:  0.9945 
## F-statistic: 440.2 on 14 and 20 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 30.52319 55.40876

AIC(model.dlmopt)

## [1] 30.52319

BIC(model.dlmopt)

## [1] 55.40876

Sama seperti Lag 2 , hanya dua koefisien yang berpengaruh nyata pada taraf 5 persen terhadap model yaitu intercept dan Xt, namun perlu diperhatikan bahwa Nilai AIC yang didapat jauh lebih kecil daripada LAG 2 , sehingga model ini secara simultan lebih baik daripada model dengan LAG 2, walaupun secara substansial kompleksitas model ini jauh lebih tinggi akibat adanya 13 waktu yang mempengaruhi peubah respon.

Melihat dari AIC, bahwa AIC model ini lebih rendah dari perbandingan 2 model sebelumnya yaitu model Koyck dan Model dengan distribsui lag 2 , sehingga untuk saat ini, model inilah yang masih menjadi acuan utama dalam hubungan produksi dan konsumsi ikan.

Forecasting

Forecasting DLM 2

fore.dlm2 <- forecast(model = model.dlm2, x=test$Production, h=13)
fore.dlm2

## $forecasts
##  [1] 24.34732 26.95435 29.78536 30.45747 33.53335 35.50991 35.20069 35.86460
##  [9] 38.32170 39.45270 41.15291 40.95053 39.82890
## 
## $call
## forecast.dlm(model = model.dlm2, x = test$Production, h = 13)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

mape.dlm2 <- MAPE(fore.dlm2$forecasts, test$Consumption)
mape.dlm2

## [1] 0.06431693

#akurasi data training
GoF(model.dlm2)

##             n       MAE          MPE       MAPE      sMAPE      MASE       MSE
## model.dlm2 46 0.3612105 -0.000960807 0.02664929 0.02665931 0.7947666 0.1940042
##                MRAE    GMRAE
## model.dlm2 25.32116 0.908793

Didapatkan MAPE yang bernilai 0.064 dalam data testing . Hal ini menunjukkan bahwa prediksi dengan Distributed Lag Model sangat bagus dengan tingkat kesalahan yang tinggi dimana dibawah nilai 10 %. Tentu saja tingkat prediksi dengan akurasi kesalahan yang mendekati 0 sangat baik dalam konteks forecasting data. Namun hal ini perlu dilihat juga dan dibandingkan dengan model dengan lag yang berbeda

Forecasting DLM 13

fore.dlmopt <- forecast(model = model.dlmopt, x=test$Production, h=13)
mape.dlmopt<- MAPE(fore.dlmopt$forecasts, test$Consumption)
mape.dlmopt

## [1] 0.04150938

#akurasi data training
GoF(model.dlmopt)

##               n       MAE           MPE       MAPE      sMAPE      MASE
## model.dlmopt 35 0.1899821 -0.0003420733 0.01296415 0.01298084 0.3957873
##                     MSE     MRAE     GMRAE
## model.dlmopt 0.05613074 75.47902 0.4540272

Ketika kita menggunakan Model DLM dengan Lag optimum (13). Didapatkan hasil MAPE yang lebih rendah daripada dengan Lag 2. Hal ini mengindikasikan , walaupun sama-sama bagus hasil MAPEnya, yaitu mape yang sangat rendah, namun perlu diperhatikan bahwa Nilai Mape dengan lag optimum lebih rendah daripada Mape dengan Lag 2 , Sehingga bisa dikatakan model dengan lag 13 lebih baik daripada model dengan Lag 2 berdasar pada nilai MAPEnya.

Autoregressive Model

Autoregressive (AR) model adalah salah satu model dalam analisis deret waktu (time series) yang digunakan untuk memprediksi nilai masa depan dari sebuah variabel berdasarkan nilai-nilai masa lalunya. Artinya dalam model ini, si peubah dependen di waktu lalu akan dikut seratakan ke dalam model .

Peubah dependen dipengaruhi oleh peubah independen pada waktu sekarang, serta dipengaruhi juga oleh peubah dependen itu sendiri pada satu waktu yang lalu maka model tersebut disebut autoregressive (Gujarati 2004).

Pada Autrorgressive, p disini sebagai panjang lag dari variabel dependen dan q sebagai ordo dari model autorgegressive yaitu panjang lag dari variabel independen.

Percobaan dengan P=2, Q=2

Pada Simulasi kali ini, saya mencoba membangkitkan Lag pada Peubah Y sebesar 2 dan Lag pada Peubah X sebesar 2. Artinya dalam model nanti akan ditambahkan dua data fish consumption dan dua data fish production dari periode sebelumnya. Maka dari itu saya ingin melihat bagaimana model autoregressive ini .

model.ardl22 <- ardlDlm(formula = Consumption ~ Production, 
                         data = train,p = 2, q = 2)
summary(model.ardl22)

## 
## Time series regression with "ts" data:
## Start = 3, End = 48
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72223 -0.14779  0.00238  0.14091  0.67025 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.576e+00  7.623e-01   3.379  0.00163 ** 
## Production.t   3.846e-06  6.877e-07   5.593 1.76e-06 ***
## Production.1  -3.996e-06  1.126e-06  -3.547  0.00101 ** 
## Production.2   1.000e-06  9.270e-07   1.079  0.28701    
## Consumption.1  1.071e+00  1.441e-01   7.432 4.73e-09 ***
## Consumption.2 -4.288e-01  1.433e-01  -2.993  0.00472 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2848 on 40 degrees of freedom
## Multiple R-squared:  0.9964, Adjusted R-squared:  0.9959 
## F-statistic:  2185 on 5 and 40 DF,  p-value: < 2.2e-16

AIC(model.ardl22)

## [1] 22.56471

BIC(model.ardl22)

## [1] 35.3652

Pada Model ARDL, Didapat hal yang menarik. Dimulai dari Modelnya, Hampir Semua Peubah berpengaruh secara signifikan terhadap Model. Dilihat dari p-valuenya yang 5 peubahnya dibawah dari alpha (5 %). Sedangkan untuk Production T-2 , tidak signfikan pada taraf 5 persen. Model tersebut dapat ditulis sebagai berikut :

\[ \hat{Y_t} = 2567 + 3.846 \times 10^{-6} X_t -3.996 \times 10^{-6} X_{t-1} +1 \times 10^{-6} x_{t-2} + 1071 Y_{t-1}- 4.228 \times 10^{-1} Y_{t-2} \]

Model Ini juga menghasilkan AIC yang jauh lebih rendah daripada perbandingan model-model yang dipaka sebelumnya ( DLM dan Koyck). sehingga dapat dikatakan sampai sekarang model ini adalah model yang paling baik dalam kasus regresi ini.

fore.ardl22 <- forecast(model = model.ardl22, x=test$Production, h=13)
fore.ardl22

## $forecasts
##  [1] 24.91718 27.61834 30.16495 31.21759 34.74109 36.18453 36.17199 36.45630
##  [9] 38.53577 39.36799 41.60887 40.90897 39.98454
## 
## $call
## forecast.ardlDlm(model = model.ardl22, x = test$Production, h = 13)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

mape.ardl22 <- MAPE(fore.ardl22$forecasts, test$Consumption)
mape.ardl22

## [1] 0.05088317

#akurasi data training
GoF(model.ardl22)

##               n       MAE           MPE      MAPE      sMAPE      MASE
## model.ardl22 46 0.2031951 -0.0002589521 0.0142199 0.01424849 0.4470874
##                     MSE    MRAE     GMRAE
## model.ardl22 0.07053207 49.5787 0.4850728

Berdasarkan akurasi di atas, terlihat bahwa nilai MAPE keduanya tidak jauh berbeda. Artinya, model regresi dengan distribusi lag ini tidak overfitted atau underfitted. Hal ini dapat dilihat pada MAPE testing dan Traininginyna yaitu di angka 0.05 dan 0,0140.

Disini batas q yang saya masukkan disamakan dengan DLM , yaitu q dan p optimum 15. Memanfaatkan Rumus Model untuk mendapatkan AIC, saya menerapkan looping sehingga didapatkan di akhir nanti kombinasi p dan q untuk mendapatkan AIC paling minimumnya

# Inisialisasi untuk menyimpan AIC minimal untuk setiap p dan q
min_aic_values <- matrix(NA, nrow = 15, ncol = 15)

# Loop untuk mencari AIC minimum dari kombinasi p dan q
for (q in 1:15) {
  for (p in 1:15) {
    # Buat model ARDL dengan parameter p dan q
    model_ardl <- ardlDlm(formula = Consumption ~ Production, 
                          data = train, p = p, q = q)
    
    # Hitung AIC untuk model
    aic_value <- AIC(model_ardl)
    
    # Simpan AIC dalam matriks
    min_aic_values[p, q] <- aic_value
  }
}

## [1] 28.38554
## [1] 29.85802
## [1] 30.7276
## [1] 31.43028
## [1] 32.40821
## [1] 25.54818
## [1] 25.81485
## [1] 28.33532
## [1] 27.99214
## [1] 28.23066
## [1] 20.82553
## [1] 17.45405
## [1] 19.15036
## [1] 22.45594
## [1] 24.50377
## [1] 21.8847
## [1] 22.56471
## [1] 23.75827
## [1] 23.8744
## [1] 23.67741
## [1] 15.91812
## [1] 16.17991
## [1] 17.63958
## [1] 13.72842
## [1] 16.17615
## [1] 13.52081
## [1] 13.74742
## [1] 15.98605
## [1] 19.45548
## [1] 17.29725
## [1] 21.13585
## [1] 21.49325
## [1] 23.48089
## [1] 24.10298
## [1] 24.83272
## [1] 17.30889
## [1] 17.33925
## [1] 17.96003
## [1] 15.39217
## [1] 16.9392
## [1] 12.25176
## [1] 10.52308
## [1] 11.39594
## [1] 15.04038
## [1] 14.49166
## [1] 21.75867
## [1] 22.05457
## [1] 24.04136
## [1] 25.80178
## [1] 26.62444
## [1] 15.86592
## [1] 14.12483
## [1] 12.81397
## [1] 12.67194
## [1] 15.22285
## [1] 12.03583
## [1] 10.84092
## [1] 12.63868
## [1] 16.35858
## [1] 15.73591
## [1] 24.15236
## [1] 24.33165
## [1] 26.00641
## [1] 27.84308
## [1] 28.62384
## [1] 17.7898
## [1] 16.0305
## [1] 14.80914
## [1] 14.02176
## [1] 16.19889
## [1] 11.12374
## [1] 11.85584
## [1] 13.87954
## [1] 17.62548
## [1] 17.70211
## [1] 17.37822
## [1] 15.14052
## [1] 16.91407
## [1] 17.41903
## [1] 17.41085
## [1] 19.09526
## [1] 17.86222
## [1] 13.86646
## [1] 13.70903
## [1] 16.64166
## [1] 10.78977
## [1] 12.68438
## [1] 15.07491
## [1] 18.72787
## [1] 19.69992
## [1] 19.75143
## [1] 17.58761
## [1] 19.48748
## [1] 19.6699
## [1] 20.25213
## [1] 21.82027
## [1] 16.04521
## [1] 12.9219
## [1] 14.66306
## [1] 16.70266
## [1] 9.35225
## [1] 11.98367
## [1] 13.04459
## [1] 16.80707
## [1] 17.16053
## [1] 17.95089
## [1] 13.70057
## [1] 15.60763
## [1] 12.22769
## [1] 12.35162
## [1] 14.33027
## [1] 7.220466
## [1] 4.360477
## [1] 4.940492
## [1] 5.698817
## [1] 6.982547
## [1] 10.58925
## [1] 9.812525
## [1] 9.320213
## [1] 10.83223
## [1] 13.73471
## [1] 11.03998
## [1] 12.23791
## [1] 7.851494
## [1] 9.732893
## [1] 11.72342
## [1] 6.823763
## [1] 4.354272
## [1] 3.348259
## [1] 4.479413
## [1] 6.430515
## [1] 9.880394
## [1] 8.115309
## [1] 6.418357
## [1] 5.255225
## [1] 16.74911
## [1] 14.09711
## [1] 15.49805
## [1] 9.715686
## [1] 11.68428
## [1] 13.67101
## [1] 9.162099
## [1] 6.279653
## [1] 5.47703
## [1] 4.545278
## [1] 7.469476
## [1] 9.797002
## [1] 7.458291
## [1] 6.459869
## [1] 6.775005
## [1] 17.19073
## [1] 14.25368
## [1] 14.80513
## [1] 11.10635
## [1] 13.07486
## [1] 14.82212
## [1] 10.90294
## [1] 8.564346
## [1] 6.09762
## [1] 3.685906
## [1] -8.626381
## [1] -6.697092
## [1] -2.580391
## [1] -9.343942
## [1] -9.271509
## [1] 18.04193
## [1] 16.18175
## [1] 17.68767
## [1] 14.10556
## [1] 16.03116
## [1] 17.7042
## [1] 14.14466
## [1] 11.3198
## [1] 7.952939
## [1] 4.828465
## [1] -5.317768
## [1] -5.552398
## [1] -1.708309
## [1] -10.30141
## [1] -13.09893
## [1] 20.20872
## [1] 17.61204
## [1] 19.1797
## [1] 15.6892
## [1] 17.27481
## [1] 18.95837
## [1] 15.93325
## [1] 11.69353
## [1] 7.774066
## [1] 4.861895
## [1] -3.419191
## [1] -2.457671
## [1] -5.75848
## [1] -11.0271
## [1] -11.11154
## [1] 23.4267
## [1] 20.98945
## [1] 22.41826
## [1] 17.55415
## [1] 17.97526
## [1] 18.93044
## [1] 17.81284
## [1] 10.08368
## [1] 7.085752
## [1] -9.585251
## [1] -18.58608
## [1] -19.0717
## [1] -17.07853
## [1] -22.97424
## [1] -21.03103
## [1] 26.38234
## [1] 24.25993
## [1] 25.59715
## [1] 19.87161
## [1] 19.70768
## [1] 20.69731
## [1] 20.35197
## [1] 8.460101
## [1] 3.025509
## [1] -22.44727
## [1] -20.48175
## [1] -21.52089
## [1] -23.50876
## [1] -21.69054
## [1] -19.70949

# Cari kombinasi p dan q yang memberikan AIC minimum
min_aic <- min(min_aic_values, na.rm = TRUE)
optimum_indices <- which(min_aic_values == min_aic, arr.ind = TRUE)

# Ambil p dan q optimum dari indeks matriks
p_opt <- optimum_indices[1, "row"]
q_opt <- optimum_indices[1, "col"]

# Hasil dalam bentuk data frame
result <- data.frame("q_optimum" = q_opt, "p_optimum" = p_opt, "AIC" = min_aic)

print(result)

##     q_optimum p_optimum       AIC
## col        15        13 -23.50876

Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=13\) dan \(q=15\), yaitu sebesar \[-23.50786\] Artinya, model autoregressive optimum didapat ketika \(p=13\) dan \(q=15\).

Selanjutnya dapat dilakukan pemodelan dengan nilai \(p\) dan \(q\) optimum seperti inisialisasi di langkah sebelumnya.

Percobaan P dan Q optimum

model.ardlopt <- ardlDlm(formula = Consumption ~ Production, 
                         data = train,p = 13, q = 15)
summary(model.ardlopt)

## 
## Time series regression with "ts" data:
## Start = 16, End = 48
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##         16         17         18         19         20         21         22 
##  0.0356177 -0.0324364 -0.0569201  0.0478216 -0.0140189  0.0601467  0.0126431 
##         23         24         25         26         27         28         29 
## -0.0585870 -0.0034234 -0.0795482  0.0030427  0.1227632  0.0119726 -0.0483270 
##         30         31         32         33         34         35         36 
##  0.0538041 -0.0170377 -0.0428077 -0.0672511  0.1092417 -0.0385621 -0.0278510 
##         37         38         39         40         41         42         43 
##  0.0969965 -0.0223928 -0.1524816  0.0826292  0.1303333 -0.1349213 -0.0003122 
##         44         45         46         47         48 
##  0.0590288 -0.0405853 -0.0157929  0.0355176 -0.0083020 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## (Intercept)     2.960e+01  1.518e+01   1.950    0.146
## Production.t    3.083e-07  4.236e-06   0.073    0.947
## Production.1   -3.499e-06  2.328e-06  -1.503    0.230
## Production.2    1.020e-05  8.034e-06   1.269    0.294
## Production.3    2.236e-06  5.723e-06   0.391    0.722
## Production.4   -1.008e-05  7.385e-06  -1.364    0.266
## Production.5    1.178e-05  9.451e-06   1.246    0.301
## Production.6    4.232e-06  6.037e-06   0.701    0.534
## Production.7   -2.112e-05  1.638e-05  -1.289    0.288
## Production.8    7.961e-06  9.836e-06   0.809    0.478
## Production.9    1.872e-05  1.341e-05   1.396    0.257
## Production.10  -9.476e-06  6.120e-06  -1.549    0.219
## Production.11  -1.367e-07  8.337e-06  -0.016    0.988
## Production.12   6.582e-06  7.802e-06   0.844    0.461
## Production.13  -9.126e-06  1.470e-05  -0.621    0.579
## Consumption.1   1.081e+00  5.781e-01   1.870    0.158
## Consumption.2  -1.515e+00  9.681e-01  -1.565    0.216
## Consumption.3  -5.349e-01  1.074e+00  -0.498    0.653
## Consumption.4   7.404e-01  8.022e-01   0.923    0.424
## Consumption.5  -2.039e+00  1.716e+00  -1.189    0.320
## Consumption.6   2.520e-01  1.097e+00   0.230    0.833
## Consumption.7   2.779e+00  2.417e+00   1.150    0.334
## Consumption.8  -2.214e+00  1.776e+00  -1.247    0.301
## Consumption.9  -1.733e+00  1.544e+00  -1.123    0.343
## Consumption.10  1.113e+00  7.610e-01   1.462    0.240
## Consumption.11 -2.999e-01  1.014e+00  -0.296    0.787
## Consumption.12 -1.067e+00  1.152e+00  -0.926    0.423
## Consumption.13  1.261e+00  1.879e+00   0.671    0.550
## Consumption.14 -4.932e-01  4.604e-01  -1.071    0.363
## Consumption.15 -3.478e-01  3.141e-01  -1.108    0.349
## 
## Residual standard error: 0.2197 on 3 degrees of freedom
## Multiple R-squared:  0.9997, Adjusted R-squared:  0.9971 
## F-statistic:   377 on 29 and 3 DF,  p-value: 0.0001931

AIC(model.ardlopt)

## [1] -23.50876

BIC(model.ardlopt)

## [1] 22.88297

Output diatas adalah model yang memakai p = 13 dan q = 15, sesuai dengan fungsi loop yang dicari dalam mencari p dan q optimumnya. Memang AIC menunjukkan angka yang sangat lebih kecil dibanding dengan fungsi sebelumnya di p=1 dan q = 2 namun ada beberapa pertimbangan seperti dalam model ini, Tidak terdapat peubah yang berpengaruh secara signfikan terhadap model. Sehingga walaupun Model dapat dikatakan lebih baik dari cara pandang AIC, namun ketika kita memandang model secara utuh, model sebelumnya jauh lebih baik

Sehingga dalam hal ini , saya memilih model ARDL dengan p= 2 dan q= 2 sebagai pembanding dengan model-model sebelumnya.

Pemodelan dengan Packages

Untuk menentukan formula model yang akan digunakan, tersedia fungsi tambahan yang memungkinkan spesifikasi dinamika (melalui d() dan L()) atau pola linier/siklus dengan mudah (melalui trend(), season(), dan harmon()). Semua fungsi formula baru mengharuskan argumennya berupa objek deret waktu (yaitu, "ts" atau "zoo").

# sama dengan model dlm q = 2
cons_lm1 <- dynlm(Consumption ~ Production + L(Production,1) + L(Production,2), data = train.ts)

# sama dengan model dlm q = 13 (dengan lag 1 sampai 13)
cons_lm2 <- dynlm(Consumption ~ Production + L(Production,1) + L(Production,2) + L(Production,3) + 
                  L(Production,4) + L(Production,5) + L(Production,6) + L(Production,7) + 
                  L(Production,8) + L(Production,9) + L(Production,10) + L(Production,11) + 
                  L(Production,12) + L(Production,13), data = train.ts)

# sama dengan ARDL p = 2, q = 2
cons_lm3 <- dynlm(Consumption ~ Production + L(Consumption,1)+L(Consumption,2) + L(Production,1)+L(Production,2), data = train.ts)



# Model dengan lag Production hingga 13 dan lag Consumption hingga 15
cons_lm4 <- dynlm(
  Consumption ~ Production + 
                L(Production, 1) + L(Production, 2) + L(Production, 3) + 
                L(Production, 4) + L(Production, 5) + L(Production, 6) + 
                L(Production, 7) + L(Production, 8) + L(Production, 9) + 
                L(Production, 10) + L(Production, 11) + L(Production, 12) + 
                L(Production, 13)  +
                L(Consumption, 1) + L(Consumption, 2) + L(Consumption, 3) + 
                L(Consumption, 4) + L(Consumption, 5) + L(Consumption, 6) + 
                L(Consumption, 7) + L(Consumption, 8) + L(Consumption, 9) + 
                L(Consumption, 10) + L(Consumption, 11) + L(Consumption, 12) + 
                L(Consumption, 13) + L(Consumption, 14) + L(Consumption, 15),
  data = train.ts
)

Summary Model

summary(cons_lm1)

## 
## Time series regression with "ts" data:
## Start = 3, End = 48
## 
## Call:
## dynlm(formula = Consumption ~ Production + L(Production, 1) + 
##     L(Production, 2), data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.77295 -0.35433 -0.00729  0.26913  1.08036 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       7.270e+00  1.348e-01  53.919  < 2e-16 ***
## Production        3.388e-06  1.097e-06   3.088  0.00357 ** 
## L(Production, 1)  4.632e-07  1.555e-06   0.298  0.76729    
## L(Production, 2) -1.427e-06  1.142e-06  -1.250  0.21837    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.461 on 42 degrees of freedom
## Multiple R-squared:   0.99,  Adjusted R-squared:  0.9892 
## F-statistic:  1381 on 3 and 42 DF,  p-value: < 2.2e-16

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 14, End = 48
## 
## Call:
## dynlm(formula = Consumption ~ Production + L(Production, 1) + 
##     L(Production, 2) + L(Production, 3) + L(Production, 4) + 
##     L(Production, 5) + L(Production, 6) + L(Production, 7) + 
##     L(Production, 8) + L(Production, 9) + L(Production, 10) + 
##     L(Production, 11) + L(Production, 12) + L(Production, 13), 
##     data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34814 -0.20061 -0.02353  0.14307  0.52561 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        6.990e+00  1.516e-01  46.114   <2e-16 ***
## Production         2.664e-06  9.699e-07   2.746   0.0124 *  
## L(Production, 1)   4.264e-07  1.358e-06   0.314   0.7568    
## L(Production, 2)  -3.983e-07  1.402e-06  -0.284   0.7792    
## L(Production, 3)  -1.108e-06  1.467e-06  -0.755   0.4591    
## L(Production, 4)   2.553e-07  1.504e-06   0.170   0.8669    
## L(Production, 5)   1.058e-06  1.527e-06   0.693   0.4964    
## L(Production, 6)   2.471e-06  1.811e-06   1.364   0.1876    
## L(Production, 7)   4.608e-08  2.011e-06   0.023   0.9819    
## L(Production, 8)  -2.292e-06  1.858e-06  -1.233   0.2319    
## L(Production, 9)  -1.044e-06  1.879e-06  -0.556   0.5845    
## L(Production, 10) -2.138e-06  1.908e-06  -1.121   0.2758    
## L(Production, 11)  1.558e-07  1.903e-06   0.082   0.9356    
## L(Production, 12)  8.787e-07  1.886e-06   0.466   0.6462    
## L(Production, 13)  1.683e-06  1.352e-06   1.245   0.2277    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3134 on 20 degrees of freedom
## Multiple R-squared:  0.9968, Adjusted R-squared:  0.9945 
## F-statistic: 440.2 on 14 and 20 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 3, End = 48
## 
## Call:
## dynlm(formula = Consumption ~ Production + L(Consumption, 1) + 
##     L(Consumption, 2) + L(Production, 1) + L(Production, 2), 
##     data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72223 -0.14779  0.00238  0.14091  0.67025 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.576e+00  7.623e-01   3.379  0.00163 ** 
## Production         3.846e-06  6.877e-07   5.593 1.76e-06 ***
## L(Consumption, 1)  1.071e+00  1.441e-01   7.432 4.73e-09 ***
## L(Consumption, 2) -4.288e-01  1.433e-01  -2.993  0.00472 ** 
## L(Production, 1)  -3.996e-06  1.126e-06  -3.547  0.00101 ** 
## L(Production, 2)   1.000e-06  9.270e-07   1.079  0.28701    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2848 on 40 degrees of freedom
## Multiple R-squared:  0.9964, Adjusted R-squared:  0.9959 
## F-statistic:  2185 on 5 and 40 DF,  p-value: < 2.2e-16

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 16, End = 48
## 
## Call:
## dynlm(formula = Consumption ~ Production + L(Production, 1) + 
##     L(Production, 2) + L(Production, 3) + L(Production, 4) + 
##     L(Production, 5) + L(Production, 6) + L(Production, 7) + 
##     L(Production, 8) + L(Production, 9) + L(Production, 10) + 
##     L(Production, 11) + L(Production, 12) + L(Production, 13) + 
##     L(Consumption, 1) + L(Consumption, 2) + L(Consumption, 3) + 
##     L(Consumption, 4) + L(Consumption, 5) + L(Consumption, 6) + 
##     L(Consumption, 7) + L(Consumption, 8) + L(Consumption, 9) + 
##     L(Consumption, 10) + L(Consumption, 11) + L(Consumption, 
##     12) + L(Consumption, 13) + L(Consumption, 14) + L(Consumption, 
##     15), data = train.ts)
## 
## Residuals:
##         16         17         18         19         20         21         22 
##  0.0356177 -0.0324364 -0.0569201  0.0478216 -0.0140189  0.0601467  0.0126431 
##         23         24         25         26         27         28         29 
## -0.0585870 -0.0034234 -0.0795482  0.0030427  0.1227632  0.0119726 -0.0483270 
##         30         31         32         33         34         35         36 
##  0.0538041 -0.0170377 -0.0428077 -0.0672511  0.1092417 -0.0385621 -0.0278510 
##         37         38         39         40         41         42         43 
##  0.0969965 -0.0223928 -0.1524816  0.0826292  0.1303333 -0.1349213 -0.0003122 
##         44         45         46         47         48 
##  0.0590288 -0.0405853 -0.0157929  0.0355176 -0.0083020 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)
## (Intercept)         2.960e+01  1.518e+01   1.950    0.146
## Production          3.083e-07  4.236e-06   0.073    0.947
## L(Production, 1)   -3.499e-06  2.328e-06  -1.503    0.230
## L(Production, 2)    1.020e-05  8.034e-06   1.269    0.294
## L(Production, 3)    2.236e-06  5.723e-06   0.391    0.722
## L(Production, 4)   -1.008e-05  7.385e-06  -1.364    0.266
## L(Production, 5)    1.178e-05  9.451e-06   1.246    0.301
## L(Production, 6)    4.232e-06  6.037e-06   0.701    0.534
## L(Production, 7)   -2.112e-05  1.638e-05  -1.289    0.288
## L(Production, 8)    7.961e-06  9.836e-06   0.809    0.478
## L(Production, 9)    1.872e-05  1.341e-05   1.396    0.257
## L(Production, 10)  -9.476e-06  6.120e-06  -1.549    0.219
## L(Production, 11)  -1.367e-07  8.337e-06  -0.016    0.988
## L(Production, 12)   6.582e-06  7.802e-06   0.844    0.461
## L(Production, 13)  -9.126e-06  1.470e-05  -0.621    0.579
## L(Consumption, 1)   1.081e+00  5.781e-01   1.870    0.158
## L(Consumption, 2)  -1.515e+00  9.681e-01  -1.565    0.216
## L(Consumption, 3)  -5.349e-01  1.074e+00  -0.498    0.653
## L(Consumption, 4)   7.404e-01  8.022e-01   0.923    0.424
## L(Consumption, 5)  -2.039e+00  1.716e+00  -1.189    0.320
## L(Consumption, 6)   2.520e-01  1.097e+00   0.230    0.833
## L(Consumption, 7)   2.779e+00  2.417e+00   1.150    0.334
## L(Consumption, 8)  -2.214e+00  1.776e+00  -1.247    0.301
## L(Consumption, 9)  -1.733e+00  1.544e+00  -1.123    0.343
## L(Consumption, 10)  1.113e+00  7.610e-01   1.462    0.240
## L(Consumption, 11) -2.999e-01  1.014e+00  -0.296    0.787
## L(Consumption, 12) -1.067e+00  1.152e+00  -0.926    0.423
## L(Consumption, 13)  1.261e+00  1.879e+00   0.671    0.550
## L(Consumption, 14) -4.932e-01  4.604e-01  -1.071    0.363
## L(Consumption, 15) -3.478e-01  3.141e-01  -1.108    0.349
## 
## Residual standard error: 0.2197 on 3 degrees of freedom
## Multiple R-squared:  0.9997, Adjusted R-squared:  0.9971 
## F-statistic:   377 on 29 and 3 DF,  p-value: 0.0001931

Setelah dilakukan Pencarian model dengan menggunakan library dynlm, ternyata model yang didapat sama dengan model ketika kita menggunakan fungsi ARDL ataupun DLM model sendiri. Hal ini mengindikasikan bahwa kita dapat memakai library dynlm ataupun langsung menerapkan linear model sesuai dengan model yang ingin kita buat baik DLM ataupun ARDL.

deviance(cons_lm1)

## [1] 8.924194

deviance(cons_lm2)

## [1] 1.964576

deviance(cons_lm3)

## [1] 3.244475

deviance(cons_lm4)

## [1] 0.1447794

Pengujian Antar Model

Pada Kasus ini dataset antar model harus sama. Hal ini menjadi pertimbangan karena saat saya membangkitkan model-model, terdapat perbedaan data set , dimana ada yang lag nya 2 dan lagnya 13. Tentunya tidak akan bisa dibandingkan secara langsung. Oleh karena itu saya mencoba membandingan data dengan lag variabel x 2 dengan lag variabel y = 2

cons_lm1compare <- dynlm(Consumption ~ Production + L(Production,1) + L(Production,2), data = train.ts)
cons_lm2compare <- dynlm(Consumption ~ Production + L(Consumption,1) + L(Consumption,2), data = train.ts)

encomptest(cons_lm1compare,cons_lm2compare)

## Encompassing test
## 
## Model 1: Consumption ~ Production + L(Production, 1) + L(Production, 2)
## Model 2: Consumption ~ Production + L(Consumption, 1) + L(Consumption, 
##     2)
## Model E: Consumption ~ Production + L(Production, 1) + L(Production, 2) + 
##     L(Consumption, 1) + L(Consumption, 2)
##           Res.Df Df       F    Pr(>F)    
## M1 vs. ME     40 -2 35.0116 1.627e-09 ***
## M2 vs. ME     40 -2  8.4557 0.0008654 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari pengujian Encomptest, dilakukan pengujian untuk Data yang dimodelkan dengan 2 peubah lag pada Produksi dan 2 Peubah lag pada Konsumsi. Pada Encompassing test dibuat lagi Model E yang merupakan gabungan dari Model 1 dan Model 2. Artinya Pada Model ini Mau dibandingkan antara Autoregressive Model Vs Distributed Lag Model . Didapatkan bahwa P-value untuk ME ( Model Autoregressive) < 0.05 yang artinya model ME / Autoregressive lebih baik dan memberikan informasi yang lebih banyak daripada M1 ataupun M2.

Pengujian Asumsi

dwtest(cons_lm1)

## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 0.50824, p-value = 1.093e-10
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm2)

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

dwtest(cons_lm3)

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

dwtest(cons_lm4)

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

Dapat dilihat bahwa model Untuk Model 1 dan 2 ( DLM lag2, DLM Lag 13 )semuanya berpotensi memiliki autokorelasi, sedangkan untuk yang Model Autoregressive, Tak Tolak H0 yang artinya kurang cukup bukti untuk menyatakan adanya autokorelasi.

Hal ini juga mengindikasikan bahwa model Autoregressive dapat mengatasi masalah autokorelasi.

bptest(cons_lm1)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 0.51323, df = 3, p-value = 0.916

bptest(cons_lm2)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 11.159, df = 14, p-value = 0.6735

bptest(cons_lm3)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm3
## BP = 10.317, df = 5, p-value = 0.06674

bptest(cons_lm4)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm4
## BP = 30.996, df = 29, p-value = 0.3656

Semua Model dengan Pengujian Breusch-Pagan tidak mengalami heteroskedastistitas. atau dalam kata lain model-model ini semuanya bersifat homoskedastis.

shapiro.test(residuals(cons_lm1))

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

shapiro.test(residuals(cons_lm2))

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

shapiro.test(residuals(cons_lm3))

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

shapiro.test(residuals(cons_lm4))

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

Semua model Normal, karena semua model memilkii p-value > 0.05

Model Comparison

akurasi <- matrix(c(mape.koyck, mape.dlm2, mape.dlmopt, mape.ardl22))
row.names(akurasi)<- c("Koyck","DLM 1","DLM OPT","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi

##                      MAPE
## Koyck          0.08192264
## DLM 1          0.06431693
## DLM OPT        0.04150938
## Autoregressive 0.05088317

Berdasarkan MAPE yang didapat, Mape terkecil didapat ketika DLM 2, disusul dengan Model Autoregressive. Ketika mendapatkan hasil ini , memang secara substansial akan langsung dikatakan MAPE terkecil didapat oleh model Distributed Lag Model 2. Namun Ketika dalam penginterpretasiannya, saya akan lebih menyarankan Model Autoregressive.

Model tidak memiliki autokorelasi dibanding dengan model- model sebelumnya yang memiliki autokorelasi.
Model memiliki peubah signfikan yang lebih banyak dan Memiliki AIC yang lebih kecil dibandingkan Model DLM 2. Sehingga Lebih sederhana dan intepretatif karena memiliki lebih banyak peubah prediktor yang signifkan terhadap model.

Sehingga Untuk Kesimpulannya model yang paling direkomendasikan adalah Model Autoregressive.

par(mfrow=c(1,1))
plot(test$Production, test$Consumption, type="b", col="black", ylim=c(5,50))
points(test$Production, fore.koyck$forecasts,col="red")
lines(test$Production, fore.koyck$forecasts,col="red")
points(test$Production, fore.dlm2$forecasts,col="blue")
lines(test$Production, fore.dlm2$forecasts,col="blue")
points(test$Production, fore.dlmopt$forecasts,col="orange")
lines(test$Production, fore.dlmopt$forecasts,col="orange")
points(test$Production, fore.ardl22$forecasts,col="green")
lines(test$Production, fore.ardl22$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)

Dapat dilihat dari grafik plot yang paling mendekati adalah model DLM 2 ( Model DLM Optimum) Namun tetap seperti kesimpulan sebelumnya, Model akan lebih direkomendasikan di Model Autoregressive. karena berpotensi DLM Optimum ini mengalami overfitted akibat autokorelasi yang belum terselesaikan.

Kesimpulan

Dari hasil di atas, dapat disimpulkan bahwa model autoregressive merupakan model terbaik dalam memprediksi konsumsi ikan di Indonesia. Hal ini dapat dilihat dari nilai MAPE ykecil dibandingkan dengan model lainnya. Selain itu, model autoregressive juga memenuhi asumsi statistik yang diperlukan, seperti uji Durbin-Watson, Breusch-Pagan dan Shapiro-Wilk. Namun, perlu diingat bahwa model autoregressive ini hanya merupakan salah satu dari banyak model yang dapat digunakan dalam memprediksi konsumsi ikan di Indonesia. Oleh karena itu, perlu dilakukan penelitian lebih lanjut untuk mengevaluasi model-model lainnya dan memilih model yang paling sesuai dengan data yang ada.