G1401221012_Tugas3_P2

Tugas 3 MPDW Individu

Ghonniyu Hiban Saputra - G1401221012

Data

Menggunakan data “Electric Motor Temperature” dengan peubah yang diambil adalah temperature stator winding (Y) dan temperature stator yoke (X), dimana pengambilan data sebanyak 200 periode (per 0.5 second).

Sumber : https://www.kaggle.com/datasets/wkirgsn/electric-motor-temperature

data <- read.csv("C:/Users/Ghonniyu/Downloads/DataMPDWIndividu3.csv")
str(data)

## 'data.frame':    200 obs. of  4 variables:
##  $ t     : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Yt    : num  51.6 51.9 51.9 51.9 52 ...
##  $ Y.t.1.: num  NA 51.6 51.9 51.9 51.9 ...
##  $ Xt    : num  29.6 29.7 29.7 29.7 29.8 ...

data

##       t       Yt   Y.t.1.       Xt
## 1     1 51.64935       NA 29.62973
## 2     2 51.88856 51.64935 29.65727
## 3     3 51.93739 51.88856 29.69390
## 4     4 51.87495 51.93739 29.73501
## 5     5 51.98421 51.87495 29.76484
## 6     6 52.00030 51.98421 29.78189
## 7     7 52.04716 52.00030 29.79709
## 8     8 51.88370 52.04716 29.81261
## 9     9 51.80456 51.88370 29.82353
## 10   10 51.93914 51.80456 29.83101
## 11   11 52.11126 51.93914 29.83247
## 12   12 52.18331 52.11126 29.82777
## 13   13 52.17953 52.18331 29.84302
## 14   14 52.30742 52.17953 29.88530
## 15   15 52.40242 52.30742 29.92845
## 16   16 52.45344 52.40242 29.96160
## 17   17 52.48754 52.45344 29.98539
## 18   18 52.51979 52.48754 30.00696
## 19   19 52.57042 52.51979 30.04093
## 20   20 52.59512 52.57042 30.07163
## 21   21 52.60328 52.59512 30.09341
## 22   22 52.63738 52.60328 30.12233
## 23   23 52.62694 52.63738 30.15196
## 24   24 52.64184 52.62694 30.17527
## 25   25 52.63043 52.64184 30.18471
## 26   26 52.63947 52.63043 30.18967
## 27   27 52.72278 52.63947 30.19959
## 28   28 52.75023 52.72278 30.21659
## 29   29 52.71454 52.75023 30.24574
## 30   30 52.71077 52.71454 30.28457
## 31   31 52.73590 52.71077 30.33292
## 32   32 52.76415 52.73590 30.37088
## 33   33 52.87951 52.76415 30.39352
## 34   34 52.95693 52.87951 30.40701
## 35   35 53.05413 52.95693 30.41308
## 36   36 53.11330 53.05413 30.43138
## 37   37 53.15305 53.11330 30.46821
## 38   38 53.13365 53.15305 30.51846
## 39   39 53.09045 53.13365 30.55644
## 40   40 53.09690 53.09045 30.55913
## 41   41 53.03942 53.09690 30.56052
## 42   42 53.05008 53.03942 30.57587
## 43   43 53.16429 53.05008 30.62524
## 44   44 53.22357 53.16429 30.67208
## 45   45 53.23690 53.22357 30.69213
## 46   46 53.26553 53.23690 30.78539
## 47   47 53.29998 53.26553 30.91174
## 48   48 53.31866 53.29998 30.99870
## 49   49 53.35300 53.31866 31.06680
## 50   50 53.48185 53.35300 31.09215
## 51   51 53.47282 53.48185 31.10514
## 52   52 53.44607 53.47282 31.14893
## 53   53 53.45237 53.44607 31.18732
## 54   54 53.48850 53.45237 31.21110
## 55   55 53.54290 53.48850 31.23894
## 56   56 53.65436 53.54290 31.26644
## 57   57 53.77165 53.65436 31.28615
## 58   58 53.83664 53.77165 31.29518
## 59   59 53.94035 53.83664 31.30069
## 60   60 54.06034 53.94035 31.31123
## 61   61 54.07041 54.06034 31.31824
## 62   62 54.09831 54.07041 31.32582
## 63   63 54.10067 54.09831 31.33406
## 64   64 54.10302 54.10067 31.34410
## 65   65 54.21091 54.10302 31.35638
## 66   66 54.26326 54.21091 31.36367
## 67   67 54.27193 54.26326 31.36326
## 68   68 54.31049 54.27193 31.35734
## 69   69 54.44696 54.31049 31.35673
## 70   70 54.54250 54.44696 31.36547
## 71   71 54.55907 54.54250 31.37760
## 72   72 54.51676 54.55907 31.38710
## 73   73 54.63171 54.51676 31.40066
## 74   74 54.73656 54.63171 31.41246
## 75   75 54.61255 54.73656 31.42257
## 76   76 54.55753 54.61255 31.43468
## 77   77 54.74929 54.55753 31.43784
## 78   78 54.87251 54.74929 31.43053
## 79   79 54.89620 54.87251 31.42800
## 80   80 54.95148 54.89620 31.43332
## 81   81 55.02780 54.95148 31.44042
## 82   82 55.19744 55.02780 31.46481
## 83   83 55.23372 55.19744 31.49809
## 84   84 55.15863 55.23372 31.52379
## 85   85 55.20731 55.15863 31.54869
## 86   86 55.37567 55.20731 31.56695
## 87   87 55.50696 55.37567 31.59035
## 88   88 55.60376 55.50696 31.62208
## 89   89 55.67872 55.60376 31.63913
## 90   90 55.72608 55.67872 31.65913
## 91   91 55.75063 55.72608 31.69776
## 92   92 55.77134 55.75063 31.73860
## 93   93 55.79393 55.77134 31.77611
## 94   94 55.77692 55.79393 31.78921
## 95   95 55.75585 55.77692 31.80034
## 96   96 55.79232 55.75585 31.81236
## 97   97 55.83733 55.79232 31.81974
## 98   98 55.83512 55.83733 31.82318
## 99   99 55.82227 55.83512 31.84167
## 100 100 55.80811 55.82227 31.89968
## 101 101 55.81616 55.80811 31.94859
## 102 102 55.83844 55.81616 31.96902
## 103 103 55.92384 55.83844 31.99486
## 104 104 55.99514 55.92384 32.02884
## 105 105 55.99227 55.99514 32.06041
## 106 106 56.00190 55.99227 32.08252
## 107 107 55.95761 56.00190 32.09864
## 108 108 55.92881 55.95761 32.11057
## 109 109 55.93888 55.92881 32.12457
## 110 110 56.01478 55.93888 32.15296
## 111 111 56.03502 56.01478 32.17449
## 112 112 56.08849 56.03502 32.18989
## 113 113 56.23805 56.08849 32.20439
## 114 114 56.27985 56.23805 32.21506
## 115 115 56.24274 56.27985 32.23447
## 116 116 56.28813 56.24274 32.23914
## 117 117 56.24395 56.28813 32.23262
## 118 118 56.22689 56.24395 32.23417
## 119 119 56.23465 56.22689 32.23948
## 120 120 56.34801 56.23465 32.24382
## 121 121 56.42891 56.34801 32.25118
## 122 122 56.45140 56.42891 32.26058
## 123 123 56.41881 56.45140 32.27291
## 124 124 56.37491 56.41881 32.28584
## 125 125 56.39643 56.37491 32.29711
## 126 126 56.48003 56.39643 32.30494
## 127 127 56.56015 56.48003 32.31563
## 128 128 56.59611 56.56015 32.35373
## 129 129 56.77263 56.59611 32.39595
## 130 130 56.83879 56.77263 32.41609
## 131 131 56.96176 56.83879 32.43014
## 132 132 57.06023 56.96176 32.44945
## 133 133 57.15655 57.06023 32.46421
## 134 134 57.16453 57.15655 32.48189
## 135 135 57.16249 57.16453 32.49259
## 136 136 57.26912 57.16249 32.50183
## 137 137 57.32336 57.26912 32.49661
## 138 138 57.29058 57.32336 32.50792
## 139 139 57.23238 57.29058 32.53847
## 140 140 57.31894 57.23238 32.58742
## 141 141 57.37175 57.31894 32.62932
## 142 142 57.43402 57.37175 32.63633
## 143 143 57.48459 57.43402 32.64416
## 144 144 57.46504 57.48459 32.66977
## 145 145 57.39945 57.46504 32.69297
## 146 146 57.40333 57.39945 32.70549
## 147 147 57.46579 57.40333 32.75422
## 148 148 57.46965 57.46579 32.78222
## 149 149 57.47084 57.46965 32.79845
## 150 150 57.50290 57.47084 32.79461
## 151 151 57.60053 57.50290 32.79073
## 152 152 57.59484 57.60053 32.83928
## 153 153 57.62471 57.59484 32.85749
## 154 154 57.69783 57.62471 32.85516
## 155 155 57.70847 57.69783 32.87476
## 156 156 57.66739 57.70847 32.90038
## 157 157 57.64862 57.66739 32.92482
## 158 158 57.64858 57.64862 32.93440
## 159 159 57.68982 57.64858 32.94460
## 160 160 57.85397 57.68982 32.94647
## 161 161 57.95845 57.85397 32.95995
## 162 162 58.04887 57.95845 32.97980
## 163 163 58.07821 58.04887 33.02584
## 164 164 58.02195 58.07821 33.06006
## 165 165 58.02968 58.02195 33.08887
## 166 166 58.08735 58.02968 33.10133
## 167 167 58.05397 58.08735 33.10937
## 168 168 58.03453 58.05397 33.11237
## 169 169 58.01578 58.03453 33.11636
## 170 170 58.05829 58.01578 33.12295
## 171 171 58.18425 58.05829 33.13206
## 172 172 58.23978 58.18425 33.15810
## 173 173 58.31623 58.23978 33.19416
## 174 174 58.32856 58.31623 33.22267
## 175 175 58.41167 58.32856 33.23521
## 176 176 58.45676 58.41167 33.25647
## 177 177 58.43987 58.45676 33.27434
## 178 178 58.54255 58.43987 33.28341
## 179 179 58.56501 58.54255 33.29035
## 180 180 58.59364 58.56501 33.30054
## 181 181 58.59679 58.59364 33.31789
## 182 182 58.59422 58.59679 33.31764
## 183 183 58.64447 58.59422 33.30633
## 184 184 58.55604 58.64447 33.30541
## 185 185 58.59396 58.55604 33.31520
## 186 186 58.57676 58.59396 33.33400
## 187 187 58.65994 58.57676 33.35766
## 188 188 58.68578 58.65994 33.42064
## 189 189 58.65026 58.68578 33.47199
## 190 190 58.73383 58.65026 33.51431
## 191 191 58.86605 58.73383 33.54534
## 192 192 58.99263 58.86605 33.57129
## 193 193 59.05439 58.99263 33.58095
## 194 194 59.06970 59.05439 33.60361
## 195 195 59.03437 59.06970 33.62875
## 196 196 58.99072 59.03437 33.63979
## 197 197 58.98935 58.99072 33.64371
## 198 198 59.02021 58.98935 33.64818
## 199 199 59.04424 59.02021 33.63995
## 200 200 59.07207 59.04424 33.65342

