Regresi Linier dengan Peubah Lag pada Data Cuaca Harian

Library

library(dLagM)

## Loading required package: nardl

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Loading required package: dynlm

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

library(dynlm) 
library(MLmetrics)

## 
## Attaching package: 'MLmetrics'

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

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

library(lmtest)
library(car)

## Loading required package: carData

library(readxl)
library(readr)
library(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

Data ini merupakan data cuaca di kota Delhi dengan rentang harian selama periode 7 bulan (1 Juni 2016 - 31 Desember 2016). Data ini memiliki 4 variabel numerik, yaitu rata-rata suhu, kelembapan, tekanan udara, dan kecepatan udara.

Analisis ini menggunakan peubah rataan suhu digunakan sebagai peubah X dan peubah tekanan udara digunakan sebagai peubah Y. Dataset bersumber dari kaggle.com.

library(readr)
MPDW <- read_csv("C:/Users/62812/Documents/SEMESTER 5/MPDW/Data MPDW Climate.csv")

## Rows: 214 Columns: 5
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## dbl  (4): Suhu, Kelembapan, Kecepatan, Tekanan
## date (1): Date
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(MPDW)

## # A tibble: 6 x 5
##   Date        Suhu Kelembapan Kecepatan Tekanan
##   <date>     <dbl>      <dbl>     <dbl>   <dbl>
## 1 2016-06-01  36         42.2      9.13   1003.
## 2 2016-06-02  37.6       35.3     11.9    1002 
## 3 2016-06-03  37.6       40.7      2.78   1002.
## 4 2016-06-04  38.2       42.1      3.21   1002.
## 5 2016-06-05  36.2       51.8      8.5    1003.
## 6 2016-06-06  35.4       45.7     11      1002.

str(MPDW)

## spec_tbl_df [214 x 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ Date      : Date[1:214], format: "2016-06-01" "2016-06-02" ...
##  $ Suhu      : num [1:214] 36 37.6 37.6 38.2 36.2 ...
##  $ Kelembapan: num [1:214] 42.2 35.3 40.7 42.1 51.8 ...
##  $ Kecepatan : num [1:214] 9.13 11.93 2.78 3.21 8.5 ...
##  $ Tekanan   : num [1:214] 1003 1002 1002 1002 1003 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   Date = col_date(format = ""),
##   ..   Suhu = col_double(),
##   ..   Kelembapan = col_double(),
##   ..   Kecepatan = col_double(),
##   ..   Tekanan = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

dim(MPDW)

## [1] 214   5

MPDW$Date <- as.Date(MPDW$Date,"%y/%m/%d")

t<-MPDW$Date
Xt<-MPDW$Suhu
Yt<-MPDW$Kecepatan

datareg1<-cbind(t, Xt, Yt)
datareg <- as.data.frame(datareg1)
View(datareg)

Model Regresi Awal

model1 <- lm(Yt ~ Xt, data = datareg)
summary(model1)

## 
## Call:
## lm(formula = Yt ~ Xt, data = datareg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8145 -2.9245 -0.2871  1.8869 13.7982 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  3.06697    1.26535   2.424  0.01620 * 
## Xt           0.14070    0.04387   3.207  0.00155 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.692 on 212 degrees of freedom
## Multiple R-squared:  0.04628,    Adjusted R-squared:  0.04178 
## F-statistic: 10.29 on 1 and 212 DF,  p-value: 0.001547

Berdasarkan output di atas, didapatkan model regresi deret waktu awal dengan persamaan: Yt = 0.873030 + 0.016550(Xt)

Pengecekan Autokorelasi Model Awal

bgtest(model1)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model1
## LM test = 49.646, df = 1, p-value = 1.841e-12

dwtest(model1)

## 
##  Durbin-Watson test
## 
## data:  model1
## DW = 1.0356, p-value = 4.183e-13
## alternative hypothesis: true autocorrelation is greater than 0

Setelah dilakukan diagnostik autokorelasi pada model regresi awal, dapat dilihat bahwa kedua uji menunjukkan p−value < 0.05 (Tolak H0). Artinya cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada model regresi awal. Oleh karena itu, akan dilakukan penanganan dengan menambahkan peubah lag melalui beberapa model, di antaranya model KOYCK, model autoregressive, dan model distributed lag.

Splitting Data

Splitting data bertujuan untuk membagi data menjadi data training dan data testing. Jumlah data training yaitu 172 observasi (80%), sementara data testing berjumlah 43 observasi (20%). Data training digunakan sebagai data sebelum inisiasi peramalan, sedangkan dan data testing digunakan sebagai data peramalannya.

train<-datareg[1:171, ]
test<-datareg[172:214, ]

#data time series
train.ts<-ts(train)
test.ts<-ts(test)
data.ts<-ts(datareg)

Model KOYCK

Model Koyck merupakan suatu model yang mengasumsikan bahwa koefisien dari variabel yang mengalami keterlambatan dapat menurun secara geometris.

Metode Koyck didasarkan asumsi bahwa semakin jauh jarak lag variabel bebas dari periode sekarang maka semakin kecil pengaruh variabel lag terhadap variabel tak bebas. Koyck mengusulkan suatu metode pendugaan model lag terdistribusi didasarkan pada asumsi bahwa koefisien menurun secara eksponensial dari waktu ke waktu.

#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 
## -9.9620 -2.2149 -0.1777  1.5050 12.4002 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.49780    2.39432  -0.208   0.8356    
## Y.1          0.42910    0.06971   6.156 5.38e-09 ***
## X.t          0.15020    0.08058   1.864   0.0641 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.199 on 167 degrees of freedom
## Multiple R-Squared: 0.2435,  Adjusted R-squared: 0.2345 
## Wald test: 25.47 on 2 and 167 DF,  p-value: 2.223e-10 
## 
## Diagnostic tests:
## NULL
## 
##                               alpha     beta       phi
## Geometric coefficients:  -0.8719619 0.150202 0.4291006

Berdasarkan output di atas, didapatkan nilai instersep model (b0) sebesar -0.49780, nilai X.t sebesar 0.15020 dan nilai Yt-1 sebesar 0.4291. Dengan demikian, dapat disimpulkan bahwa pemodelan koyck yang dibentuk dari data adalah sebagai berikut: Yt = -0.4978 + 0.1502(Xt) + 0.4291(Yt-1)

AIC dan BIC

AIC dan BIC adalah suatu metrik yang digunakan untuk mengukur kebaikan model. Semakin kecil nilai AIC dan BIC yang diperoleh, maka semakin bagus pula suatu model.

AIC(model.koyck)

## [1] 882.7284

BIC(model.koyck)

## [1] 895.2716

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

## $forecasts
##  [1] 5.944775 5.395099 5.052358 5.061570 5.102644 4.909230 5.153519 5.362601
##  [9] 5.353522 5.172115 4.967587 4.618241 4.450377 4.296987 4.531571 4.285515
## [17] 4.140334 4.063331 3.986962 3.751419 4.028356 3.700756 4.100493 4.252107
## [25] 4.188127 4.086407 3.984347 3.590082 3.671239 3.493277 3.683344 4.056722
## [33] 3.954086 3.795247 3.466381 3.590849 3.145861 3.426980 3.503620 3.591690
## [41] 3.332186 3.049174 3.071535
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 43)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.6798613
## 
## $MAPE_training.MAPE
## [1] 0.4941523

