##Library yang digunakan

library(dLagM)
## 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.numeric
library(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':
## 
##     Recall
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.1.3
library(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.1
library(readxl)

##Data

data <- read.csv("DailyDelhiClimateTest.csv")
View(data)

#Split Data
train<-data[1:1170,]
test<-data[1171:1462,]

#Data Time Series
train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(data)

##Model Koyck

model.koyck = dLagM::koyckDlm(x = train$wind_speed, y = train$meanpressure)
summary(model.koyck)
## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.4762  -3.0153   0.5225   3.7603  11.7870 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.017e+03  6.237e+00 163.029   <2e-16 ***
## Y.1          9.471e-04  6.371e-03   0.149   0.8821    
## X.t         -6.488e-01  3.497e-01  -1.855   0.0662 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.804 on 110 degrees of freedom
## Multiple R-Squared: -0.01488,    Adjusted R-squared: -0.03333 
## Wald test: 1.801 on 2 and 110 DF,  p-value: 0.1699 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha       beta          phi
## Geometric coefficients:  1017.726 -0.6488312 0.0009471302
AIC(model.koyck)
## [1] 723.0849
BIC(model.koyck)
## [1] 733.9945

#Ramalan

(fore.koyck <- forecast(model = model.koyck, x=test$wind_speed, h=292))
## $forecasts
##   [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [51] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [101] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [126] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [151] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [176] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [201] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [226] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [251] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [276] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$wind_speed, h = 292)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

#MAPE Data Testing

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

#Akurasi Data Training

mape_train <- dLagM::GoF(model.koyck)["MAPE"]

c("MAPE_testing" = mape.koyck, "MAPE_training" = mape_train)
## $MAPE_testing
## [1] NA
## 
## $MAPE_training.MAPE
## [1] 0.004472461

##Regresi dengan Distribusi Lag

model.dlm = dLagM::dlm(x = train$wind_speed,y = train$meanpressure , q = 2)
summary(model.dlm)
## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.178  -3.509   0.127   4.128  10.605 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1016.0760     1.7212 590.336   <2e-16 ***
## x.t           -0.1974     0.1686  -1.171    0.244    
## x.1           -0.1251     0.1848  -0.677    0.500    
## x.2           -0.1339     0.1672  -0.801    0.425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.651 on 108 degrees of freedom
##   (1056 observations deleted due to missingness)
## Multiple R-squared:  0.04651,    Adjusted R-squared:  0.02003 
## F-statistic: 1.756 on 3 and 108 DF,  p-value: 0.1599
## 
## AIC and BIC values for the model:
##       AIC      BIC
## 1 711.689 725.2815
AIC(model.dlm)
## [1] 711.689
BIC(model.dlm)
## [1] 725.2815
#Ramalan
(fore.dlm <- forecast(model = model.dlm, x=test$wind_speed, h=292)) #meramal 292 periode ke depan
## $forecasts
##   [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [51] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [101] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [126] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [151] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [176] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [201] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [226] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [251] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [276] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 
## $call
## forecast.dlm(model = model.dlm, x = test$wind_speed, h = 292)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

#MAPE Data Testing

mape.dlm <- MAPE(fore.dlm$forecasts, test$meanpressure)

#Akurasi Data Training

mape_train <- GoF(model.dlm)["MAPE"]
c("MAPE_testing" = mape.dlm, "MAPE_training" = mape_train)
## $MAPE_testing
## [1] NA
## 
## $MAPE_training.MAPE
## [1] 0.004459747

#Penentuan Lag Optimum

finiteDLMauto(formula = meanpressure ~ wind_speed,
              data = data.frame(train), q.min = 1, q.max = 20 ,
              model.type = "dlm", error.type = "AIC", trace = TRUE)
##    q - k    MASE      AIC      BIC   GMRAE     MBRAE R.Adj.Sq    Ljung-Box
## 20    20 1.85801 548.6189 607.1147 2.40450   0.79779  0.49402 0.000000e+00
## 19    19 1.90639 554.6560 610.8413 2.44062   0.62222  0.50069 0.000000e+00
## 18    18 1.93839 559.5167 613.3680 2.41813   0.99817  0.51571 0.000000e+00
## 17    17 1.95582 562.8001 614.2943 2.59537   0.71161  0.53414 0.000000e+00
## 16    16 1.91118 566.2377 615.3521 2.42372   0.74214  0.54218 0.000000e+00
## 15    15 1.89365 573.9711 620.6833 2.11928   0.92212  0.52768 0.000000e+00
## 14    14 1.99932 586.3576 630.6455 2.30098   0.86201  0.49133 1.110223e-16
## 13    13 2.10370 601.6183 643.4603 2.26769   0.64950  0.43846 5.551115e-16
## 12    12 2.25251 619.1454 658.5200 2.63701   0.99496  0.36419 4.440892e-16
## 11    11 2.38785 637.1845 674.0707 2.55637 -35.97878  0.27568 1.110223e-16
## 10    10 2.60652 654.5040 688.8811 3.29621   1.01042  0.18322 0.000000e+00
## 9      9 2.71584 670.1054 701.9529 3.31891  -1.44335  0.09598 0.000000e+00
## 8      8 2.80405 679.2424 708.5403 3.23197   0.98873  0.05810 0.000000e+00
## 7      7 2.88748 685.4482 712.1765 3.44514   0.65168  0.04490 0.000000e+00
## 6      6 2.90046 690.7741 714.9133 3.51397   0.36190  0.03673 0.000000e+00
## 5      5 2.97863 696.2313 717.7621 3.68160   0.73000  0.02713 0.000000e+00
## 4      4 2.98148 701.6980 720.6013 3.77703  -0.40375  0.01847 0.000000e+00
## 3      3 2.99105 706.6783 722.9355 3.76636   0.84307  0.01609 0.000000e+00
## 2      2 2.99104 711.6890 725.2815 3.36729   4.99124  0.02003 0.000000e+00
## 1      1 3.02548 716.0901 726.9996 3.68938   0.93271  0.02870 0.000000e+00

##Model dlm dengan Lag Optimum

model.dlm2 = dLagM::dlm(x = train$wind_speed,y = train$meanpressure , q = 2) #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 
## -14.178  -3.509   0.127   4.128  10.605 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1016.0760     1.7212 590.336   <2e-16 ***
## x.t           -0.1974     0.1686  -1.171    0.244    
## x.1           -0.1251     0.1848  -0.677    0.500    
## x.2           -0.1339     0.1672  -0.801    0.425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.651 on 108 degrees of freedom
##   (1056 observations deleted due to missingness)
## Multiple R-squared:  0.04651,    Adjusted R-squared:  0.02003 
## F-statistic: 1.756 on 3 and 108 DF,  p-value: 0.1599
## 
## AIC and BIC values for the model:
##       AIC      BIC
## 1 711.689 725.2815
AIC(model.dlm2)
## [1] 711.689
BIC(model.dlm2)
## [1] 725.2815
#Ramalan
(fore.dlm2 <- forecast(model = model.dlm2, x=test$wind_speed, h=292))
## $forecasts
##   [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [51] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [101] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [126] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [151] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [176] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [201] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [226] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [251] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [276] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 
## $call
## forecast.dlm(model = model.dlm2, x = test$wind_speed, h = 292)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"
#akurasi testing
mape.dlm2 <- MAPE(fore.dlm2$forecasts, test$meanpressure)

#akurasi data training
mape_train <- GoF(model.dlm2)["MAPE"]
c("MAPE_testing" = mape.dlm2, "MAPE_training" = mape_train)
## $MAPE_testing
## [1] NA
## 
## $MAPE_training.MAPE
## [1] 0.004459747

##Model Autoregressive

model.ardl = ardlDlm(x = train$wind_speed, y = train$meanpressure, p = 1 , q = 2) #p:lag x, q:lag y
summary(model.ardl)
## 
## Time series regression with "ts" data:
## Start = 3, End = 114
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5163 -1.2938 -0.1179  1.2545  4.9277 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 59.5763471 33.5019697   1.778   0.0782 .  
## X.t         -0.0715131  0.0579386  -1.234   0.2198    
## X.1          0.0186724  0.0581577   0.321   0.7488    
## Y.1          0.9419248  0.0329084  28.623   <2e-16 ***
## Y.2         -0.0004643  0.0020553  -0.226   0.8217    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.935 on 107 degrees of freedom
##   (0 observations deleted due to missingness)
## Multiple R-squared:  0.8892, Adjusted R-squared:  0.8851 
## F-statistic: 214.8 on 4 and 107 DF,  p-value: < 2.2e-16
AIC(model.ardl)
## [1] 472.5776
BIC(model.ardl)
## [1] 488.8886
#Ramalam
(fore.ardl <- forecast(model = model.ardl, x=test$wind_speed, h=292))
## $forecasts
##   [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [51] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##  [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [101] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [126] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [151] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [176] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [201] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [226] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [251] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [276] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$wind_speed, h = 292)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

#Akurasi Testing

mape.ardl <- MAPE(fore.ardl$forecasts, test$meanpressure)

#Akurasi Data Training

mape_train <- GoF(model.ardl)["MAPE"]
c("MAPE_testing" = mape.ardl, "MAPE_training" = mape_train)
## $MAPE_testing
## [1] NA
## 
## $MAPE_training.MAPE
## [1] 0.001425909
ardlBoundOrders(data = data.frame(data) , formula = meanpressure ~ wind_speed )
## $p
##   wind_speed
## 1         15
## 
## $q
## [1] 2
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  473.5991 460.7036 451.3858 447.9351 446.4951 442.2671 440.2490 432.0854
## p = 2  460.5934 460.6039 451.8635 448.0050 446.5470 442.9552 441.1002 432.2163
## p = 3  451.9482 451.9482 453.8521 449.9607 448.4549 444.9035 443.0825 434.2162
## p = 4  451.3463 452.1143 452.1143 451.6035 450.1688 446.7385 444.8841 435.9404
## p = 5  449.8774 448.2091 450.1686 450.1686 452.1659 448.7043 446.8802 437.9330
## p = 6  451.9393 447.7393 449.5916 450.7988 450.7988 450.5786 448.7511 439.5049
## p = 7  447.4751 443.0913 444.9407 446.8074 448.3906 448.3906 449.8099 441.2011
## p = 8  439.0353 436.4102 438.4095 440.4080 441.1878 441.3884 441.3884 441.4833
## p = 9  425.5120 422.2572 424.2417 426.2314 427.4237 428.7796 428.7683 428.7683
## p = 10 426.8926 421.5270 423.5088 425.4984 427.2061 429.1897 425.7288 427.6503
## p = 11 423.3295 415.9815 417.9207 419.9021 421.4330 423.3641 420.3817 422.1188
## p = 12 420.3664 414.7203 416.7060 418.7031 420.3577 422.3076 419.0020 420.8891
## p = 13 419.2491 414.0452 416.0388 417.8335 419.6142 421.6133 417.3752 419.3688
## p = 14 413.9923 409.6999 411.6644 413.6117 415.4139 417.4134 414.9856 416.9121
## p = 15 410.8162 407.1743 409.1688 411.1347 412.9626 414.9604 411.7240 413.5707
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  430.7683 429.7617 426.1571 424.8081 423.8523 418.9567 417.1553
## p = 2  430.8664 429.9994 426.2612 424.8154 423.9391 418.9764 417.1931
## p = 3  432.8627 431.9976 428.2209 426.7665 425.8855 420.9632 419.1786
## p = 4  434.5332 433.6460 429.5783 428.2528 427.3860 422.5940 420.7018
## p = 5  436.5073 435.6241 431.5763 430.2482 429.3836 424.5870 422.6849
## p = 6  438.0351 436.9871 433.1083 431.6076 430.7787 425.2356 423.3752
## p = 7  439.6069 438.5454 435.0971 433.6019 432.7720 427.2275 425.3715
## p = 8  439.8903 438.9000 435.4771 434.4580 433.6425 427.6642 426.1569
## p = 9  429.7311 427.9020 422.7188 421.8916 419.4369 414.2641 410.7542
## p = 10 427.6503 429.0586 423.0306 422.2958 419.8671 411.8350 408.5194
## p = 11 421.0320 421.0320 422.5983 421.9689 419.9512 411.6340 409.4523
## p = 12 421.4903 421.3863 421.3863 423.3750 420.8079 411.4564 409.3368
## p = 13 421.1008 418.6826 420.0560 420.0560 420.8020 412.6027 410.3545
## p = 14 418.3704 417.3584 419.2114 412.9575 412.9575 414.6003 412.2951
## p = 15 414.7440 413.8658 415.5867 411.8044 412.4660 412.4660 412.5168
## 
## $min.Stat
## [1] 407.1743

##Pemodelan DLM dan ARDL dengan Library dynLm

#sama dengan model dlm p=1
cons_lm1 <- dynlm(meanpressure ~ wind_speed+L(wind_speed),data = train.ts)

#sama dengan model ardl p=0 q=1
cons_lm2 <- dynlm(meanpressure ~ wind_speed+L(meanpressure),data = train.ts)

#sama dengan ardl p=1 q=1
cons_lm3 <- dynlm(meanpressure ~ wind_speed+L(wind_speed)+L(meanpressure),data = train.ts)

#sama dengan dlm p=2
cons_lm4 <- dynlm(meanpressure ~ wind_speed+L(wind_speed)+L(wind_speed,2),data = train.ts)

#Ringkasan Model
summary(cons_lm1)
## 
## Time series regression with "ts" data:
## Start = 2, End = 114
## 
## Call:
## dynlm(formula = meanpressure ~ wind_speed + L(wind_speed), data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.1313  -3.7613   0.2305   4.0520  10.4525 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1015.6749     1.5174 669.350   <2e-16 ***
## wind_speed      -0.2049     0.1672  -1.225    0.223    
## L(wind_speed)   -0.1971     0.1664  -1.184    0.239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.627 on 110 degrees of freedom
##   (0 observations deleted due to missingness)
## Multiple R-squared:  0.04604,    Adjusted R-squared:  0.0287 
## F-statistic: 2.655 on 2 and 110 DF,  p-value: 0.07483
summary(cons_lm2)
## 
## Time series regression with "ts" data:
## Start = 2, End = 114
## 
## Call:
## dynlm(formula = meanpressure ~ wind_speed + L(meanpressure), 
##     data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5427  -3.7704   0.5197   4.4160  11.0536 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.016e+03  6.018e+00 168.787   <2e-16 ***
## wind_speed      -2.915e-01  1.513e-01  -1.927   0.0566 .  
## L(meanpressure) -9.225e-04  6.006e-03  -0.154   0.8782    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.663 on 110 degrees of freedom
##   (0 observations deleted due to missingness)
## Multiple R-squared:  0.03409,    Adjusted R-squared:  0.01652 
## F-statistic: 1.941 on 2 and 110 DF,  p-value: 0.1485
summary(cons_lm3)
## 
## Time series regression with "ts" data:
## Start = 2, End = 114
## 
## Call:
## dynlm(formula = meanpressure ~ wind_speed + L(wind_speed) + L(meanpressure), 
##     data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.1310  -3.7662   0.2414   4.0646  10.4549 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.016e+03  6.013e+00 168.963   <2e-16 ***
## wind_speed      -2.041e-01  1.685e-01  -1.211    0.229    
## L(wind_speed)   -1.963e-01  1.677e-01  -1.170    0.244    
## L(meanpressure) -3.542e-04  6.016e-03  -0.059    0.953    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.653 on 109 degrees of freedom
##   (0 observations deleted due to missingness)
## Multiple R-squared:  0.04607,    Adjusted R-squared:  0.01982 
## F-statistic: 1.755 on 3 and 109 DF,  p-value: 0.1601
summary(cons_lm4)
## 
## Time series regression with "ts" data:
## Start = 3, End = 114
## 
## Call:
## dynlm(formula = meanpressure ~ wind_speed + L(wind_speed) + L(wind_speed, 
##     2), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.178  -3.509   0.127   4.128  10.605 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1016.0760     1.7212 590.336   <2e-16 ***
## wind_speed         -0.1974     0.1686  -1.171    0.244    
## L(wind_speed)      -0.1251     0.1848  -0.677    0.500    
## L(wind_speed, 2)   -0.1339     0.1672  -0.801    0.425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.651 on 108 degrees of freedom
##   (0 observations deleted due to missingness)
## Multiple R-squared:  0.04651,    Adjusted R-squared:  0.02003 
## F-statistic: 1.756 on 3 and 108 DF,  p-value: 0.1599

#SSE

deviance(cons_lm1)
## [1] 3483.506
deviance(cons_lm2)
## [1] 3527.166
deviance(cons_lm3)
## [1] 3483.395
deviance(cons_lm4)
## [1] 3448.526

#Uji Diagnostik Model #UJi Non Autokorelasi #uji model

if(require("lmtest")) encomptest(cons_lm1, cons_lm2)
## Encompassing test
## 
## Model 1: meanpressure ~ wind_speed + L(wind_speed)
## Model 2: meanpressure ~ wind_speed + L(meanpressure)
## Model E: meanpressure ~ wind_speed + L(wind_speed) + L(meanpressure)
##           Res.Df Df      F Pr(>F)
## M1 vs. ME    109 -1 0.0035 0.9532
## M2 vs. ME    109 -1 1.3696 0.2444

#Diagnostik #Durbin Watson

dwtest(cons_lm1)
## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 0.12968, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(cons_lm2)
## 
##  Durbin-Watson test
## 
## data:  cons_lm2
## DW = 0.1446, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(cons_lm3)
## 
##  Durbin-Watson test
## 
## data:  cons_lm3
## DW = 0.12966, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(cons_lm4)
## 
##  Durbin-Watson test
## 
## data:  cons_lm4
## DW = 0.13477, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

#Uji Heterogenitas

bptest(cons_lm1)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 0.8666, df = 2, p-value = 0.6484
bptest(cons_lm2)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 0.78429, df = 2, p-value = 0.6756
bptest(cons_lm3)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm3
## BP = 1.0519, df = 3, p-value = 0.7887
bptest(cons_lm4)
## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm4
## BP = 2.0749, df = 3, p-value = 0.557

#Uji Normalitas #Shapiro Wilk

shapiro.test(residuals(cons_lm1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm1)
## W = 0.97648, p-value = 0.04351
shapiro.test(residuals(cons_lm2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.97922, p-value = 0.07572
shapiro.test(residuals(cons_lm3))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm3)
## W = 0.97664, p-value = 0.04502
shapiro.test(residuals(cons_lm4))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm4)
## W = 0.97573, p-value = 0.03888

##Perbandingan Ramalan Keakuratan

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            NA
## DLM 1            NA
## DLM 2            NA
## Autoregressive   NA