Packages

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

## 
## Attaching package: 'MLmetrics'

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

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

## Loading required package: carData

Pembagian Data

#SPLIT DATA
train<-data[1:160,]
test<-data[161:200,]

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

Pemodelan

Pemodelan model Koyck dengan R dapat menggunakan dLagM::koyckDlm() . Fungsi umum dari koyckDlm adalah sebagai berikut.

#MODEL KOYCK
model.koyck <- koyckDlm(x = train$Xt, y = train$Yt)
summary(model.koyck)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.199265 -0.037980 -0.007312  0.043446  0.200607 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.005562   0.176788   0.031    0.975    
## Y.1         0.976676   0.014901  65.545   <2e-16 ***
## X.t         0.041734   0.028664   1.456    0.147    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06238 on 156 degrees of freedom
## Multiple R-Squared: 0.9989,  Adjusted R-squared: 0.9989 
## Wald test: 6.913e+04 on 2 and 156 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                              alpha       beta       phi
## Geometric coefficients:  0.2384473 0.04173386 0.9766759

AIC(model.koyck)

## [1] -426.1257

BIC(model.koyck)

## [1] -413.8501

Dari hasil tersebut, didapat bahwa peubah \(y_{t-1}\) memiliki nilai \(P-Value<0.05\). Hal ini menunjukkan bahwa peubah \(y_{t-1}\) berpengaruh signifikan terhadap \(y\). Adapun model keseluruhannya adalah sebagai berikut

