Tugas Individu 2 STA1341 MPDW

Library

library(dLagM) 
library(dynlm)
library(MLmetrics)
library(lmtest)
library(car)
library(readxl)

Input Data

Data di bawah ini merupakan data yang bersumber dari kaggle.com, membahas mengenai cuaca di Kota New Delhi, India dari tanggal 1 Januari 2013 hingga 24 April 2017.

train.data <- read.csv("D:/IPB/Semester 5/STA1341 - Metode Peramalan Deret Waktu/4. Peubah Lag/DailyDelhiClimateTrain.csv")
test.data <- read.csv("D:/IPB/Semester 5/STA1341 - Metode Peramalan Deret Waktu/4. Peubah Lag/DailyDelhiClimateTest.csv")

str(c(train.data, test.data))

## List of 10
##  $ date        : chr [1:1462] "2013-01-01" "2013-01-02" "2013-01-03" "2013-01-04" ...
##  $ meantemp    : num [1:1462] 10 7.4 7.17 8.67 6 ...
##  $ humidity    : num [1:1462] 84.5 92 87 71.3 86.8 ...
##  $ wind_speed  : num [1:1462] 0 2.98 4.63 1.23 3.7 ...
##  $ meanpressure: num [1:1462] 1016 1018 1019 1017 1016 ...
##  $ date        : chr [1:114] "2017-01-01" "2017-01-02" "2017-01-03" "2017-01-04" ...
##  $ meantemp    : num [1:114] 15.9 18.5 17.1 18.7 18.4 ...
##  $ humidity    : num [1:114] 85.9 77.2 81.9 70 74.9 ...
##  $ wind_speed  : num [1:114] 2.74 2.89 4.02 4.54 3.3 ...
##  $ meanpressure: num [1:114] 59 1018 1018 1016 1014 ...

head(train.data)

##         date  meantemp humidity wind_speed meanpressure
## 1 2013-01-01 10.000000 84.50000   0.000000     1015.667
## 2 2013-01-02  7.400000 92.00000   2.980000     1017.800
## 3 2013-01-03  7.166667 87.00000   4.633333     1018.667
## 4 2013-01-04  8.666667 71.33333   1.233333     1017.167
## 5 2013-01-05  6.000000 86.83333   3.700000     1016.500
## 6 2013-01-06  7.000000 82.80000   1.480000     1018.000

head(test.data)

##         date meantemp humidity wind_speed meanpressure
## 1 2017-01-01 15.91304 85.86957   2.743478       59.000
## 2 2017-01-02 18.50000 77.22222   2.894444     1018.278
## 3 2017-01-03 17.11111 81.88889   4.016667     1018.333
## 4 2017-01-04 18.70000 70.05000   4.545000     1015.700
## 5 2017-01-05 18.38889 74.94444   3.300000     1014.333
## 6 2017-01-06 19.31818 79.31818   8.681818     1011.773

Peubah yang dipilih untuk dianalisis lebih lanjut adalah rata-rata temperatur (meantemp) dan kelembapan (humidity)

train.data <- train.data[c(1,2,3)]
test.data <- test.data[c(1,2,3)]
colnames(train.data) <- c("t", "Xt", "Yt")
colnames(test.data) <- c("t", "Xt", "Yt")

Data diubah dalam bentuk data time series

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

Model KOYCK

Metode Koyck dapat digunakan untuk menentukan estimasi model dinamis terdistribusi lag yang panjang beda kala (lag) tidak diketahui. Pada persamaan Koyck diakhiri dengan model autoregresif karena muncul variabel bebas Yt (Pratami et al. 2016).

model.koyck = dLagM::koyckDlm(x = train.data$Xt, y = train.data$Yt)
summary(model.koyck)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -35.2141  -4.7702  -0.4399   4.1195  44.4010 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.73159    1.59267   6.738  2.3e-11 ***
## Y.1          0.85937    0.01486  57.828  < 2e-16 ***
## X.t         -0.08542    0.03480  -2.454   0.0142 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.902 on 1458 degrees of freedom
## Multiple R-Squared: 0.7781,  Adjusted R-squared: 0.7778 
## Wald test:  2539 on 2 and 1458 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha        beta       phi
## Geometric coefficients:  76.31293 -0.08541751 0.8593739