Berdasarkan output di atas, terlihat bahwa nilai MAPE Testing lebih besar daripada nilai MAPE Training. Namun, nilai MAPE Testing dan MAPE Training memiliki perbedaan tidak cukup jauh. Hal ini mengindikasikan bahwa model ini tidak underfitted atau overfitted.

Regression with Distributed Lag

Distributed Lag adalah model regresi yang memuat variabel takbebas yang dipengaruhi oleh variabel bebas pada waktu t, serta dipengaruhi juga oleh variabel bebas. Disebut distributed lag sebab pengaruh dari satu atau beberapa variabel bebas (X) terhadap variabel takbebas (Y) menyebar ke beberapa periode waktu. ### Lag = 2

model.dlm = dLagM::dlm(x = train$Xt,y = train$Yt , q = 2)
summary(model.dlm)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6678 -2.6921 -0.1816  2.0256 13.3907 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -0.4318     2.5355  -0.170   0.8650  
## x.t           0.3063     0.1811   1.691   0.0927 .
## x.1           0.2713     0.2254   1.204   0.2303  
## x.2          -0.3283     0.1826  -1.798   0.0740 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.514 on 165 degrees of freedom
## Multiple R-squared:  0.08882,    Adjusted R-squared:  0.07225 
## F-statistic: 5.361 on 3 and 165 DF,  p-value: 0.00151
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 910.3184 925.9679