\[ \hat{Y_t}=0.0055+0.041X_t+0.97Y_{t-1} \]

Peramalan dan Akurasi

Berikut adalah hasil peramalan y untuk 40 periode kedepan menggunakan model koyck

fore.koyck <- forecast(model = model.koyck, x=test$Xt, h=40) #h=40 sesuai banyaknya data test
fore.koyck

## $forecasts
##  [1] 57.88569 57.91749 57.95048 57.98412 58.01819 58.05198 58.08531 58.11799
##  [9] 58.15008 58.18169 58.21295 58.24457 58.27695 58.30976 58.34234 58.37504
## [17] 58.40773 58.44003 58.47186 58.50339 58.53489 58.56566 58.59523 58.62408
## [25] 58.65266 58.68136 58.71038 58.74135 58.77374 58.80714 58.84105 58.87526
## [33] 58.90908 58.94305 58.97728 59.01117 59.04443 59.07710 59.10867 59.14007
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 40)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

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

## [1] 0.001178709

#akurasi data training
GoF(model.koyck) #mape data training

##               n        MAE           MPE         MAPE        sMAPE      MASE
## model.koyck 159 0.04856078 -1.263817e-06 0.0008874618 0.0008875074 0.8784775
##                     MSE     MRAE     GMRAE
## model.koyck 0.003817241 8.043174 0.9216455