Pada pemodelan koyck didapatkan hasil yang signifikan untuk peubah Xt dan Y(t-1) terhadap nilai Y (Kelembapan). Persamaan model koyck berikut adalah:

\[Y_t = 10.73519 - 0.08542X_t + 0.85937Y_{t-1}\]

AIC(model.koyck)

## [1] 10191.43

BIC(model.koyck)

## [1] 10212.58

Peramalan menggunakan Metode Koyck

#ramalan 114 periode ke depan
(fore.koyck <- forecast(model = model.koyck, x=test.data$Xt, h=114))

## $forecasts
##   [1] 95.30973 91.05806 87.52292 84.34920 81.64835 79.24794 77.57886 76.06113
##   [9] 74.85188 74.02284 73.40530 72.80716 72.16964 71.62478 70.88024 70.39285
##  [17] 70.22409 69.96635 69.61012 69.24897 68.92766 68.39113 67.95758 67.56533
##  [25] 67.00170 66.92917 66.83939 66.90200 66.88722 66.80807 66.76735 66.80709
##  [33] 66.68400 66.70246 66.45663 66.25122 66.22356 66.35045 66.41038 66.48964
##  [41] 66.61826 66.64694 66.61820 66.58638 66.51276 66.39002 66.05573 65.67877
##  [49] 65.36977 64.99842 64.59290 64.37609 64.42111 64.50251 64.52976 64.56385
##  [57] 64.61450 64.56193 64.22135 63.83230 63.55874 63.60125 63.75524 63.83418
##  [65] 63.88067 63.69640 63.63025 63.63831 63.71768 63.90583 64.16646 64.38451
##  [73] 64.52445 64.48454 64.09790 64.03043 63.64851 63.32250 63.15622 62.87096
##  [81] 62.43362 61.99373 61.53739 61.35165 60.96995 60.60776 60.26328 59.87234
##  [89] 59.68280 59.40552 59.10317 58.98214 58.81407 58.63262 58.62211 58.44298
##  [97] 58.64969 58.94478 59.07025 59.11553 59.03547 58.96667 58.87196 58.71938
## [105] 58.52657 58.37985 58.12090 57.77498 57.52041 57.21623 56.97617 56.88519
## [113] 56.80914 56.81852
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test.data$Xt, h = 114)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

MAPE testing dan nilai akurasi data training

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

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

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

## $MAPE_testing
## [1] 0.3675654
## 
## $MAPE_training.MAPE
## [1] 0.1075979

Output menunjukkan bahwa nilai MAPE testing lebih besar dari nilai MAPE training. Model terindikasi tidak underfitted atau overfitted karena nilai dari kedua MAPE tersebut memiliki selisih yang kecil.

Regresi with Distributed Lag

Model distributed lag adalah suatu bentuk dari model regresi yang memperlihatkan hubungan antara peubah tak bebas dan peubah bebas dalam periode tertentu. Model dapat digunakan untuk melihat seberapa besar dampak yang diberikan oleh peubah bebas terhadap peubah tak bebas dari waktu ke waktu dan meramalkannya juga (Virgantari dan Rahayu 2021).

Regression with Distributed Lag `(Lag=2)`

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

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.730  -8.569  -0.331   9.204  33.463 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  92.6051     1.2938  71.574  < 2e-16 ***
## x.t          -2.7561     0.2155 -12.788  < 2e-16 ***
## x.1           0.5628     0.2783   2.023   0.0433 *  
## x.2           0.9445     0.2159   4.376  1.3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.51 on 1456 degrees of freedom
## Multiple R-squared:  0.3511, Adjusted R-squared:  0.3498 
## F-statistic: 262.6 on 3 and 1456 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 11750.52 11776.96