AIC(model.dlm)

## [1] 910.3184

BIC(model.dlm)

## [1] 925.9679

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

## $forecasts
##  [1] 4.661822 4.858106 5.419047 5.024974 5.616557 4.911681 5.116784 6.381825
##  [9] 5.653431 4.885103 4.521989 4.147605 3.915307 4.288609 4.793527 4.806900
## [17] 3.443319 4.099514 4.071128 3.611465 4.110763 4.326299 3.795646 5.706708
## [25] 4.226782 3.885745 3.914462 3.256501 3.261586 4.045772 3.657654 5.199101
## [33] 4.607938 3.261247 3.096632 3.418190 3.451277 2.933267 4.785510 3.786938
## [41] 3.376455 2.368736 3.001623
## 
## $call
## forecast.dlm(model = model.dlm, x = test$Xt, h = 43)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

MAPE Testing dan Training

#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.6579371
## 
## $MAPE_training.MAPE
## [1] 0.5914658

Terlihat bahwa nilai MAPE Testing lebih besar dibandingkan nilai MAPE Training. MAPE Testing dan MAPE Training memiliki perbedaan nilai yang juga tidak jauh berbeda, yaitu hanya sebesar 0.067. Hal ini mengindikasikan bahwa model ini tidak underfitted atau overfitted.

Lag Optimum

Lag optimum dapat ditentukan dengan melihat masing-masing nilai AIC, BIC, serta MAPE yang optimum di antara model-model tersebut.

#penentuan lag optimum 
finiteDLMauto(formula = Yt ~ Xt,
              data = data.frame(train), q.min = 1, q.max = 4 ,
              model.type = "dlm", error.type = "AIC", trace = TRUE) ##q max lag maksimum

##   q - k    MASE      AIC      BIC   GMRAE   MBRAE R.Adj.Sq    Ljung-Box
## 4     4 1.01846 896.7649 918.5909 1.04896 0.46139  0.09501 1.027108e-08
## 3     3 1.02402 902.9247 921.6685 1.05299 0.35948  0.08632 2.069213e-08
## 2     2 1.04759 910.3184 925.9679 1.09395 0.52686  0.07225 7.995897e-09
## 1     1 1.03271 916.5817 929.1249 1.02055 0.11469  0.06578 1.492963e-08

#model dlm dengan lag optimum
model.dlm2 = dLagM::dlm(x = train$Xt,y = train$Yt , q = 4) #terdapat lag yang tidak signifikan sehingga dapat dikurangi jumlah lagnya 
summary(model.dlm2)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1903 -2.5602 -0.4377  1.8816 13.7363 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.95308    2.66309  -0.358   0.7209  
## x.t          0.38009    0.18135   2.096   0.0377 *
## x.1          0.31163    0.22309   1.397   0.1644  
## x.2         -0.17517    0.22311  -0.785   0.4335  
## x.3         -0.15563    0.22293  -0.698   0.4861  
## x.4         -0.09047    0.18236  -0.496   0.6205  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.464 on 161 degrees of freedom
## Multiple R-squared:  0.1223, Adjusted R-squared:  0.09501 
## F-statistic: 4.486 on 5 and 161 DF,  p-value: 0.0007457
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 896.7649 918.5909

Dapat terlihat berdasarkan output summary model distributed lag, hanya terdapat 1 peubah yang signifikan terhadap model, yaitu peubah Xt.

AIC(model.dlm2)

## [1] 896.7649

BIC(model.dlm2)

## [1] 918.5909

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

## $forecasts
##  [1] 4.467751 4.786537 5.199376 5.261377 5.654064 4.982262 5.231773 6.288149
##  [9] 6.069073 5.081326 4.202507 3.524292 3.327295 3.627392 4.587061 4.603991
## [17] 3.529218 3.551945 3.636262 3.300260 3.925745 3.870516 4.106920 5.428650
## [25] 4.665552 3.942207 3.485971 2.710194 2.839337 3.334231 3.673335 5.164459
## [33] 4.863720 3.539662 2.485451 2.847068 2.718175 2.908250 4.239444 4.093724
## [41] 3.318291 1.963964 2.287738
## 
## $call
## forecast.dlm(model = model.dlm2, x = test$Xt, h = 43)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.6332332
## 
## $MAPE_training.MAPE
## [1] 0.5647357