Regression with Distributed Lag

Fungsi dlm() akan menerapkan model lag terdistribusi dengan satu atau lebih prediktor. Nilai x dan y tidak perlu sebagai objek time series (ts). \(q\) adalah integer yang mewakili panjang lag yang terbatas.

Pemodelan (Lag=2)

model.dlm <- dlm(x = train$Xt,y = train$Yt , q = 2) #q adalah lag, 
summary(model.dlm)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.92696 -0.17066  0.05731  0.24593  0.48539 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.3393     0.8938  -4.855 2.93e-06 ***
## x.t           2.2364     1.9008   1.177   0.2412    
## x.1          -3.7499     3.4501  -1.087   0.2788    
## x.2           3.3983     1.8974   1.791   0.0753 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3308 on 154 degrees of freedom
## Multiple R-squared:  0.9682, Adjusted R-squared:  0.9675 
## F-statistic:  1561 on 3 and 154 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 104.7959 120.1089

AIC(model.dlm)

## [1] 104.7959

BIC(model.dlm)

## [1] 120.1089

Dari hasil diatas, didapat bahwa \(P-value\) dari intercept \(<0.05\). Hal ini menunjukkan bahwa intercept berpengaruh signifikan terhadap \(y\). Adapun model keseluruhan yang terbentuk adalah sebagai berikut

\[ \hat{Y_t}=-4.3393+2.2364X_t-3.7499X_{t-1}+3.3983X_{t-2} \]

Peramalan dan Akurasi

Berikut merupakan hasil peramalan \(y\) untuk 40 periode kedepan

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

## $forecasts
##  [1] 57.78307 57.78325 57.85763 57.82894 57.92152 57.95762 58.02683 58.04569
##  [9] 58.07071 58.08065 58.08989 58.13636 58.15032 58.16733 58.21104 58.30843
## [17] 58.31128 58.33684 58.37904 58.40665 58.43079 58.39985 58.43442 58.47394
## [25] 58.46082 58.46307 58.47872 58.59475 58.55383 58.66993 58.75513 58.84066
## [33] 58.87036 58.97303 58.97707 58.98454 59.03732 59.07013 59.04828 59.12448
## 
## $call
## forecast.dlm(model = model.dlm, x = test$Xt, h = 40)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

mape.dlm <- MAPE(fore.dlm$forecasts, test$Yt)
#akurasi data training
GoF(model.dlm)

##             n       MAE         MPE        MAPE       sMAPE     MASE      MSE
## model.dlm 158 0.2531926 -3.3024e-05 0.004647291 0.004641742 4.576917 0.106682
##               MRAE    GMRAE
## model.dlm 15.98012 4.631467

Lag Optimum

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

##    q - k    MASE       AIC       BIC   GMRAE  MBRAE R.Adj.Sq    Ljung-Box
## 78    78 0.04652 -559.3937 -364.4494 0.05024 -0.015  0.99947 5.448975e-13

#qmin, qmax dicarinya dari coba-coba
#tracing ubah trace jadi TRUE dan liat

Berdasarkan output tersebut, lag optimum didapatkan ketika lag=78. Selanjutnya dilakukan pemodelan untuk lag=78