AIC(model.dlm)

## [1] 11750.52

BIC(model.dlm)

## [1] 11776.96

Pada pemodelan regresi with distributed lag didapatkan hasil yang signifikan untuk peubah X.t, X.1, dan X.2 terhadap nilai Y (Kelembapan) pada taraf nyata 5%.

#ramalan
(fore.dlm <- forecast(model = model.dlm, x=test.data$Xt, h=114)) #meramalkan 114 periode ke depan

## $forecasts
##   [1] 68.59372 60.01954 70.88848 68.17116 68.61102 67.37549 80.30971 75.90333
##   [9] 75.16532 82.24156 82.86807 77.74301 73.15309 74.80979 67.24028 73.94654
##  [17] 84.06653 77.09643 70.72630 71.08515 72.58438 64.86248 67.58783 69.65161
##  [25] 62.16713 77.16543 76.07099 76.20782 73.37722 70.13998 72.21664 75.18304
##  [33] 69.31517 73.53364 66.01020 66.54921 74.20737 78.02331 73.79978 73.35307
##  [41] 75.66547 72.31860 70.46639 71.49366 70.63803 69.10220 62.62347 61.89498
##  [49] 65.49826 62.97776 60.72222 66.71031 74.26223 72.65979 68.44223 68.59588
##  [57] 69.68681 66.32747 57.19312 57.09276 62.95530 72.57767 73.86600 68.29983
##  [65] 66.66385 60.16022 64.84907 68.81477 69.70523 72.38021 73.98118 71.81157
##  [73] 69.22542 64.43607 54.64734 67.28381 58.98340 58.26042 65.50271 60.13461
##  [81] 53.61213 54.38567 54.40659 62.29141 54.63130 52.70960 54.31871 51.85311
##  [89] 57.57048 53.96395 51.03149 57.19722 54.92397 52.61693 58.26160 52.22408
##  [97] 63.37351 66.66771 57.77170 55.29937 53.28943 54.76491 54.93215 52.86994
## [105] 51.74202 53.54808 49.81796 46.50259 50.16287 48.48971 49.26929 53.79381
## [113] 52.86676 53.98912
## 
## $call
## forecast.dlm(model = model.dlm, x = test.data$Xt, h = 114)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.3313371
## 
## $MAPE_training.MAPE
## [1] 0.2190723

Regression with Distributed Lag `Optimum`

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

##   q - k    MASE      AIC      BIC   GMRAE   MBRAE R.Adj.Sq Ljung-Box
## 4     4 1.76196 11710.00 11746.99 2.19782 0.78909  0.36108         0
## 3     3 1.76152 11729.37 11761.08 2.15908 1.07937  0.35538         0
## 2     2 1.76723 11750.52 11776.96 2.16998 0.56111  0.34977         0
## 1     1 1.77625 11774.77 11795.92 2.14709 1.84338  0.34319         0

Berdasarkan output di atas, lag optimum berada di q = 4karena beberapa parameter menunjukkan nilai AIC dan BIC yang paling minimum serta nilai Adj R-square yang paling besar dibandingkan dengan q lainnya.

#model dlm dengan lag optimum
model.dlm2 = dLagM::dlm(x = train.data$Xt,y = train.data$Yt , q = 4)
summary(model.dlm2)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.336  -8.704  -0.278   9.291  32.775 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  91.9328     1.2948  71.003  < 2e-16 ***
## x.t          -2.9559     0.2170 -13.618  < 2e-16 ***
## x.1           0.5185     0.2764   1.876 0.060874 .  
## x.2           0.1833     0.2771   0.661 0.508490    
## x.3           0.2261     0.2768   0.817 0.414293    
## x.4           0.8063     0.2174   3.709 0.000216 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.38 on 1452 degrees of freedom
## Multiple R-squared:  0.3633, Adjusted R-squared:  0.3611 
## F-statistic: 165.7 on 5 and 1452 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##     AIC      BIC
## 1 11710 11746.99