Model Autoregressive

Model autoregressive (AR) adalah model regresi time series yang menghubungkan nilai pengamatan aktual dengan nilai pengamatan sebelumnya. Autoregressive sendiri adalah variabel takbebas dipengaruhi oleh variabel bebas pada waktu t, serta dipengaruhi juga oleh variabel takbebas itu sendiri.

Konsep dasar model ini yaitu meregresikan pengamatan aktual dengan nilai pengamatan sebelumnya untuk melakukan peramalan nilai ke depan.

#MODEL AUTOREGRESSIVE 
#library(dLagM)
model.ardl = ardlDlm(x = train$Xt, y = train$Yt, p = 1 , q = 1) #p:lag x, q:lag y
#model untuk p=1, q=1: yt=b0+b1yt-1+b2xt+b3xt-1
#model untuk p=2, q=3: yt=b0+b1yt-1+b2yt-2+b3xt+b4xt-1+b5xt-2

summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 2, End = 171
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0624  -2.1056  -0.2035   1.6505  12.3906 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.93773    2.23214  -0.420   0.6750    
## X.t          0.32240    0.16198   1.990   0.0482 *  
## X.1         -0.15830    0.16631  -0.952   0.3425    
## Y.1          0.43305    0.07029   6.161 5.28e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.197 on 166 degrees of freedom
## Multiple R-squared:  0.2486, Adjusted R-squared:  0.2351 
## F-statistic: 18.31 on 3 and 166 DF,  p-value: 2.598e-10

AIC(model.ardl)

## [1] 883.575

BIC(model.ardl)

## [1] 899.2539

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

## $forecasts
##  [1] 5.654258 5.412615 4.827922 5.022810 5.022170 4.529788 5.241474 5.428518
##  [9] 5.187576 4.806344 4.556415 4.020230 4.025180 3.871619 4.535665 3.762419
## [17] 3.707986 3.694585 3.611281 3.185631 4.026385 3.033834 4.234264 4.141915
## [25] 3.845940 3.694357 3.581609 2.842082 3.428537 2.961930 3.556003 4.158842
## [33] 3.548141 3.314304 2.774437 3.385478 2.300401 3.369043 3.239876 3.348747
## [41] 2.699611 2.363376 2.707351
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 43)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.6606996
## 
## $MAPE_training.MAPE
## [1] 0.4943607

ardlBoundOrders(data = data.frame(datareg) , formula = Yt ~ Xt )