#model dlm dengan lag optimum
model.dlm2 <- dlm(x = train$Xt,y = train$Yt , q = 78)
summary(model.dlm2)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##          1          2          3          4          5          6          7 
## -6.987e-04  7.329e-04 -3.441e-03  6.360e-03 -3.938e-03  3.883e-03 -3.807e-03 
##          8          9         10         11         12         13         14 
##  2.559e-03 -1.173e-03 -3.281e-03  4.451e-03 -1.990e-03  3.268e-03 -4.286e-03 
##         15         16         17         18         19         20         21 
##  8.974e-04 -1.759e-04  7.523e-04 -6.252e-04  1.201e-03  3.213e-04 -3.042e-03 
##         22         23         24         25         26         27         28 
##  1.606e-03 -2.088e-03  2.417e-03 -1.596e-03  5.447e-03 -6.070e-03  4.083e-03 
##         29         30         31         32         33         34         35 
## -3.842e-03  4.124e-03 -3.678e-03  3.133e-03 -3.498e-03  4.071e-03 -4.609e-03 
##         36         37         38         39         40         41         42 
##  2.749e-03 -1.843e-03  1.820e-03 -3.949e-04 -3.691e-04  1.303e-03 -2.467e-03 
##         43         44         45         46         47         48         49 
##  3.211e-03 -3.854e-03  5.954e-03 -4.431e-03  1.683e-03 -1.261e-03  3.234e-05 
##         50         51         52         53         54         55         56 
## -2.049e-03  2.759e-03 -9.609e-04 -3.225e-03  6.729e-03 -1.158e-03 -9.695e-04 
##         57         58         59         60         61         62         63 
##  1.128e-03 -5.090e-03  3.376e-03 -3.424e-04  4.456e-04 -1.541e-04  1.299e-04 
##         64         65         66         67         68         69         70 
## -5.132e-05 -3.240e-03  3.434e-03 -3.132e-03  5.984e-03 -3.557e-03 -6.156e-04 
##         71         72         73         74         75         76         77 
##  1.911e-03 -5.305e-04 -2.987e-03  2.959e-03 -6.412e-04  2.744e-03  5.648e-05 
##         78         79         80         81         82 
##  2.242e-04 -2.704e-04 -2.437e-03  1.713e-03 -1.783e-03 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 13.14204    3.53182   3.721   0.0652 .
## x.t         -3.07437    0.76866  -4.000   0.0572 .
## x.1         -1.23744    1.05915  -1.168   0.3631  
## x.2         -2.74079    1.28060  -2.140   0.1657  
## x.3         -1.67364    1.42927  -1.171   0.3622  
## x.4         -1.38479    1.36363  -1.016   0.4167  
## x.5          0.08722    1.47092   0.059   0.9581  
## x.6          1.18627    1.60314   0.740   0.5364  
## x.7          1.42644    1.63290   0.874   0.4745  
## x.8         -0.04152    1.63504  -0.025   0.9820  
## x.9          3.19314    1.75747   1.817   0.2109  
## x.10        -5.90813    1.88989  -3.126   0.0889 .
## x.11         0.96255    1.87652   0.513   0.6590  
## x.12         1.20703    1.78984   0.674   0.5696  
## x.13        -0.95404    1.80491  -0.529   0.6499  
## x.14         2.99483    2.07466   1.444   0.2857  
## x.15        -2.86336    2.11461  -1.354   0.3084  
## x.16         1.34398    2.10598   0.638   0.5887  
## x.17         0.49706    2.11511   0.235   0.8361  
## x.18        -1.31718    2.11659  -0.622   0.5972  
## x.19         0.24454    2.10400   0.116   0.9181  
## x.20         0.21653    2.16146   0.100   0.9293  
## x.21         2.72185    2.15427   1.263   0.3338  
## x.22        -1.07867    1.99929  -0.540   0.6436  
## x.23         2.53744    2.04540   1.241   0.3406  
## x.24        -5.78306    2.36916  -2.441   0.1347  
## x.25         1.23847    2.39101   0.518   0.6561  
## x.26         0.93721    2.29978   0.408   0.7231  
## x.27        -2.62714    2.28258  -1.151   0.3688  
## x.28         4.65880    2.19265   2.125   0.1675  
## x.29         2.12983    2.23170   0.954   0.4406  
## x.30         1.66073    2.24331   0.740   0.5362  
## x.31        -1.73827    2.39062  -0.727   0.5427  
## x.32         1.75447    2.46625   0.711   0.5506  
## x.33        -2.73417    2.50064  -1.093   0.3883  
## x.34         0.39557    2.47325   0.160   0.8876  
## x.35        -2.13802    2.20021  -0.972   0.4337  
## x.36         2.13438    2.13169   1.001   0.4222  
## x.37        -0.84298    2.19214  -0.385   0.7376  
## x.38         0.01486    2.30293   0.006   0.9954  
## x.39         1.06058    2.46996   0.429   0.7095  
## x.40         0.21985    2.57854   0.085   0.9398  
## x.41         0.88026    2.70698   0.325   0.7759  
## x.42        -0.79897    2.84107  -0.281   0.8050  
## x.43         2.93783    2.88642   1.018   0.4159  
## x.44        -0.57424    2.96275  -0.194   0.8642  
## x.45         0.35202    2.94055   0.120   0.9157  
## x.46        -0.70268    2.98768  -0.235   0.8359  
## x.47         1.26767    3.10755   0.408   0.7228  
## x.48        -2.00659    3.25832  -0.616   0.6008  
## x.49         1.87908    3.31130   0.567   0.6276  
## x.50        -1.60478    3.38742  -0.474   0.6824  
## x.51         2.91858    3.37667   0.864   0.4785  
## x.52        -3.40189    3.35264  -1.015   0.4170  
## x.53         0.40503    3.45146   0.117   0.9173  
## x.54         0.38261    3.41526   0.112   0.9210  
## x.55         0.50449    3.43169   0.147   0.8966  
## x.56         0.72482    3.51199   0.206   0.8556  
## x.57         1.25141    3.56737   0.351   0.7592  
## x.58         2.86393    3.71633   0.771   0.5215  
## x.59        -3.31293    3.89966  -0.850   0.4851  
## x.60         4.66440    4.01986   1.160   0.3657  
## x.61        -4.95514    4.01866  -1.233   0.3428  
## x.62         3.68702    3.98692   0.925   0.4527  
## x.63        -4.12049    4.00555  -1.029   0.4118  
## x.64         5.40929    4.08691   1.324   0.3167  
## x.65        -3.99632    4.13038  -0.968   0.4353  
## x.66         1.08434    3.99971   0.271   0.8117  
## x.67        -0.13958    3.83216  -0.036   0.9743  
## x.68        -0.53703    3.79664  -0.141   0.9005  
## x.69         1.19078    3.68527   0.323   0.7773  
## x.70        -2.04344    3.48348  -0.587   0.6169  
## x.71         1.19052    3.48209   0.342   0.7650  
## x.72         0.96670    3.43755   0.281   0.8050  
## x.73        -2.66694    3.20193  -0.833   0.4925  
## x.74         1.18663    2.78714   0.426   0.7117  
## x.75        -4.27205    2.34799  -1.819   0.2105  
## x.76         3.57124    1.94554   1.836   0.2078  
## x.77        -2.29771    1.46676  -1.567   0.2577  
## x.78         2.97245    0.76929   3.864   0.0609 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01905 on 2 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:  0.9995 
## F-statistic:  1926 on 79 and 2 DF,  p-value: 0.000519
## 
## AIC and BIC values for the model:
##         AIC       BIC
## 1 -559.3937 -364.4494