AIC(model.dlm2)

## [1] 11710

BIC(model.dlm2)

## [1] 11746.99

Pada pemodelan regresi with distributed lag optimum didapatkan hasil yang signifikan untuk peubah X.t dan X.4 terhadap nilai Y (Kelembapan) pada taraf nyata 5%.

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

## $forecasts
##   [1] 68.31144 62.09964 68.25973 60.58003 67.42120 66.57599 79.86634 75.97260
##   [9] 78.88227 85.46372 83.77229 80.94698 75.42567 74.18950 64.15715 72.06686
##  [17] 81.54616 76.50369 74.19335 71.24827 69.43023 62.02353 65.99798 66.25745
##  [25] 59.27070 77.32477 74.15294 79.05804 76.94328 70.94607 72.33000 73.75549
##  [33] 68.62503 74.27353 64.34562 65.58046 72.41306 76.23058 75.75164 76.23808
##  [41] 76.46489 71.91317 70.82296 70.65417 68.85634 68.00786 61.07344 59.67586
##  [49] 61.82259 59.43002 59.44770 65.56308 73.09348 74.11192 72.34532 71.02721
##  [57] 69.49635 65.72835 56.41830 54.92547 58.30018 69.43342 74.72560 72.31843
##  [65] 70.34318 59.89699 63.52754 66.25498 68.96539 74.40932 75.90757 74.04009
##  [73] 71.22037 64.58490 52.54767 64.37489 53.18566 57.10386 64.35178 57.08714
##  [81] 54.10802 53.04109 50.43448 60.04298 53.25627 53.79105 53.51819 49.14809
##  [89] 57.04940 52.81581 51.54387 57.48982 53.48838 53.54793 59.17999 51.34950
##  [97] 65.25308 67.45767 60.61461 59.62744 52.44978 52.53662 52.93515 51.59312
## [105] 51.08976 52.36800 48.31142 45.64279 48.68681 45.90886 48.65540 53.75105
## [113] 52.90208 55.97786
## 
## $call
## forecast.dlm(model = model.dlm2, x = test.data$Xt, h = 114)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

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

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

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

## $MAPE_testing
## [1] 0.331841
## 
## $MAPE_training.MAPE
## [1] 0.2172468

Model Autoregressive / Dynamic Regression

Model Autoregressive Distributed Lag (ARDL) adalah model regresi yang memasukkan nilai variabel yang menjelaskan baik nilai masa kini atau nilai masa lalu (lag) dari variabel independen sebagai tambahan pada model yang memasukkan nilai lag dari variabel dependen sebagai salah satu variabel penjelas (Islamiyahm D 2013).

model.ardl = ardlDlm(x = train.data$Xt, y = train.data$Yt, p = 1 , q = 1) #p:lag x, q:lag y
#model untuk p=1, q=1: yt=b0+b1yt-1+b2xt+b3xt-1

summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 2, End = 1462
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.722  -3.254   0.135   3.497  32.680 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.13099    1.21889   8.312   <2e-16 ***
## X.t         -3.17306    0.09522 -33.322   <2e-16 ***
## X.1          3.04419    0.09746  31.236   <2e-16 ***
## Y.1          0.88751    0.01150  77.183   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.025 on 1457 degrees of freedom
## Multiple R-squared:  0.8711, Adjusted R-squared:  0.8708 
## F-statistic:  3282 on 3 and 1457 DF,  p-value: < 2.2e-16

Pada pemodelan ardl didapatkan hasil yang signifikan untuk peubah X.t, X.1, dan Y.1 terhadap nilai Y (Kelembapan) pada taraf nyata 5%.

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