## $p
##   Xt
## 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  1101.215 1091.007 1088.315 1081.051 1077.106 1057.279 1050.946 1048.127
## p = 2  1092.054 1092.016 1089.151 1081.596 1077.454 1054.873 1049.317 1046.440
## p = 3  1089.011 1089.011 1090.828 1083.326 1079.180 1056.444 1051.219 1048.295
## p = 4  1086.541 1087.282 1087.282 1084.986 1080.691 1057.781 1052.385 1049.731
## p = 5  1078.408 1079.058 1080.076 1080.076 1081.811 1059.455 1054.020 1051.298
## p = 6  1058.745 1057.211 1059.180 1060.068 1060.068 1061.270 1055.698 1052.988
## p = 7  1055.678 1053.665 1055.658 1055.350 1055.512 1055.512 1057.512 1054.742
## p = 8  1051.657 1049.578 1051.560 1051.311 1052.499 1053.020 1053.020 1054.944
## p = 9  1049.036 1047.372 1049.367 1049.234 1050.537 1050.970 1052.949 1052.949
## p = 10 1046.493 1044.320 1046.320 1045.711 1046.866 1047.721 1049.582 1051.215
## p = 11 1039.810 1036.968 1038.968 1037.704 1038.594 1039.254 1040.933 1042.401
## p = 12 1036.130 1033.077 1035.061 1034.212 1035.725 1036.841 1038.644 1040.064
## p = 13 1034.158 1031.299 1033.294 1032.600 1034.370 1034.977 1036.665 1038.314
## p = 14 1027.469 1024.705 1026.584 1025.751 1027.683 1028.403 1030.307 1031.628
## p = 15 1024.599 1022.246 1023.987 1023.249 1025.140 1025.909 1027.891 1029.349
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  1045.257 1040.135 1037.253 1034.431 1029.805 1027.381 1025.041
## p = 2  1043.320 1038.180 1035.089 1031.907 1027.620 1025.150 1022.503
## p = 3  1045.230 1039.962 1036.846 1033.603 1029.067 1026.635 1023.939
## p = 4  1046.512 1041.066 1038.111 1034.879 1030.574 1028.235 1025.498
## p = 5  1047.924 1042.648 1039.752 1036.621 1032.273 1029.874 1027.025
## p = 6  1049.757 1044.290 1041.450 1038.282 1033.948 1031.522 1028.584
## p = 7  1051.405 1045.758 1043.091 1040.006 1035.668 1033.288 1030.371
## p = 8  1051.696 1046.504 1043.957 1040.372 1035.967 1033.630 1030.957
## p = 9  1053.571 1048.350 1045.749 1042.200 1037.957 1035.623 1032.956
## p = 10 1051.215 1050.280 1047.696 1044.169 1039.925 1037.525 1034.815
## p = 11 1044.394 1044.394 1045.321 1041.366 1036.746 1034.279 1032.216
## p = 12 1041.754 1041.479 1041.479 1043.302 1038.718 1036.269 1034.199
## p = 13 1040.073 1038.559 1039.987 1039.987 1040.299 1037.815 1035.747
## p = 14 1033.371 1031.840 1032.819 1034.800 1034.800 1036.618 1034.639
## p = 15 1031.133 1028.833 1029.941 1031.791 1033.401 1033.401 1035.355
## 
## $min.Stat
## [1] 1022.246

#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)
summary(cons_lm1)

## 
## Time series regression with "ts" data:
## Start = 2, End = 171
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.428 -2.332 -0.270  2.020 12.935 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.71679    2.46283  -0.697    0.487
## Xt           0.26960    0.17876   1.508    0.133
## L(Xt)        0.02126    0.18095   0.117    0.907
## 
## Residual standard error: 3.533 on 167 degrees of freedom
## Multiple R-squared:  0.07684,    Adjusted R-squared:  0.06578 
## F-statistic:  6.95 on 2 and 167 DF,  p-value: 0.001261

summary(cons_lm1)

## 
## Time series regression with "ts" data:
## Start = 2, End = 171
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.428 -2.332 -0.270  2.020 12.935 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.71679    2.46283  -0.697    0.487
## Xt           0.26960    0.17876   1.508    0.133
## L(Xt)        0.02126    0.18095   0.117    0.907
## 
## Residual standard error: 3.533 on 167 degrees of freedom
## Multiple R-squared:  0.07684,    Adjusted R-squared:  0.06578 
## F-statistic:  6.95 on 2 and 167 DF,  p-value: 0.001261

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 171
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.9400 -2.1866 -0.2348  1.5070 12.3084 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.48588    2.15597  -0.689   0.4917    
## Xt           0.18434    0.07210   2.557   0.0115 *  
## L(Yt)        0.42133    0.06918   6.090 7.52e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.196 on 167 degrees of freedom
## Multiple R-squared:  0.2445, Adjusted R-squared:  0.2355 
## F-statistic: 27.03 on 2 and 167 DF,  p-value: 6.77e-11

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 171
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0624  -2.1056  -0.2035   1.6505  12.3906 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.93773    2.23214  -0.420   0.6750    
## Xt           0.32240    0.16198   1.990   0.0482 *  
## L(Xt)       -0.15830    0.16631  -0.952   0.3425    
## L(Yt)        0.43305    0.07029   6.161 5.28e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.197 on 166 degrees of freedom
## Multiple R-squared:  0.2486, Adjusted R-squared:  0.2351 
## F-statistic: 18.31 on 3 and 166 DF,  p-value: 2.598e-10

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 171
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6678 -2.6921 -0.1816  2.0256 13.3907 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -0.4318     2.5355  -0.170   0.8650  
## Xt            0.3063     0.1811   1.691   0.0927 .
## L(Xt)         0.2713     0.2254   1.204   0.2303  
## L(Xt, 2)     -0.3283     0.1826  -1.798   0.0740 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.514 on 165 degrees of freedom
## Multiple R-squared:  0.08882,    Adjusted R-squared:  0.07225 
## F-statistic: 5.361 on 3 and 165 DF,  p-value: 0.00151

