library(readxl)## Warning: package 'readxl' was built under R version 4.1.3library(dLagM) #bisa otomatis timeseries datanya## Warning: package 'dLagM' was built under R version 4.1.3## Loading required package: nardl## Warning: package 'nardl' was built under R version 4.1.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Loading required package: dynlm## Warning: package 'dynlm' was built under R version 4.1.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.1.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(dynlm) #data harus timeseries
library(MLmetrics) #MAPE## Warning: package 'MLmetrics' was built under R version 4.1.3##
## Attaching package: 'MLmetrics'## The following object is masked from 'package:dLagM':
##
## MAPE## The following object is masked from 'package:base':
##
## Recalllibrary(lmtest)## Warning: package 'lmtest' was built under R version 4.1.3library(car)## Warning: package 'car' was built under R version 4.1.2## Loading required package: carData## Warning: package 'carData' was built under R version 4.1.2library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## 'data.frame': 252 obs. of 2 variables:
## $ Yt: num 154 173 154 185 184 ...
## $ Xt: num 67.8 72.2 66.2 72.2 71.2 ...t <- seq(1,252, by=1)
df <- cbind(t, data)
str(df)## 'data.frame': 252 obs. of 3 variables:
## $ t : num 1 2 3 4 5 6 7 8 9 10 ...
## $ Yt: num 154 173 154 185 184 ...
## $ Xt: num 67.8 72.2 66.2 72.2 71.2 ...#Model Awal
model_awal <- lm(data);model_awal;summary(model_awal)##
## Call:
## lm(formula = data)
##
## Coefficients:
## (Intercept) Xt
## 5.412 2.473##
## Call:
## lm(formula = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.109 -18.810 -4.477 14.102 179.028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4118 33.9096 0.160 0.873
## Xt 2.4735 0.4827 5.124 5.99e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.01 on 250 degrees of freedom
## Multiple R-squared: 0.09504, Adjusted R-squared: 0.09142
## F-statistic: 26.25 on 1 and 250 DF, p-value: 5.99e-07bgtest(model_awal)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_awal
## LM test = 18.374, df = 1, p-value = 1.815e-05Didapatkan P-value sebesar 1.815e-05 dimana P-value < 0.05 pada signifikansi 95%. maka dapat dinyatakan bahwa cukup bukti bahwa terdapat autokorelasi pada model.
#SPLIT DATA
train<-df[1:176,]
test<-df[177:252,]train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(df)#MODEL KOYCK
model.koyck = dLagM::koyckDlm(x = train$Xt, y = train$Yt)
summary(model.koyck)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -262.931 -27.226 -2.827 27.602 198.034
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 650.95737 476.54157 1.366 0.174
## Y.1 0.08412 0.20481 0.411 0.682
## X.t -6.95352 6.35981 -1.093 0.276
##
## Residual standard error: 45.41 on 172 degrees of freedom
## Multiple R-Squared: -1.41, Adjusted R-squared: -1.438
## Wald test: 3.128 on 2 and 172 DF, p-value: 0.04633
##
## Diagnostic tests:
## NULL
##
## alpha beta phi
## Geometric coefficients: 710.7438 -6.953517 0.08411805AIC(model.koyck)## [1] 1837.082BIC(model.koyck)## [1] 1849.741#Ramalan
(fore.koyck <- forecast(model = model.koyck, x=test$Xt, h=76))## $forecasts
## [1] 197.2899 170.3766 183.7581 148.3778 147.1400 190.4954 199.3575 182.7192
## [9] 150.0288 166.4011 148.6561 162.8088 160.5226 189.8827 185.3989 138.0855
## [17] 172.3499 145.6798 157.3433 189.6153 187.1148 169.5207 159.3488 194.9991
## [25] 178.8758 184.4730 169.2985 150.6382 178.6210 158.3759 208.8243 164.3933
## [33] 176.3013 173.8262 192.7402 149.1333 164.5874 172.8408 157.8897 219.2137
## [41] 184.3894 179.7217 167.1604 183.4876 176.1691 163.3848 195.3386 199.7649
## [49] 189.7069 202.7679 193.4363 150.9303 180.3839 189.8150 209.7305 169.6847
## [57] 168.0545 197.4698 198.2058 198.2677 165.2437 181.5880 182.9628 209.1541
## [65] 211.3573 194.1588 166.6364 159.1062 188.0252 185.2427 176.3167 199.9032
## [73] 182.7651 207.3991 178.1805 179.1994
##
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 76)
##
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"fore.koyck$forecasts## [1] 197.2899 170.3766 183.7581 148.3778 147.1400 190.4954 199.3575 182.7192
## [9] 150.0288 166.4011 148.6561 162.8088 160.5226 189.8827 185.3989 138.0855
## [17] 172.3499 145.6798 157.3433 189.6153 187.1148 169.5207 159.3488 194.9991
## [25] 178.8758 184.4730 169.2985 150.6382 178.6210 158.3759 208.8243 164.3933
## [33] 176.3013 173.8262 192.7402 149.1333 164.5874 172.8408 157.8897 219.2137
## [41] 184.3894 179.7217 167.1604 183.4876 176.1691 163.3848 195.3386 199.7649
## [49] 189.7069 202.7679 193.4363 150.9303 180.3839 189.8150 209.7305 169.6847
## [57] 168.0545 197.4698 198.2058 198.2677 165.2437 181.5880 182.9628 209.1541
## [65] 211.3573 194.1588 166.6364 159.1062 188.0252 185.2427 176.3167 199.9032
## [73] 182.7651 207.3991 178.1805 179.1994test$Yt## [1] 151.00 241.25 187.25 234.75 219.25 118.50 145.75 159.25 170.50 167.50
## [11] 232.75 210.50 202.25 185.00 153.00 244.25 193.50 224.75 162.75 180.00
## [21] 156.25 168.00 167.25 170.75 178.25 150.00 200.50 184.00 223.00 208.75
## [31] 166.00 195.00 160.50 159.75 140.50 216.25 168.25 194.75 172.75 219.00
## [41] 149.25 154.50 199.25 154.50 153.25 230.00 161.75 142.25 179.75 126.50
## [51] 169.50 198.50 174.50 167.75 147.75 182.25 175.50 161.75 157.75 168.75
## [61] 191.50 219.15 155.25 189.75 127.50 224.50 234.25 227.75 199.50 155.50
## [71] 215.50 134.25 201.00 186.75 190.75 207.50str(test$Yt)## num [1:76] 151 241 187 235 219 ...#mape data testing
mape.koyck <- MAPE(fore.koyck$forecasts, test$Yt)
#akurasi data training
mape_train_koyck <- dLagM::GoF(model.koyck)["MAPE"]
c("MAPE_testing" = mape.koyck, "MAPE_training" = mape_train_koyck)## $MAPE_testing
## [1] 0.1970092
##
## $MAPE_training.MAPE
## [1] 0.1865072#REGRESSION WITH DISTRIBUTED LAG -> estimasi parameter menggunakan least square
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
## -47.94 -17.80 -4.09 13.77 180.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.66349 59.74855 1.534 0.126851
## x.t 1.84743 0.53638 3.444 0.000721 ***
## x.1 -0.57839 0.53558 -1.080 0.281701
## x.2 -0.03925 0.53567 -0.073 0.941673
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.39 on 170 degrees of freedom
## Multiple R-squared: 0.06877, Adjusted R-squared: 0.05234
## F-statistic: 4.185 on 3 and 170 DF, p-value: 0.006892
##
## AIC and BIC values for the model:
## AIC BIC
## 1 1664.168 1679.963AIC(model.dlm)## [1] 1664.168BIC(model.dlm)## [1] 1679.963#ramalan
(fore.dlm <- forecast(model = model.dlm, x=test$Xt, h=76)) ## $forecasts
## [1] 172.5712 182.2841 175.6130 186.4367 183.0266 171.4187 173.6579 178.9556
## [9] 185.8525 178.0712 184.5655 178.9262 181.0433 172.9909 177.2769 189.3354
## [17] 175.2312 185.9989 180.0618 172.7384 176.8052 181.1666 182.0005 171.4804
## [25] 179.5483 176.3163 180.9435 184.2999 174.9141 183.2782 167.7096 184.8519
## [33] 176.9409 178.8912 173.6099 187.1891 178.4571 177.9451 182.7882 164.8992
## [41] 180.6391 178.1285 180.9912 175.3506 179.0372 181.6453 171.8088 174.0256
## [49] 177.1278 172.5752 176.4442 186.7397 174.4352 175.0996 170.7840 183.4984
## [57] 179.6742 171.7317 174.6615 174.6837 183.4492 175.6214 177.0255 170.2056
## [65] 172.3746 177.1404 182.6222 181.7407 173.3082 177.1224 179.1648 171.9463
## [73] 179.0021 170.6210 180.9955 177.6162
##
## $call
## forecast.dlm(model = model.dlm, x = test$Xt, h = 76)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"#mape data testing
mape.dlm <- MAPE(fore.dlm$forecasts, test$Yt)
#akurasi data training
mape_train <- GoF(model.dlm)["MAPE"]
c("MAPE_testing" = mape.dlm, "MAPE_training" = mape_train)## $MAPE_testing
## [1] 0.1294158
##
## $MAPE_training.MAPE
## [1] 0.1133169lag_optimum <-finiteDLMauto(formula=Yt~Xt,data=data.frame(train),model.type = "dlm",error.type = "AIC")
lag_optimum## q - k MASE AIC BIC GMRAE MBRAE R.Adj.Sq Ljung-Box
## 10 10 0.73151 1594.147 1634.603 0.69959 0.37386 0.08926 1.838307e-05model.dlm2 = dLagM::dlm(x = train$Xt,y = train$Yt , q = 10) #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
## -48.555 -16.917 -2.199 11.825 180.749
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -183.87074 102.54614 -1.793 0.07493 .
## x.t 1.67334 0.54618 3.064 0.00258 **
## x.1 -0.70990 0.54726 -1.297 0.19651
## x.2 -0.34410 0.54609 -0.630 0.52955
## x.3 0.61878 0.54539 1.135 0.25832
## x.4 0.43710 0.54558 0.801 0.42427
## x.5 0.79312 0.54490 1.456 0.14756
## x.6 -0.07687 0.54631 -0.141 0.88828
## x.7 1.13744 0.54611 2.083 0.03892 *
## x.8 0.92632 0.54475 1.700 0.09107 .
## x.9 0.06828 0.54388 0.126 0.90026
## x.10 0.62793 0.54234 1.158 0.24873
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.27 on 154 degrees of freedom
## Multiple R-squared: 0.15, Adjusted R-squared: 0.08926
## F-statistic: 2.47 on 11 and 154 DF, p-value: 0.007095
##
## AIC and BIC values for the model:
## AIC BIC
## 1 1594.147 1634.603AIC(model.dlm2)## [1] 1594.147BIC(model.dlm2)## [1] 1634.603#ramalan
(fore.dlm2 <- forecast(model = model.dlm2, x=test$Xt, h=76))## $forecasts
## [1] 167.4753 182.8617 169.3054 179.8663 180.7052 167.3103 177.4610 179.7931
## [9] 189.1872 176.6186 180.6871 186.9300 182.3284 176.1495 179.5158 196.7402
## [17] 178.7340 189.2157 188.5753 178.6150 186.4575 183.2178 192.0393 175.4871
## [25] 182.6733 189.1862 181.8032 184.1563 172.7941 184.8922 168.8513 185.0149
## [33] 182.6610 174.4270 180.7885 185.9220 184.4499 174.4458 181.8244 169.8849
## [41] 182.9202 182.2051 181.2774 180.2705 173.9239 188.7963 168.2019 169.3319
## [49] 179.5749 171.2219 176.1017 182.0646 173.0023 168.6329 165.2624 180.8748
## [57] 175.2066 159.6394 174.2629 174.9132 180.8207 169.8735 168.1753 170.7686
## [65] 167.4032 175.5403 174.9583 174.0144 165.0817 170.5741 176.2113 164.4887
## [73] 172.6455 167.5219 180.5674 176.6186
##
## $call
## forecast.dlm(model = model.dlm2, x = test$Xt, h = 76)
##
## attr(,"class")
## [1] "forecast.dlm" "dLagM"#akurasi testing
mape.dlm2 <- MAPE(fore.dlm2$forecasts, test$Yt)
#akurasi data training
mape_train_dlm <- GoF(model.dlm2)["MAPE"]
c("MAPE_testing" = mape.dlm2, "MAPE_training" = mape_train_dlm)## $MAPE_testing
## [1] 0.1315448
##
## $MAPE_training.MAPE
## [1] 0.1047278model.ardl = ardlDlm(x = train$Xt, y = train$Yt, p = 1 , q = 1) #p:lag x, q:lag y
summary(model.ardl)##
## Time series regression with "ts" data:
## Start = 2, End = 176
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -58.472 -15.329 -3.389 12.154 170.736
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.24868 48.36296 0.687 0.4927
## X.t 2.40244 0.50753 4.734 4.61e-06 ***
## X.1 -1.29190 0.51438 -2.512 0.0129 *
## Y.1 0.37490 0.07279 5.150 7.11e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.35 on 171 degrees of freedom
## Multiple R-squared: 0.1933, Adjusted R-squared: 0.1792
## F-statistic: 13.66 on 3 and 171 DF, p-value: 4.964e-08AIC(model.ardl)## [1] 1647.568BIC(model.ardl)## [1] 1663.392(fore.ardl <- forecast(model = model.ardl, x=test$Xt, h=76))## $forecasts
## [1] 159.2020 178.1504 174.0350 188.0118 185.8685 170.3728 170.8360 177.9847
## [9] 188.2459 179.6726 186.6179 180.2635 181.9892 171.7799 175.8454 192.2941
## [17] 176.5270 187.9317 181.9119 171.4279 175.1129 181.5315 183.7111 170.3006
## [25] 178.6622 175.8418 181.4818 186.6956 175.2101 184.2027 165.3571 184.7981
## [33] 177.6378 179.0614 172.3425 188.9921 180.2300 178.0954 183.9925 161.6746
## [41] 178.7487 178.3220 182.0433 175.1715 178.8281 182.7883 170.6006 171.5674
## [49] 175.8564 170.7217 174.9841 189.0588 175.1725 173.7006 167.8340 183.6019
## [57] 181.1613 170.3590 172.4003 172.8427 184.4201 176.0172 176.4197 167.5615
## [65] 169.0852 175.6625 183.9077 183.9560 172.7949 175.9030 179.1023 170.2783
## [73] 178.0986 168.4685 180.5138 178.0148
##
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 76)
##
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"#akurasi testing
mape.ardl <- MAPE(fore.ardl$forecasts, test$Yt) #data testing
#akurasi data training
mape_train_ardl <- GoF(model.ardl)["MAPE"]
c("MAPE_testing" = mape.ardl, "MAPE_training" = mape_train_ardl)## $MAPE_testing
## [1] 0.1258184
##
## $MAPE_training.MAPE
## [1] 0.0982388ardlBoundOrders(data = data.frame(df) , formula = Yt ~ Xt )## $p
## Xt
## 1 15
##
## $q
## [1] 15
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7 q = 8
## p = 1 2349.429 2338.060 2330.653 2324.208 2313.239 2306.886 2298.410 2288.703
## p = 2 2341.537 2339.871 2332.554 2326.093 2315.141 2308.789 2300.387 2290.703
## p = 3 2332.528 2332.528 2332.588 2325.959 2315.631 2309.095 2300.725 2290.337
## p = 4 2328.082 2325.268 2325.268 2327.238 2317.229 2310.628 2302.680 2292.299
## p = 5 2321.371 2319.396 2319.617 2319.617 2319.175 2312.598 2304.614 2293.227
## p = 6 2310.422 2307.296 2308.331 2310.110 2310.110 2312.101 2304.560 2293.678
## p = 7 2303.530 2302.061 2301.993 2303.941 2304.134 2304.134 2303.874 2293.678
## p = 8 2292.480 2290.581 2290.898 2292.387 2293.761 2295.737 2295.737 2295.245
## p = 9 2282.605 2282.088 2282.699 2284.604 2285.703 2287.577 2287.938 2287.938
## p = 10 2275.147 2274.168 2274.325 2276.081 2276.510 2278.333 2279.559 2280.086
## p = 11 2269.576 2269.473 2269.895 2271.842 2272.384 2274.374 2275.036 2275.824
## p = 12 2263.657 2263.653 2263.561 2265.502 2265.522 2267.518 2267.349 2268.065
## p = 13 2255.271 2255.112 2256.000 2258.000 2258.303 2260.283 2259.613 2259.751
## p = 14 2245.170 2244.582 2245.343 2247.288 2247.091 2248.995 2249.073 2249.313
## p = 15 2236.417 2235.432 2236.723 2238.631 2239.130 2241.112 2241.422 2241.844
## q = 9 q = 10 q = 11 q = 12 q = 13 q = 14 q = 15
## p = 1 2279.723 2270.989 2263.813 2257.365 2248.507 2241.483 2234.092
## p = 2 2281.715 2272.979 2265.794 2259.346 2250.457 2243.457 2236.078
## p = 3 2280.860 2272.452 2265.423 2258.949 2250.069 2243.150 2235.373
## p = 4 2282.669 2274.090 2267.103 2260.606 2251.771 2244.833 2237.048
## p = 5 2283.531 2275.331 2268.437 2261.974 2253.023 2246.090 2238.382
## p = 6 2285.069 2276.894 2269.856 2263.428 2254.588 2247.697 2239.940
## p = 7 2285.395 2278.151 2271.128 2264.760 2255.691 2248.878 2241.139
## p = 8 2286.799 2279.444 2271.934 2265.569 2256.934 2250.236 2242.375
## p = 9 2286.436 2279.380 2272.106 2265.425 2256.815 2249.793 2241.473
## p = 10 2280.086 2281.296 2274.057 2267.348 2258.780 2251.763 2243.472
## p = 11 2274.805 2274.805 2275.905 2269.146 2260.672 2253.762 2245.471
## p = 12 2267.666 2269.126 2269.126 2270.912 2262.545 2255.669 2247.426
## p = 13 2258.900 2260.411 2261.508 2261.508 2262.948 2255.816 2247.343
## p = 14 2248.815 2250.737 2252.462 2254.327 2254.327 2256.153 2247.062
## p = 15 2241.371 2243.179 2244.922 2246.658 2247.914 2247.914 2248.447
##
## $min.Stat
## [1] 2234.092mape <- matrix(c(mape_train_koyck, mape_train_dlm, mape_train_ardl))
row.names(mape)<- c("Koyck","DLM","ARDLM")
colnames(mape) <- c("MAPE")
mape## MAPE
## Koyck 0.1865072
## DLM 0.1047278
## ARDLM 0.0982388mapetest <- matrix(c(mape.koyck, mape.dlm2, mape.ardl))
row.names(mapetest)<- c("Koyck","DLM","ARDLM")
colnames(mapetest) <- c("MAPE")
mapetest## MAPE
## Koyck 0.1970092
## DLM 0.1315448
## ARDLM 0.1258184MAPE terkecil terdapat pada ARDLM dimana hal tersebut menunjukkan bahwa tingkat akurasi dari ARDLM sebesar 90%+. Dapat disimpulkan, bahwa model pada metode ARDLM merupakan model yang terbaik.
#PEMODELAN DLM dan ARDL dengan library dynlm
#library(dynlm)
#sama dengan model dlm p=1
cons_lm1 <- dynlm(Yt ~ Xt+L(Xt),data = train.ts)
#sama dengan model ardl p=0 q=1
cons_lm2 <- dynlm(Yt ~ Xt+L(Yt),data = train.ts)
#sama dengan ardl p=1 q=1
cons_lm3 <- dynlm(Yt ~ Xt+L(Xt)+L(Yt),data = train.ts)
#sama dengan dlm p=2
cons_lm4 <- dynlm(Yt ~ Xt+L(Xt)+L(Xt,2),data = train.ts)bgtest(cons_lm1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: cons_lm1
## LM test = 23.384, df = 1, p-value = 1.327e-06#tidak ada autokol pada model dengan ardl p = 0, q = 1
bgtest(cons_lm2)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: cons_lm2
## LM test = 3.0832, df = 1, p-value = 0.0791bgtest(cons_lm3)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: cons_lm3
## LM test = 15.867, df = 1, p-value = 6.796e-05Tidak terdapat autokolerasi metode ardlm pada model p = 0, q =1, setelah didapat P-value (0.0791) > 0.05 pada tingkat signifikasi 95%.