## $forecasts
##   [1] 78.83051 69.83431 74.13234 68.67718 69.65975 66.63600 81.40869 77.38977
##   [9] 80.32463 87.34857 89.61831 85.74523 80.12348 79.64748 68.85341 74.78413
##  [17] 83.90449 78.89947 73.46899 71.40746 71.13201 61.60394 63.45138 63.37931
##  [25] 55.60055 72.06926 70.98810 76.19933 73.29487 70.69377 71.76661 74.52274
##  [33] 68.46383 73.30596 63.45435 64.25438 70.28579 75.93056 73.74620 74.58922
##  [41] 76.59564 73.17919 71.08125 70.85752 69.19389 67.15653 58.96915 56.52516
##  [49] 58.10779 55.04938 52.89655 58.95648 68.23530 69.81217 68.10302 68.50486
##  [57] 69.27424 65.62633 54.83999 52.21500 55.57753 66.70820 71.05608 68.74593
##  [65] 67.80501 59.40169 63.36016 65.99768 68.71715 73.00408 76.20967 75.30811
##  [73] 72.95656 66.61127 53.59876 64.44384 52.59871 53.71194 58.85264 54.06768
##  [81] 47.75525 46.62257 44.99180 54.00324 46.39418 46.28677 46.17515 43.73431
##  [89] 50.39002 46.82875 45.36535 51.50594 49.62915 48.87164 54.92083 48.78760
##  [97] 62.80405 66.75346 61.31323 58.72814 54.24301 54.50255 53.41079 51.06603
## [105] 49.23022 50.50463 46.02331 42.19820 44.78710 42.38959 44.09925 49.13918
## [113] 49.58181 52.67498
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test.data$Xt, h = 114)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