AIC(model.dlm2)

## [1] -559.3937

BIC(model.dlm2)

## [1] -364.4494

Dari hasil tersebut terdapat beberapa peubah yang berpengaruh signifikan terhadap taraf nyata 1% yaitu \(x_t\) , \(x_{t-10}\) , \(x_{t-78}\). Adapun keseluruhan model yang terbentuk adalah

\[ \hat{Y_t}=13.14204-3.07437X_t+...+2.97245X_{t-78} \]

Adapun hasil peramalan 40 periode kedepan menggunakan model tersebut adalah sebagai berikut

#peramalan dan akurasi
fore.dlm2 <- forecast(model = model.dlm2, x=test$Xt, h=40)
mape.dlm2<- MAPE(fore.dlm2$forecasts, test$Yt)
mape.dlm2 #data test

## [1] 0.007521861

#akurasi data training
GoF(model.dlm2)

##             n         MAE           MPE         MAPE       sMAPE       MASE
## model.dlm2 82 0.002430579 -2.913616e-09 4.313069e-05 4.31308e-05 0.04652476
##                     MSE      MRAE      GMRAE
## model.dlm2 8.847629e-06 0.9329747 0.05024046

Model tersebut merupakan model yang sangat baik dengan nilai MAPE yang kurang dari 10%.

Model Autoregressive

Peubah dependen dipengaruhi oleh peubah independen pada waktu sekarang, serta dipengaruhi juga oleh peubah dependen itu sendiri pada satu waktu yang lalu maka model tersebut disebut autoregressive (Gujarati 2004).

Pemodelan

Pemodelan Autoregressive dilakukan menggunakan fungsi dLagM::ardlDlm() . Fungsi tersebut akan menerapkan autoregressive berordo \((p,q)\) dengan satu prediktor.

Dengan \(p\) adalah integer yang mewakili panjang lag yang terbatas dan \(q\) adalah integer yang merepresentasikan ordo dari proses autoregressive.

model.ardl <- ardlDlm(x = train$Xt, y = train$Yt, p = 1 , q = 1)
summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 2, End = 160
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.202509 -0.038196 -0.002897  0.040807  0.200508 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.04509    0.17887   0.252    0.801    
## X.t         -0.23296    0.27208  -0.856    0.393    
## X.1          0.27325    0.27215   1.004    0.317    
## Y.1          0.97689    0.01488  65.652   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06237 on 155 degrees of freedom
## Multiple R-squared:  0.9989, Adjusted R-squared:  0.9989 
## F-statistic: 4.609e+04 on 3 and 155 DF,  p-value: < 2.2e-16

AIC(model.ardl)

## [1] -425.1679

BIC(model.ardl)