SSE

deviance(cons_lm1)

## [1] 2084.99

deviance(cons_lm2)

## [1] 1706.229

deviance(cons_lm3)

## [1] 1696.967

deviance(cons_lm4)

## [1] 2037.237

Diagnostik Model

Uji Non Autokorelasi

if(require("lmtest")) encomptest(cons_lm1, cons_lm2)

## Encompassing test
## 
## Model 1: Yt ~ Xt + L(Xt)
## Model 2: Yt ~ Xt + L(Yt)
## Model E: Yt ~ Xt + L(Xt) + L(Yt)
##           Res.Df Df       F    Pr(>F)    
## M1 vs. ME    166 -1 37.9571 5.284e-09 ***
## M2 vs. ME    166 -1  0.9061    0.3425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#durbin watson
dwtest(cons_lm1)

## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 1.1306, p-value = 4.301e-09
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm2)

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

dwtest(cons_lm3)

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

dwtest(cons_lm4)

## 
##  Durbin-Watson test
## 
## data:  cons_lm4
## DW = 1.0947, p-value = 1.303e-09
## alternative hypothesis: true autocorrelation is greater than 0

Uji Heterogenitas

bptest(cons_lm1)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 1.916, df = 2, p-value = 0.3837

bptest(cons_lm2)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm2
## BP = 2.5692, df = 2, p-value = 0.2768

bptest(cons_lm3)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm3
## BP = 3.7691, df = 3, p-value = 0.2875

bptest(cons_lm4)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm4
## BP = 3.6115, df = 3, p-value = 0.3066

Normalitas

shapiro.test(residuals(cons_lm1))

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

shapiro.test(residuals(cons_lm2))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.96035, p-value = 9.398e-05

shapiro.test(residuals(cons_lm3))

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

shapiro.test(residuals(cons_lm4))

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

Perbandingan Keakuratan Ramalan

#PERBANDINGAN
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.6798613
## DLM 1          0.6579371
## DLM 2          0.6332332
## Autoregressive 0.6606996

Berdasarkan perbandingan keempat model di atas, diperoleh bahwa DLM 2 (Distributed Lag Model) memiliki nilai keakuratan MAPE terkecil. Oleh karena itu, dapat dikatakan bahwa model peramalan terbaik ialah DLM 2.

#PLOT
par(mfrow=c(1,1))
plot(test$Xt, test$Yt, type="b", col="black")
points(test$Xt, fore.koyck$forecasts,col="red")
lines(test$Xt, fore.koyck$forecasts,col="red")
points(test$Xt, fore.dlm$forecasts,col="blue")
lines(test$Xt, fore.dlm$forecasts,col="blue")
points(test$Xt, fore.dlm2$forecasts,col="orange")
lines(test$Xt, fore.dlm2$forecasts,col="orange")
points(test$Xt, fore.ardl$forecasts,col="green")
lines(test$Xt, fore.ardl$forecasts,col="green")
legend("topleft",c("aktual", "koyck","DLM 1","DLM 2", "autoregressive"), lty=1, col=c("black","red","blue","orange","green"), cex=0.8)

Kesimpulan

Metode yang paling cocok untuk metode peramalan terbaik pada data ini, yaitu Distributed Lag Model. Kemudian, hasil uji diagnostik menunjukkan bahwa dengan Distributed Lag Model, autokorelasi pada model regresi deret waktu berhasil ditangani dengan menggunakan taraf nyata 5%.

Daftar Pustaka

https://www.kaggle.com/datasets/sumanthvrao/daily-climate-time-series-data

Aqibah M, Suciptawati NLP, Sumarjaya IW. 2020. Model dinamis Autoregressive Distributed Lag. (Studi kasus: pengaruh kurs Dolar Amerika dan inflasi terhadap harga saham tahun 2014-2018. E-Jurnal Matematika. 9(4):240-250. Gujarati, D. N. (2004). Basic Econometrics. Fourth Edition. New York: The McGrawHill.

Gujarati DN. 2004. Basic Econometrics. Fourth Edition. New York: The McGrawHill.

Gujarati DN. 2006. Essentials of Econometrics. Third Edition. New York: The McGraw-Hill.