#akurasi testing
mape.ardl <- MAPE(fore.ardl$forecasts, test.data$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.2971956
## 
## $MAPE_training.MAPE
## [1] 0.08092876

#penentuan lag optimum
ardlBoundOrders(data = data.frame(train.data) , formula = Yt ~ Xt )

## $p
##   Xt
## 1 15
## 
## $q
## [1] 10
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  9383.034 9365.121 9336.041 9321.435 9316.466 9312.353 9304.178 9297.530
## p = 2  9378.219 9347.557 9321.881 9306.246 9300.889 9296.794 9289.079 9282.262
## p = 3  9341.140 9341.140 9307.657 9294.615 9288.827 9284.657 9276.650 9270.052
## p = 4  9300.781 9302.728 9302.728 9294.548 9289.091 9284.885 9276.683 9269.936
## p = 5  9294.374 9292.917 9290.489 9290.489 9286.736 9282.401 9273.839 9266.811
## p = 6  9298.911 9290.558 9280.338 9281.964 9281.964 9283.736 9275.309 9268.160
## p = 7  9312.367 9298.427 9279.023 9277.481 9277.649 9277.649 9275.821 9269.080
## p = 8  9309.688 9295.104 9273.755 9270.901 9271.303 9272.285 9272.285 9269.815
## p = 9  9298.481 9285.432 9264.807 9261.946 9262.103 9263.344 9265.340 9265.340
## p = 10 9282.842 9270.586 9252.673 9250.461 9250.762 9252.340 9254.304 9256.228
## p = 11 9272.651 9259.127 9240.145 9238.649 9238.970 9240.743 9242.088 9243.949
## p = 12 9279.070 9262.333 9240.270 9237.617 9237.723 9239.649 9240.206 9240.709
## p = 13 9282.063 9263.465 9238.772 9235.168 9234.658 9236.591 9236.813 9236.633
## p = 14 9275.126 9257.862 9233.548 9229.679 9229.046 9230.958 9231.333 9231.143
## p = 15 9260.295 9242.059 9217.011 9211.299 9209.790 9211.773 9211.443 9211.271
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  9287.868 9281.382 9273.784 9268.417 9263.160 9258.179 9245.452
## p = 2  9272.544 9265.065 9257.877 9252.081 9246.781 9242.180 9229.796
## p = 3  9260.153 9252.612 9243.826 9238.055 9232.348 9227.954 9215.079
## p = 4  9260.197 9252.576 9243.666 9237.348 9231.658 9227.446 9214.281
## p = 5  9256.713 9249.210 9240.170 9233.710 9227.103 9222.840 9208.753
## p = 6  9257.884 9250.261 9241.135 9234.606 9227.928 9223.872 9209.911
## p = 7  9258.540 9250.737 9241.559 9235.007 9228.236 9224.227 9209.934
## p = 8  9259.693 9251.724 9242.567 9235.915 9229.120 9225.143 9210.759
## p = 9  9260.353 9252.805 9243.636 9236.898 9230.011 9226.068 9211.661
## p = 10 9256.228 9254.401 9245.367 9238.569 9231.639 9227.708 9213.177
## p = 11 9245.571 9245.571 9247.289 9240.556 9233.610 9229.680 9215.159
## p = 12 9240.345 9241.366 9241.366 9242.514 9235.476 9231.541 9216.931
## p = 13 9234.522 9235.107 9237.072 9237.072 9237.099 9233.109 9218.613
## p = 14 9228.780 9229.201 9231.182 9233.115 9233.115 9235.099 9220.573
## p = 15 9207.961 9207.519 9209.392 9211.180 9211.587 9211.587 9213.423
## 
## $min.Stat
## [1] 9207.519

Syntax di atas melakukan looping dari p=1 hingga p=15 dengan masing-masing p terdapat looping q=1 hingga q=15, didapatkan nilai AIC terkecil sebesar 323.0135 pada p=15 dan q=10.

Pemodelan DLM dan ARDL

#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)

#Ringkasan Model
summary(cons_lm1)

## 
## Time series regression with "ts" data:
## Start = 2, End = 1462
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.586  -9.171  -0.383   9.077  35.719 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  93.1386     1.2936   72.00  < 2e-16 ***
## Xt           -2.6185     0.2141  -12.23  < 2e-16 ***
## L(Xt)         1.3489     0.2141    6.30 3.94e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.59 on 1458 degrees of freedom
## Multiple R-squared:  0.3441, Adjusted R-squared:  0.3432 
## F-statistic: 382.4 on 2 and 1458 DF,  p-value: < 2.2e-16

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 1462
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.717  -4.675  -0.409   4.050  41.787 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.61327    1.52483  12.863   <2e-16 ***
## Xt          -0.30785    0.03301  -9.325   <2e-16 ***
## L(Yt)        0.80656    0.01447  55.737   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.782 on 1458 degrees of freedom
## Multiple R-squared:  0.7848, Adjusted R-squared:  0.7845 
## F-statistic:  2658 on 2 and 1458 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 1462
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.722  -3.254   0.135   3.497  32.680 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.13099    1.21889   8.312   <2e-16 ***
## Xt          -3.17306    0.09522 -33.322   <2e-16 ***
## L(Xt)        3.04419    0.09746  31.236   <2e-16 ***
## L(Yt)        0.88751    0.01150  77.183   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.025 on 1457 degrees of freedom
## Multiple R-squared:  0.8711, Adjusted R-squared:  0.8708 
## F-statistic:  3282 on 3 and 1457 DF,  p-value: < 2.2e-16

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 1462
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.730  -8.569  -0.331   9.204  33.463 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  92.6051     1.2938  71.574  < 2e-16 ***
## Xt           -2.7561     0.2155 -12.788  < 2e-16 ***
## L(Xt)         0.5628     0.2783   2.023   0.0433 *  
## L(Xt, 2)      0.9445     0.2159   4.376  1.3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.51 on 1456 degrees of freedom
## Multiple R-squared:  0.3511, Adjusted R-squared:  0.3498 
## F-statistic: 262.6 on 3 and 1456 DF,  p-value: < 2.2e-16

#SSE
deviance(cons_lm1)

## [1] 269121.4

deviance(cons_lm2)

## [1] 88301.73

deviance(cons_lm3)

## [1] 52886.73

deviance(cons_lm4)

## [1] 265606.2

Uji Diagnostik Model

Uji Non Autokorelasi

\[H_0: \text{Tidak terdeteksi autokorelasi}\] \[H_1: \text{Terdeteksi autokorelasi}\]

#Diagnostik
#durbin watson
dwtest(cons_lm1)

## 
##  Durbin-Watson test
## 
## data:  cons_lm1
## DW = 0.22587, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm2)

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

dwtest(cons_lm3)

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

dwtest(cons_lm4)

## 
##  Durbin-Watson test
## 
## data:  cons_lm4
## DW = 0.22418, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Hanya model ke-2 yang tak tolak H0, artinya tidak terdeteksi autokorelasi pada model tersebut.

Uji Heterogenitas

\[H_0: \text{Ragam galat bersifat Homoskedastisitas}\] \[H_1: \text{Ragam galat bersifat Heteroskedastisitas}\]

bptest(cons_lm1)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm1
## BP = 104.94, df = 2, p-value < 2.2e-16

bptest(cons_lm2)

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

bptest(cons_lm3)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm3
## BP = 62.914, df = 3, p-value = 1.401e-13

bptest(cons_lm4)

## 
##  studentized Breusch-Pagan test
## 
## data:  cons_lm4
## BP = 117.14, df = 3, p-value < 2.2e-16

Hanya model ke-2 yang tak tolak H0, artinya ragam galat bersifat homoskedastisitas pada model tersebut.

Uji Normalitas

\[H_0: \text{Sisaan menyebar normal}\] \[H_1: \text{Sisaan tidak menyebar normal}\]

shapiro.test(residuals(cons_lm1))

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

shapiro.test(residuals(cons_lm2))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.98199, p-value = 1.506e-12

shapiro.test(residuals(cons_lm3))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm3)
## W = 0.97509, p-value = 3.272e-15

