Tugas 3 MPDW
Farid Yafi Suwandi
2024-09-09
Packages
## 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.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## 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## Warning: package 'lmtest' was built under R version 4.3.3## Loading required package: carDataImpor Data
data <- read.csv("D:/STATATISTIKA/SEM5/MPDW/Data MPDW Tugas 3.csv")
str(data)## 'data.frame': 56 obs. of 3 variables:
## $ Periode : chr "1/1/2010" "4/1/2010" "7/1/2010" "10/1/2010" ...
## $ Tingkat.Pengangguran....: num 9.83 9.63 9.47 9.5 9.03 ...
## $ PDB.Riil..Milyar.USD. : num 16583 16743 16872 16961 16921 ...data## Periode Tingkat.Pengangguran.... PDB.Riil..Milyar.USD.
## 1 1/1/2010 9.833333 16582.71
## 2 4/1/2010 9.633333 16743.16
## 3 7/1/2010 9.466667 16872.27
## 4 10/1/2010 9.500000 16960.86
## 5 1/1/2011 9.033333 16920.63
## 6 4/1/2011 9.066667 17035.11
## 7 7/1/2011 9.000000 17031.31
## 8 10/1/2011 8.633333 17222.58
## 9 1/1/2012 8.266667 17367.01
## 10 4/1/2012 8.200000 17444.53
## 11 7/1/2012 8.033333 17469.65
## 12 10/1/2012 7.800000 17489.85
## 13 1/1/2013 7.733333 17662.40
## 14 4/1/2013 7.533333 17709.67
## 15 7/1/2013 7.233333 17860.45
## 16 10/1/2013 6.933333 18016.15
## 17 1/1/2014 6.666667 17953.97
## 18 4/1/2014 6.200000 18185.91
## 19 7/1/2014 6.066667 18406.94
## 20 10/1/2014 5.700000 18500.03
## 21 1/1/2015 5.533333 18666.62
## 22 4/1/2015 5.433333 18782.24
## 23 7/1/2015 5.100000 18857.42
## 24 10/1/2015 5.033333 18892.21
## 25 1/1/2016 4.900000 19001.69
## 26 4/1/2016 4.933333 19062.71
## 27 7/1/2016 4.900000 19197.94
## 28 10/1/2016 4.766667 19304.35
## 29 1/1/2017 4.566667 19398.34
## 30 4/1/2017 4.366667 19506.95
## 31 7/1/2017 4.333333 19660.77
## 32 10/1/2017 4.166667 19882.35
## 33 1/1/2018 4.033333 20044.08
## 34 4/1/2018 3.933333 20150.48
## 35 7/1/2018 3.766667 20276.15
## 36 10/1/2018 3.833333 20304.87
## 37 1/1/2019 3.866667 20415.15
## 38 4/1/2019 3.633333 20584.53
## 39 7/1/2019 3.600000 20817.58
## 40 10/1/2019 3.600000 20951.09
## 41 1/1/2020 3.833333 20665.55
## 42 4/1/2020 13.000000 19034.83
## 43 7/1/2020 8.800000 20511.78
## 44 10/1/2020 6.733333 20724.13
## 45 1/1/2021 6.233333 20990.54
## 46 4/1/2021 5.933333 21309.54
## 47 7/1/2021 5.066667 21483.08
## 48 10/1/2021 4.166667 21847.60
## 49 1/1/2022 3.800000 21738.87
## 50 4/1/2022 3.633333 21708.16
## 51 7/1/2022 3.533333 21851.13
## 52 10/1/2022 3.566667 21989.98
## 53 1/1/2023 3.500000 22112.33
## 54 4/1/2023 3.566667 22225.35
## 55 7/1/2023 3.700000 22490.69
## 56 10/1/2023 3.733333 22679.26Pembagian Data
#SPLIT DATA
train<-data[1:45,]
test<-data[46:56,]#data time series
train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(data)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.
Model kyock merupakan jenis paling umum dari model infinite distributed lag dan juga dikenal sebagai geometric lag
\[ y_t=a(1-\lambda)+\beta_0X_t+\beta_1Z_t+\lambda Y_{t-1}+V_t \]
dengan \[V_t=u_t-\lambda u_{t-1}\]
Pemodelan
Pemodelan model Koyck dengan R dapat menggunakan
dLagM::koyckDlm() . Fungsi umum dari koyckDlm
adalah sebagai berikut.
koyckDlm(x , y , intercept)Fungsi koyckDlm() akan menerapkan model lag
terdistribusi dengan transformasi Koyck satu prediktor. Nilai
x dan y tidak perlu sebagai objek time
series (ts). intercept dapat dibuat
TRUE untuk memasukkan intersep ke dalam model.
#MODEL KOYCK
model.koyck <- koyckDlm(x = data$Tingkat.Pengangguran...., y = data$PDB.Riil..Milyar.USD.)
summary(model.koyck)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -3144.33 -116.78 78.39 221.49 912.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.706e+03 1.946e+03 -1.904 0.0625 .
## Y.1 1.146e+00 7.979e-02 14.365 <2e-16 ***
## X.t 1.692e+02 7.630e+01 2.217 0.0310 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 494.9 on 52 degrees of freedom
## Multiple R-Squared: 0.9224, Adjusted R-squared: 0.9194
## Wald test: 326.3 on 2 and 52 DF, p-value: < 2.2e-16
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: 25358.77 169.1856 1.146128AIC(model.koyck)## [1] 843.4784BIC(model.koyck)## [1] 851.5078Dari 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\). Adapun model keseluruhannya adalah sebagai berikut
\[ \hat{Y_t}=-3706+ 169.2X_t+1.146Y_{t-1} \]
Peramalan dan Akurasi
Berikut adalah hasil peramalan y untuk 5 periode kedepan menggunakan model koyck
fore.koyck <- forecast(model = model.koyck, x=test$Tingkat.Pengangguran...., h=11)
fore.koyck## $forecasts
## [1] 23291.54 23846.66 24330.64 24823.31 25359.77 25957.71 26648.66 27429.30
## [9] 28335.30 29396.24 30617.85
##
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Tingkat.Pengangguran....,
## h = 11)
##
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"mape.koyck <- MAPE(fore.koyck$forecasts, test$PDB.Riil..Milyar.USD.)
#akurasi data training
GoF(model.koyck)## n MAE MPE MAPE sMAPE MASE MSE
## model.koyck 55 240.6466 -0.0003744553 0.01253489 0.01233398 1.266573 231570.5
## MRAE GMRAE
## model.koyck 3.283851 1.256227Regression with Distributed Lag
Pemodelan model Regression with Distributed Lag dengan R
dapat menggunakan dLagM::dlm() . Fungsi umum dari
dlm adalah sebagai berikut.
dlm(formula , data , x , y , q , remove )Fungsi dlm() akan menerapkan model lag terdistribusi
dengan satu atau lebih prediktor. Nilai x dan
y tidak perlu sebagai objek time series
(ts). \(q\) adalah integer
yang mewakili panjang lag yang terbatas.
Pemodelan (Lag=2)
model.dlm <- dlm(x = train$Tingkat.Pengangguran....,y = train$PDB.Riil..Milyar.USD. , q = 2)
summary(model.dlm)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -626.9 -492.2 -398.2 191.9 3349.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21952.905 455.387 48.207 < 2e-16 ***
## x.t -280.325 94.898 -2.954 0.00529 **
## x.1 -2.947 110.689 -0.027 0.97890
## x.2 -204.994 92.045 -2.227 0.03178 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 897.5 on 39 degrees of freedom
## Multiple R-squared: 0.566, Adjusted R-squared: 0.5326
## F-statistic: 16.96 on 3 and 39 DF, p-value: 3.317e-07
##
## AIC and BIC values for the model:
## AIC BIC
## 1 712.595 721.401AIC(model.dlm)## [1] 712.595BIC(model.dlm)## [1] 721.401Dari hasil diatas, didapat bahwa \(P-value\) dari intercept dan \(x_{t-1}<0.05\). Hal ini menunjukkan bahwa intercept dan \(x_{t-1}\) berpengaruh signifikan terhadap \(y\). Adapun model keseluruhan yang terbentuk adalah sebagai berikut
\[ \hat{Y_t}=21952.905-280.325X_t-2.947X_{t-1}-204.994X_{t-2} \]
Peramalan dan Akurasi
Berikut merupakan hasil peramalan \(y\) untuk 5 periode kedepan
fore.dlm <- forecast(model = model.dlm, x=test$Tingkat.Pengangguran...., h=11)
fore.dlm## $forecasts
## [1] 18890.99 19237.32 19553.66 19836.76 20069.05 20172.74 20197.86 20236.95
## [9] 20211.62 20187.72 20164.31
##
## $call
## forecast.dlm(model = model.dlm, x = test$Tingkat.Pengangguran....,
## h = 11)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"mape.dlm <- MAPE(fore.dlm$forecasts, test$PDB.Riil..Milyar.USD.)
#akurasi data training
GoF(model.dlm)## n MAE MPE MAPE sMAPE MASE MSE
## model.dlm 43 594.9688 -0.001989484 0.0309996 0.03152953 3.061138 730541.1
## MRAE GMRAE
## model.dlm 8.683747 3.718007Lag Optimum
#penentuan lag optimum
finiteDLMauto(formula = PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran....,
data = data.frame(train), q.min = 1, q.max = 6,
model.type = "dlm", error.type = "AIC", trace = FALSE)## q - k MASE AIC BIC GMRAE MBRAE R.Adj.Sq Ljung-Box
## 6 6 0.73549 541.422 556.394 0.66697 0.58874 0.96626 5.787971e-08finiteDLMauto## function (formula, data, x, y, q.min = 1, q.max = 10, k.order = NULL,
## model.type = c("dlm", "poly"), error.type = c("MASE", "AIC",
## "BIC", "GMRAE", "MBRAE", "radj"), trace = FALSE, type)
## UseMethod("finiteDLMauto")
## <bytecode: 0x0000024150060fc8>
## <environment: namespace:dLagM>Berdasarkan output tersebut, lag optimum didapatkan ketika lag=6. Selanjutnya dilakukan pemodelan untuk lag=6
#model dlm dengan lag optimum
model.dlm2 <- dlm(x = train$Tingkat.Pengangguran....,y = train$PDB.Riil..Milyar.USD. , q = 6)
summary(model.dlm2)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -401.55 -84.93 16.34 134.28 463.71
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22794.71 134.59 169.367 < 2e-16 ***
## x.t -156.19 25.19 -6.201 6.93e-07 ***
## x.1 85.16 28.72 2.965 0.00578 **
## x.2 13.96 28.28 0.494 0.62500
## x.3 44.61 29.37 1.519 0.13895
## x.4 31.47 294.06 0.107 0.91546
## x.5 -161.71 385.83 -0.419 0.67801
## x.6 -454.17 284.24 -1.598 0.12023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 222.8 on 31 degrees of freedom
## Multiple R-squared: 0.9725, Adjusted R-squared: 0.9663
## F-statistic: 156.4 on 7 and 31 DF, p-value: < 2.2e-16
##
## AIC and BIC values for the model:
## AIC BIC
## 1 541.422 556.394AIC(model.dlm2)## [1] 541.422BIC(model.dlm2)## [1] 556.394Dari hasil tersebut terdapat beberapa peubah yang berpengaruh signifikan terhadap taraf nyata 5% yaitu \(x_t\) , \(x_{t-1}\). Adapun keseluruhan model yang terbentuk adalah
\[ \hat{Y_t}=22794.71-156.19X_t-3.85.16X_{t-1} \]
Adapun hasil peramalan 5 periode kedepan menggunakan model tersebut adalah sebagai berikut
#peramalan dan akurasi
fore.dlm2 <- forecast(model = model.dlm2, x=test$Tingkat.Pengangguran...., h=11)
mape.dlm2<- MAPE(fore.dlm2$forecasts, test$PDB.Riil..Milyar.USD.)
#akurasi data training
GoF(model.dlm2)## n MAE MPE MAPE sMAPE MASE MSE
## model.dlm2 39 153.2149 -8.681193e-05 0.007937706 0.007929914 0.7354851 39466.17
## MRAE GMRAE
## model.dlm2 1.619134 0.6669745Model tersebut merupakan model yang sangat baik dengan nilai MAPE yang kurang dari 10%.
Model Autoregressive
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).
Pemodelan
Pemodelan Autoregressive dilakukan menggunakan fungsi
dLagM::ardlDlm() . Fungsi tersebut akan menerapkan
autoregressive berordo \((p,q)\) dengan satu prediktor. Fungsi umum
dari ardlDlm() adalah sebagai berikut.
ardlDlm(formula = NULL , data = NULL , x = NULL , y = NULL , p = 1 , q = 1 ,
remove = NULL )Dengan \(p\) adalah integer yang mewakili panjang lag yang terbatas dan \(q\) adalah integer yang merepresentasikan ordo dari proses autoregressive.
#model.ardl <- ardlDlm(x = train$Tingkat.Pengangguran...., y = train$PDB.Riil..Milyar.USD., p = 1 , q = 1)
#summary(model.ardl)
#AIC(model.ardl)
#BIC(model.ardl)model.ardl <- ardlDlm(formula = PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran....,
data = train,p = 1 , q = 1)
summary(model.ardl)##
## Time series regression with "ts" data:
## Start = 2, End = 45
##
## Call:
## dynlm(formula = as.formula(model.text), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -409.22 -35.73 14.69 46.85 387.13
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -804.72732 494.97673 -1.626 0.112
## Tingkat.Pengangguran.....t -197.53374 13.69436 -14.424 <2e-16 ***
## Tingkat.Pengangguran.....1 222.39356 15.36315 14.476 <2e-16 ***
## PDB.Riil..Milyar.USD..1 1.03897 0.02262 45.927 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 128.9 on 40 degrees of freedom
## Multiple R-squared: 0.9913, Adjusted R-squared: 0.9907
## F-statistic: 1528 on 3 and 40 DF, p-value: < 2.2e-16AIC(model.ardl)## [1] 558.2844BIC(model.ardl)## [1] 567.2053Hasil di atas menunjukkan bahwa selain peubah \(x_{t-1}\), hasil uji t menunjukkan nilai-p pada peubah \(\ge0.05\) Hal ini menunjukkan bahwa peubah \(x_{t-1}\) berpengaruh signifikan terhadap \(y_t\), sementara \(x_t\) dan \(y_{t-1}\) berpengaruh signifikan terhadap \(y_t\). Model keseluruhannya adalah sebagai berikut:
\[ \hat{Y}=-804.72732-197.53374X_t+222.39356X_{t-1}+1.03897Y_{t-1} \]
Peramalan dan Akurasi
fore.ardl <- forecast(model = model.ardl, x=test$Tingkat.Pengangguran...., h=11)
fore.ardl## $forecasts
## [1] 21217.98 21558.76 21897.86 22122.44 22307.16 22481.76 22634.34 22813.45
## [9] 22971.54 23124.28 23306.04
##
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Tingkat.Pengangguran....,
## h = 11)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"Data di atas merupakan hasil peramalan untuk 5 periode ke depan menggunakan Model Autoregressive dengan \(p=1\) dan \(q=1\).
mape.ardl <- MAPE(fore.ardl$forecasts, test$PDB.Riil..Milyar.USD.)
mape.ardl## [1] 0.02132807#akurasi data training
GoF(model.ardl)## n MAE MPE MAPE sMAPE MASE MSE
## model.ardl 44 78.23707 -4.780592e-05 0.00409586 0.00409115 0.4057006 15110.28
## MRAE GMRAE
## model.ardl 1.391476 0.3721362Berdasarkan akurasi di atas, terlihat bahwa nilai MAPE keduanya tidak
jauh berbeda. Artinya, model regresi dengan distribusi lag ini
tidakoverfittedatauunderfitted
Lag Optimum
#penentuan lag optimum
model.ardl.opt <- ardlBoundOrders(data = data.frame(data), ic = "AIC",
formula = PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran.... )
min_p=c()
for(i in 1:6){
min_p[i]=min(model.ardl.opt$Stat.table[[i]])
}
q_opt=which(min_p==min(min_p, na.rm = TRUE))
p_opt=which(model.ardl.opt$Stat.table[[q_opt]] ==
min(model.ardl.opt$Stat.table[[q_opt]], na.rm = TRUE))
data.frame("q_optimum" = q_opt, "p_optimum" = p_opt,
"AIC"=model.ardl.opt$min.Stat)## q_optimum p_optimum AIC
## 1 4 15 461.2204Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat
ketika \(p=15\) dan \(q=4\), yaitu sebesar 461.2204.
Artinya, model autoregressive optimum didapat ketika \(p=15\) dan \(q=4\).
Selanjutnya dapat dilakukan pemodelan dengan nilai \(p\) dan \(q\) optimum seperti inisialisasi di langkah sebelumnya.
Pemodelan DLM & ARDL dengan Library dynlm
Pemodelan regresi dengan peubah lag tidak hanya dapat
dilakukan dengan fungsi pada packages dLagM ,
tetapi terdapat packages dynlm yang dapat
digunakan. Fungsi dynlm secara umum adalah sebagai
berikut.
dynlm(formula, data, subset, weights, na.action, method = "qr",
model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE,
contrasts = NULL, offset, start = NULL, end = NULL, ...)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=1
cons_lm1 <- dynlm(PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran....+L(Tingkat.Pengangguran....),data = train.ts)
#sama dengan model ardl p=1 q=0
cons_lm2 <- dynlm(PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran....+L(PDB.Riil..Milyar.USD.),data = train.ts)
#sama dengan ardl p=1 q=1
cons_lm3 <- dynlm(PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran....+L(Tingkat.Pengangguran....)+L(PDB.Riil..Milyar.USD.),data = train.ts)
#sama dengan dlm p=2
cons_lm4 <- dynlm(PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran....+L(Tingkat.Pengangguran....)+L(Tingkat.Pengangguran....,2),data = train.ts)Ringkasan Model
summary(cons_lm1)##
## Time series regression with "ts" data:
## Start = 2, End = 45
##
## Call:
## dynlm(formula = PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran.... +
## L(Tingkat.Pengangguran....), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -684.6 -527.7 -459.5 243.2 3433.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21748.38 449.96 48.334 < 2e-16 ***
## Tingkat.Pengangguran.... -313.76 97.44 -3.220 0.00251 **
## L(Tingkat.Pengangguran....) -146.84 94.79 -1.549 0.12904
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 933.5 on 41 degrees of freedom
## Multiple R-squared: 0.5351, Adjusted R-squared: 0.5124
## F-statistic: 23.59 on 2 and 41 DF, p-value: 1.517e-07summary(cons_lm2)##
## Time series regression with "ts" data:
## Start = 2, End = 45
##
## Call:
## dynlm(formula = PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran.... +
## L(PDB.Riil..Milyar.USD.), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -863.25 -90.30 -27.21 35.09 1644.73
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3156.51698 1017.56564 3.102 0.00347 **
## Tingkat.Pengangguran.... -91.38064 28.53311 -3.203 0.00263 **
## L(PDB.Riil..Milyar.USD.) 0.86760 0.04756 18.243 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 318.1 on 41 degrees of freedom
## Multiple R-squared: 0.946, Adjusted R-squared: 0.9434
## F-statistic: 359.3 on 2 and 41 DF, p-value: < 2.2e-16summary(cons_lm3)##
## Time series regression with "ts" data:
## Start = 2, End = 45
##
## Call:
## dynlm(formula = PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran.... +
## L(Tingkat.Pengangguran....) + L(PDB.Riil..Milyar.USD.), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -409.22 -35.73 14.69 46.85 387.13
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -804.72732 494.97673 -1.626 0.112
## Tingkat.Pengangguran.... -197.53374 13.69436 -14.424 <2e-16 ***
## L(Tingkat.Pengangguran....) 222.39356 15.36315 14.476 <2e-16 ***
## L(PDB.Riil..Milyar.USD.) 1.03897 0.02262 45.927 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 128.9 on 40 degrees of freedom
## Multiple R-squared: 0.9913, Adjusted R-squared: 0.9907
## F-statistic: 1528 on 3 and 40 DF, p-value: < 2.2e-16summary(cons_lm4)##
## Time series regression with "ts" data:
## Start = 3, End = 45
##
## Call:
## dynlm(formula = PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran.... +
## L(Tingkat.Pengangguran....) + L(Tingkat.Pengangguran....,
## 2), data = train.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -626.9 -492.2 -398.2 191.9 3349.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21952.905 455.387 48.207 < 2e-16 ***
## Tingkat.Pengangguran.... -280.325 94.898 -2.954 0.00529 **
## L(Tingkat.Pengangguran....) -2.947 110.689 -0.027 0.97890
## L(Tingkat.Pengangguran...., 2) -204.994 92.045 -2.227 0.03178 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 897.5 on 39 degrees of freedom
## Multiple R-squared: 0.566, Adjusted R-squared: 0.5326
## F-statistic: 16.96 on 3 and 39 DF, p-value: 3.317e-07SSE
deviance(cons_lm1)## [1] 35724616deviance(cons_lm2)## [1] 4147818deviance(cons_lm3)## [1] 664852.5deviance(cons_lm4)## [1] 31413269Uji Diagnostik
#uji model
if(require("lmtest")) encomptest(cons_lm1, cons_lm2)## Encompassing test
##
## Model 1: PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran.... + L(Tingkat.Pengangguran....)
## Model 2: PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran.... + L(PDB.Riil..Milyar.USD.)
## Model E: PDB.Riil..Milyar.USD. ~ Tingkat.Pengangguran.... + L(Tingkat.Pengangguran....) +
## L(PDB.Riil..Milyar.USD.)
## Res.Df Df F Pr(>F)
## M1 vs. ME 40 -1 2109.33 < 2.2e-16 ***
## M2 vs. ME 40 -1 209.55 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Autokorelasi
#durbin watson
dwtest(cons_lm1)##
## Durbin-Watson test
##
## data: cons_lm1
## DW = 0.15092, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm2)##
## Durbin-Watson test
##
## data: cons_lm2
## DW = 2.0845, p-value = 0.5139
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm3)##
## Durbin-Watson test
##
## data: cons_lm3
## DW = 2.5299, p-value = 0.9363
## alternative hypothesis: true autocorrelation is greater than 0dwtest(cons_lm4)##
## Durbin-Watson test
##
## data: cons_lm4
## DW = 0.13439, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0Heterogenitas
bptest(cons_lm1)##
## studentized Breusch-Pagan test
##
## data: cons_lm1
## BP = 7.7704, df = 2, p-value = 0.02054bptest(cons_lm2)##
## studentized Breusch-Pagan test
##
## data: cons_lm2
## BP = 10.9, df = 2, p-value = 0.004296bptest(cons_lm3)##
## studentized Breusch-Pagan test
##
## data: cons_lm3
## BP = 27.532, df = 3, p-value = 4.554e-06bptest(cons_lm4)##
## studentized Breusch-Pagan test
##
## data: cons_lm4
## BP = 7.192, df = 3, p-value = 0.06602Kenormalan
shapiro.test(residuals(cons_lm1))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm1)
## W = 0.68699, p-value = 2.126e-08shapiro.test(residuals(cons_lm2))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm2)
## W = 0.58827, p-value = 6.792e-10shapiro.test(residuals(cons_lm3))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm3)
## W = 0.85556, p-value = 6.072e-05shapiro.test(residuals(cons_lm4))##
## Shapiro-Wilk normality test
##
## data: residuals(cons_lm4)
## W = 0.67017, p-value = 1.465e-08Perbandingan Model
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.19991016
## DLM 1 0.09391342
## DLM 2 0.11740808
## Autoregressive 0.02132807Berdasarkan nilai MAPE, model paling optimum didapat pada Model Koyck karena memiliki nilai MAPE yang terkecil.
Plot
par(mfrow=c(1,1))
plot(test$Tingkat.Pengangguran...., test$PDB.Riil..Milyar.USD., type="b", col="black", ylim= c(14000,33000))
points(test$Tingkat.Pengangguran...., fore.koyck$forecasts,col="red")
lines(test$Tingkat.Pengangguran...., fore.koyck$forecasts,col="red")
points(test$Tingkat.Pengangguran...., fore.dlm$forecasts,col="blue")
lines(test$Tingkat.Pengangguran...., fore.dlm$forecasts,col="blue")
points(test$Tingkat.Pengangguran...., fore.dlm2$forecasts,col="orange")
lines(test$Tingkat.Pengangguran...., fore.dlm2$forecasts,col="orange")
points(test$Tingkat.Pengangguran...., fore.ardl$forecasts,col="green")
lines(test$Tingkat.Pengangguran...., fore.ardl$forecasts,col="green")
legend("topright",c("aktual", "koyck","DLM 1","DLM 2", "autoregressive"), lty=1, col=c("black","red","blue","orange","green"), cex=0.8)
Berdasarkan plot tersebut, terlihat bahwa plot yang paling mendekati data aktualnya adalah Model koyck, sehingga dapat disimpulkan model terbaik dalam hal ini adalah model regresi koyck