## [1] -409.8234

model.ardl <- ardlDlm(formula = Yt ~ Xt, 
                         data = train,p = 1 , q = 1)
summary(model.ardl)

## 
## Time series regression with "ts" data:
## Start = 2, End = 160
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.202509 -0.038196 -0.002897  0.040807  0.200508 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.04509    0.17887   0.252    0.801    
## Xt.t        -0.23296    0.27208  -0.856    0.393    
## Xt.1         0.27325    0.27215   1.004    0.317    
## Yt.1         0.97689    0.01488  65.652   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06237 on 155 degrees of freedom
## Multiple R-squared:  0.9989, Adjusted R-squared:  0.9989 
## F-statistic: 4.609e+04 on 3 and 155 DF,  p-value: < 2.2e-16

AIC(model.ardl)

## [1] -425.1679

BIC(model.ardl)

## [1] -409.8234

Hasil di atas menunjukkan \(x_t\) dan \(x_{t-1}\) tidak berpengaruh signifikan terhadap \(y_t\) dimana hasil uji t dengan nilai-p \(\ge0.05\). Sementara peubah \(y_{t-1}\) berpengaruh signifikan terhadap \(y_t\).

\[ \hat{Y}=0.04509-0.23296X_t+0.27325X_{t-1}+0.97689Y_{t-1} \]

Peramalan dan Akurasi

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

## $forecasts
##  [1] 57.88617 57.91669 57.94120 57.96975 58.00028 58.03507 58.07059 58.10679
##  [9] 58.14205 58.17604 58.20892 58.23747 58.26408 58.29328 58.32667 58.35777
## [17] 58.38979 58.42385 58.45797 58.49084 58.52168 58.55661 58.59331 58.62627
## [25] 58.65594 58.68323 58.70950 58.72697 58.74927 58.77523 58.80493 58.83638
## [33] 58.87194 58.90403 58.93573 58.97098 59.00753 59.04326 59.08130 59.11307
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 40)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

Data di atas merupakan hasil peramalan untuk 40 periode ke depan menggunakan Model Autoregressive dengan \(p=1\) dan \(q=1\).

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

## [1] 0.001248034

#akurasi data training
GoF(model.ardl)

##              n        MAE           MPE         MAPE        sMAPE      MASE
## model.ardl 159 0.04800685 -1.255336e-06 0.0008772722 0.0008773073 0.8684568
##                    MSE     MRAE     GMRAE
## model.ardl 0.003792302 8.304738 0.8728073

Berdasarkan akurasi di atas, terlihat bahwa nilai MAPE keduanya tidak jauh berbeda. Artinya, model regresi dengan distribusi lag ini overfitted atau underfitted

Lag Optimum

#penentuan lag optimum
model.ardl.opt <- ardlBoundOrders(data = data.frame(data), ic = "AIC", 
                                  formula = Yt ~ Xt )
min_p=c()
for(i in 1:10){
  min_p[i]=min(model.ardl.opt$Stat.table[[i]])
}
q_opt=which(min_p==min(min_p, na.rm = TRUE))
p_opt=which(model.ardl.opt$Stat.table[[q_opt]] == 
              min(model.ardl.opt$Stat.table[[q_opt]], na.rm = TRUE))
data.frame("q_optimum" = q_opt, "p_optimum" = p_opt, 
           "AIC"=model.ardl.opt$min.Stat)

##   q_optimum p_optimum       AIC
## 1         7         1 -579.2507

Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=1\) dan \(q=7\), yaitu sebesar -579.2507. Artinya, model autoregressive optimum didapat ketika \(p=1\) dan \(q=7\).

Selanjutnya dapat dilakukan pemodelan dengan nilai \(p\) dan \(q\) optimum seperti inisialisasi di langkah sebelumnya.

Pemodelan DLM & ARDL dengan Library `dynlm`

Pemodelan regresi dengan peubah lag tidak hanya dapat dilakukan dengan fungsi pada packages dLagM , tetapi terdapat packages dynlm yang dapat digunakan.

Untuk menentukan formula model yang akan digunakan, tersedia fungsi tambahan yang memungkinkan spesifikasi dinamika (melalui d() dan L()) atau pola linier/siklus dengan mudah (melalui trend(), season(), dan harmon()). Semua fungsi formula baru mengharuskan argumennya berupa objek deret waktu (yaitu, "ts" atau "zoo").