shapiro.test(residuals(cons_lm4))

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

Tidak ada model yang sisaannya menyebar normal karena nilai p-value seluruh model kurang taraf nyata 5%.

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.3675654
## DLM 1          0.3313371
## DLM 2          0.3318410
## Autoregressive 0.2971956

Nilai MAPE terkecil terdapat pada model autoregressive sebesar 29.7%, artinya model tersebut merupakan model terbaik diantara model lainnya. Keempat model memiliki nilai MAPE pada kisaran 29.7% sampai 36.7%. Menurut Lewis (1982) dalam Kasemset et al. (2014), saat nilai MAPE berada di kisaran 21% hingga 50%, tingkat keakuratan forecast masih berada di tingkat wajar.

#PLOT
par(mfrow=c(1,1))
plot(test.data$Xt, test.data$Yt, type="b", col="black")
points(test.data$Xt, fore.koyck$forecasts,col="red")
lines(test.data$Xt, fore.koyck$forecasts,col="red")
points(test.data$Xt, fore.dlm$forecasts,col="blue")
lines(test.data$Xt, fore.dlm$forecasts,col="blue")
points(test.data$Xt, fore.dlm2$forecasts,col="orange")
lines(test.data$Xt, fore.dlm2$forecasts,col="orange")
points(test.data$Xt, fore.ardl$forecasts,col="green")
lines(test.data$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)

Secara eksploratif, model yang mendekati data aktual adalah model autoregressive, didukung dengan nilai MAPE terkecil dimiliki juga oleh model autoregressive.

Daftar Pustaka

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

Islamiyahm S. 2013. Penerapan Autoregressive Distributed Lag (ARDL) dalam Memodelkan Pengaruh Harga Minyak Dunia dan Jumlah Uang Beredar Terhadap Inflasi Di Indonesia. Jurnal Mahasiswa Statistik. 18(6): 45-48.

Pratami FR, Sudarno, Ispriyanti D. 2016. Peramalan dinamis produksi padi di Jawa Tengah menggunakan metode koyck dan almon. Jurnal Gaussian. 5(1): 91-97.

Virgantari F, Rahayu W. 2021. Pendugaan parameter model distributed lag pola polinimoal menggunakan metode almon. Jurnal Ilmu Matematika dan Terapan. 15(4): 761-772.

Tugas Individu 2 STA1341 MPDW