#sama dengan model dlm q=1
cons_lm1 <- dynlm(Yt ~ Xt+L(Xt),data = train.ts)
#sama dengan model ardl p=1 q=0
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 = 160
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93350 -0.15772  0.06111  0.25153  0.45356 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.29182    0.88931  -4.826 3.29e-06 ***
## Xt           0.01894    1.45551   0.013    0.990    
## L(Xt)        1.86367    1.45024   1.285    0.201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3337 on 156 degrees of freedom
## Multiple R-squared:  0.9677, Adjusted R-squared:  0.9673 
## F-statistic:  2340 on 2 and 156 DF,  p-value: < 2.2e-16

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 160
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Yt), data = train.ts)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.199900 -0.038392 -0.005869  0.043529  0.200117 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.01584    0.17649   0.090    0.929    
## Xt           0.03872    0.02851   1.358    0.176    
## L(Yt)        0.97822    0.01482  66.002   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06237 on 156 degrees of freedom
## Multiple R-squared:  0.9989, Adjusted R-squared:  0.9989 
## F-statistic: 6.914e+04 on 2 and 156 DF,  p-value: < 2.2e-16

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 160
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.202509 -0.038196 -0.002897  0.040807  0.200508 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.04509    0.17887   0.252    0.801    
## Xt          -0.23296    0.27208  -0.856    0.393    
## L(Xt)        0.27325    0.27215   1.004    0.317    
## L(Yt)        0.97689    0.01488  65.652   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06237 on 155 degrees of freedom
## Multiple R-squared:  0.9989, Adjusted R-squared:  0.9989 
## F-statistic: 4.609e+04 on 3 and 155 DF,  p-value: < 2.2e-16

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 160
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2), data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.92696 -0.17066  0.05731  0.24593  0.48539 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.3393     0.8938  -4.855 2.93e-06 ***
## Xt            2.2364     1.9008   1.177   0.2412    
## L(Xt)        -3.7499     3.4501  -1.087   0.2788    
## L(Xt, 2)      3.3983     1.8974   1.791   0.0753 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3308 on 154 degrees of freedom
## Multiple R-squared:  0.9682, Adjusted R-squared:  0.9675 
## F-statistic:  1561 on 3 and 154 DF,  p-value: < 2.2e-16

SSE

deviance(cons_lm1)

## [1] 17.37057

deviance(cons_lm2)

## [1] 0.6068978

deviance(cons_lm3)

## [1] 0.6029761

deviance(cons_lm4)

## [1] 16.85575

Uji Diagnostik

#uji model
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    155 -1 4310.2505 <2e-16 ***
## M2 vs. ME    155 -1    1.0081 0.3169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Autokorelasi

#durbin watson
dwtest(cons_lm1)

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

dwtest(cons_lm2)

## 
##  Durbin-Watson test
## 
## data:  cons_lm2
## DW = 1.3881, p-value = 2.622e-05
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm3)

## 
##  Durbin-Watson test
## 
## data:  cons_lm3
## DW = 1.3935, p-value = 2.445e-05
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm4)

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

Heterogenitas

bptest(cons_lm1)

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

bptest(cons_lm2)

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

bptest(cons_lm3)

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

bptest(cons_lm4)

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

Kenormalan

shapiro.test(residuals(cons_lm1))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm1)
## W = 0.91007, p-value = 2.48e-08

shapiro.test(residuals(cons_lm2))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm2)
## W = 0.98813, p-value = 0.1984

shapiro.test(residuals(cons_lm3))

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

shapiro.test(residuals(cons_lm4))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cons_lm4)
## W = 0.92172, p-value = 1.501e-07

Kesimpulan

Perbandingan Model

akurasi <- matrix(c(mape.koyck, mape.dlm, mape.dlm2, mape.ardl))
row.names(akurasi)<- c("Koyck","DLM 1","DLM 2","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi

##                       MAPE
## Koyck          0.001178709
## DLM 1          0.002087447
## DLM 2          0.007521861
## Autoregressive 0.001248034

Berdasarkan nilai MAPE, model paling optimum didapat pada Model Koyck karena memiliki nilai MAPE yang terkecil.

par(mfrow=c(1,1))
plot(test$Xt, test$Yt, type="b", col="black", ylim=c(56,62),main="Plot Perbandingan")
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)

Berdasarkan plot tersebut, terlihat bahwa plot yang paling mendekati data aktualnya adalah Model koyck, sehingga dapat disimpulkan model terbaik dalam kasus ini adalah model regresi koyck.

par(mfrow=c(1,1))
plot(test$Xt, test$Yt, type="b", col="black", ylim=c(57,60),main="Aktual vs Koycke")
points(test$Xt, fore.koyck$forecasts,col="red")
lines(test$Xt, fore.koyck$forecasts,col="red")
legend("topleft",c("aktual", "koyck"), lty=1, col=c("black","red"), cex=0.8)

G1401221012_Tugas3_P2