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

Impor Data

Peubah respon : Y Peubah penjelas : X

data <- import("https://raw.githubusercontent.com/khenihikmah/praktikummpdw/main/Pertemuan3/salesmonthly_P3.csv")
Yt = data$Y
Xt = data$X
data = data.frame(Yt,Xt)
str(data)

## 'data.frame':    70 obs. of  2 variables:
##  $ Yt: num  354 347 232 209 270 323 348 420 399 472 ...
##  $ Xt: int  50 31 20 18 23 23 21 29 14 30 ...

data

##       Yt Xt
## 1  354.0 50
## 2  347.0 31
## 3  232.0 20
## 4  209.0 18
## 5  270.0 23
## 6  323.0 23
## 7  348.0 21
## 8  420.0 29
## 9  399.0 14
## 10 472.0 30
## 11 489.0 19
## 12 492.0 25
## 13 463.0 24
## 14 243.0  9
## 15 208.0 13
## 16 192.0  5
## 17 194.0 10
## 18 217.0 12
## 19 203.0  6
## 20 265.5 15
## 21 243.5 11
## 22 222.0  8
## 23 228.0 18
## 24 286.0 28
## 25 248.0 24
## 26 239.0 20
## 27 250.0 13
## 28 318.0 18
## 29 275.0 18
## 30 311.0 20
## 31 240.0  8
## 32 275.5 12
## 33 307.0 18
## 34 312.0 11
## 35 246.0 27
## 36 257.0 18
## 37   1.0  0
## 38 144.0  7
## 39 165.0  9
## 40 132.0  9
## 41 148.0 23
## 42 163.0  8
## 43 219.0 15
## 44 239.0 12
## 45 223.0 23
## 46 226.0 15
## 47 192.0 15
## 48 226.0  6
## 49 229.0 11
## 50 268.0 12
## 51 381.0 42
## 52 289.0 21
## 53 259.0 13
## 54 248.0 18
## 55 283.0 19
## 56 253.0 20
## 57 263.0 12
## 58 287.0 25
## 59 252.2 22
## 60 254.0 27
## 61 295.2 23
## 62 249.4 12
## 63 301.4 19
## 64 299.4 22
## 65 265.8 26
## 66 193.0 25
## 67 250.6 20
## 68 237.0 26
## 69 227.8 16
## 70  86.0  7

Pembagian Data

#Split data dengan proporsi: 80% training dan 20% testing
train<-data[1:50,] # 80% data pertama
test<-data[51:64,]# 20% data sisanya

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

Model Koyck

Model Koyck didasarkan pada asumsi bahwa semakin jauh jarak lag peubah independen dari periode sekarang maka semakin kecil pengaruh peubah lag terhadap peubah dependen. Koyck mengusulkan suatu metode untuk menduga model dinamis distributed lag dengan mengasumsikan bahwa semua koefisien \(\beta\) mempunyai tanda sama.

Model kyock merupakan jenis paling umum dari model infinite distributed lag dan juga dikenal sebagai geometric lag

\[ y_t=a(1-\lambda)+\beta_0X_t+\beta_1Z_t+\lambda Y_{t-1}+V_t \]

dengan \[V_t=u_t-\lambda u_{t-1}\]

Pemodelan

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

koyckDlm(x , y , intercept)

Fungsi koyckDlm() akan menerapkan model lag terdistribusi dengan transformasi Koyck satu prediktor. Nilai x dan y tidak perlu sebagai objek time series (ts). intercept dapat dibuat TRUE untuk memasukkan intersep ke dalam model.

#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 
## -313.880  -27.133    1.345   52.619  143.246 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  90.1133    68.1314   1.323  0.19249   
## Y.1           0.8746     0.2906   3.009  0.00424 **
## X.t          -3.6772     8.0847  -0.455  0.65137   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 78.78 on 46 degrees of freedom
## Multiple R-Squared: 0.3399,  Adjusted R-squared: 0.3112 
## Wald test: 20.29 on 2 and 46 DF,  p-value: 4.817e-07 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha      beta       phi
## Geometric coefficients:  718.4886 -3.677212 0.8745794

AIC(model.koyck)

## [1] 571.8911

BIC(model.koyck)

## [1] 579.4584

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}=90.1133-3.6772 X_t+ 0.8746_{t-1} \]

Peramalan dan Akurasi Model Koyck

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

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

## $forecasts
##  [1] 170.0576 161.6207 183.6596 184.5484 181.6485 175.4350 199.4186 172.5903
##  [9] 160.1585 130.8999 120.0197 150.9535 152.2671 142.3842
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 14)
## 
## attr(,"class")
## [1] "forecast.koyckDlm" "dLagM"

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

#akurasi data training
mape.koyck.train <- GoF(model.koyck)["MAPE"]

c("MAPE Testing"=mape.koyck,"MAPE Training"=mape.koyck.train)

## $`MAPE Testing`
## [1] 0.4073759
## 
## $`MAPE Training.MAPE`
## [1] 6.582258

Regression with Distributed Lag

Pemodelan model Regression with Distributed Lag dengan R dapat menggunakan dLagM::dlm() . Fungsi umum dari dlm adalah sebagai berikut.

dlm(formula , data , x , y , q , remove )

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)```{r}

model.dlm <- 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 
## -148.741  -56.489    5.546   39.075  160.020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  96.6542    32.0225   3.018  0.00422 ** 
## x.t           7.6265     1.6093   4.739 2.27e-05 ***
## x.1           3.4174     1.6187   2.111  0.04048 *  
## x.2          -0.6756     1.3417  -0.504  0.61709    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 72.7 on 44 degrees of freedom
## Multiple R-squared:  0.4529, Adjusted R-squared:  0.4156 
## F-statistic: 12.14 on 3 and 44 DF,  p-value: 6.427e-06
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 553.5303 562.8863

AIC(model.dlm)

## [1] 553.5303

BIC(model.dlm)

## [1] 562.8863

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

\[ \hat{Y_t}= 96.6542+7.6265X_t+3.4174X_{t-1}+ -0.6756X_{t-2} \]

Peramalan dan Akurasi

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

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

## $forecasts
##  [1] 450.5461 392.2345 239.1893 264.1705 294.2887 301.9546 243.6840 314.8145
##  [9] 341.7656 360.8634 349.4709 248.5314 267.0284 321.2613
## 
## $call
## forecast.dlm(model = model.dlm, x = test$Xt, h = 14)
## 
## attr(,"class")
## [1] "forecast.dlm" "dLagM"

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

## [1] 0.1596735

mape_train <- GoF(model.dlm)["MAPE"]

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

## $MAPE_testing
## [1] 0.1596735
## 
## $MAPE_training.MAPE
## [1] 3.09374

#akurasi data training
GoF(model.dlm)

##            n      MAE       MPE    MAPE     sMAPE     MASE      MSE     MRAE
## model.dlm 48 54.75584 -2.931796 3.09374 0.2293388 1.305695 4844.842 4.134931
##              GMRAE
## model.dlm 1.328697

Lag Optimum

#penentuan lag optimum 
finiteDLMauto(formula = Yt ~ Xt,
              data = data.frame(train),
              model.type = "dlm", error.type = "AIC")

##    q - k    MASE      AIC      BIC   GMRAE    MBRAE R.Adj.Sq   Ljung-Box
## 10    10 0.94499 448.6853 470.6408 1.23383 -6.22705  0.59259 0.005845444

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

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

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -127.098  -30.418    4.495   38.131   94.663 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -91.0862    49.2543  -1.849   0.0750 .  
## x.t           8.1712     1.6434   4.972 2.99e-05 ***
## x.1           3.2053     1.5495   2.069   0.0479 *  
## x.2          -0.9618     1.5517  -0.620   0.5404    
## x.3           3.0833     1.4823   2.080   0.0468 *  
## x.4           2.6489     1.5013   1.764   0.0886 .  
## x.5           0.3288     1.5268   0.215   0.8310    
## x.6          -0.8076     1.5229  -0.530   0.6001    
## x.7           1.4702     1.5488   0.949   0.3506    
## x.8           1.9763     1.5773   1.253   0.2206    
## x.9           1.2692     1.5546   0.816   0.4212    
## x.10          1.5810     1.2610   1.254   0.2203    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 57 on 28 degrees of freedom
## Multiple R-squared:  0.7075, Adjusted R-squared:  0.5926 
## F-statistic: 6.157 on 11 and 28 DF,  p-value: 4.881e-05
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 448.6853 470.6408

AIC(model.dlm2)

## [1] 448.6853

BIC(model.dlm2)

## [1] 470.6408

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

\[ \hat{Y_t}= -91.0862 + 8.1712X_t+...+1.5810X_{t-10} \]

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

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

mape.dlm.train = GoF(model.dlm)["MAPE"]

c("MAPE Testing"=mape.dlm,"MAPE Training"=mape.dlm.train)

## $`MAPE Testing`
## [1] 0.1596735
## 
## $`MAPE Training.MAPE`
## [1] 3.09374

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

Model Autoregressive Distributed Lag (ARDL)

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. Fungsi umum dari ardlDlm() adalah sebagai berikut.

ardlDlm(formula = NULL , data = NULL , x = NULL , y = NULL , p = 1 , q = 1 , 
         remove = NULL )

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 = 50
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -169.13  -32.02   11.48   30.98   93.68 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.1359    23.7405   1.185    0.242    
## X.t           5.3093     1.1774   4.509 4.62e-05 ***
## X.1          -1.8098     1.1123  -1.627    0.111    
## Y.1           0.6792     0.1025   6.629 3.62e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 52.58 on 45 degrees of freedom
## Multiple R-squared:  0.7123, Adjusted R-squared:  0.6931 
## F-statistic: 37.14 on 3 and 45 DF,  p-value: 3.098e-12

AIC(model.ardl)

## [1] 533.1938

BIC(model.ardl)

## [1] 542.6529

Hasil di atas menunjukkan bahwa peubah \(x_{t-1}\), hasil uji t menunjukkan nilai-p pada peubah \(\ge0.05\) Hal ini menunjukkan bahwa peubah \(x_{t-1}\) tidak berpengaruh signifikan terhadap \(y_t\), sementara \(x_t\) dan \(y_{t-1}\) berpengaruh signifikan terhadap \(y_t\). Model keseluruhannya adalah sebagai berikut:

\[ \hat{Y}=28.1359+ 5.3093X_t -1.80982X_{t-1}+ 0.6792 Y_{t-1} \]

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:14){
  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        13         1 568.2968
## 2        14         1 568.2968

Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=13\) dan \(p=14\) serta \(q=1\), yaitu sebesar 568.2968. Artinya, model autoregressive optimum didapat ketika \(p=13\) dan \(p=14\) serta \(q=3\).

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

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

## 
## Time series regression with "ts" data:
## Start = 14, End = 50
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -107.954  -16.416    1.382   24.068   53.566 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 110.6781    57.5116   1.924  0.06795 . 
## X.t           4.5962     1.4952   3.074  0.00576 **
## X.1          -0.5402     1.8372  -0.294  0.77162   
## X.2          -2.5799     1.3940  -1.851  0.07834 . 
## X.3           2.6835     1.3416   2.000  0.05858 . 
## X.4          -0.2327     1.4173  -0.164  0.87116   
## X.5          -1.2504     1.2901  -0.969  0.34343   
## X.6          -0.1291     1.2776  -0.101  0.92047   
## X.7           0.8923     1.3147   0.679  0.50471   
## X.8           1.0029     1.3702   0.732  0.47229   
## X.9          -0.1689     1.3628  -0.124  0.90256   
## X.10          1.3475     1.3364   1.008  0.32477   
## X.11         -1.9013     1.3488  -1.410  0.17329   
## X.12         -0.9977     1.3290  -0.751  0.46117   
## X.13         -1.2990     1.0534  -1.233  0.23113   
## Y.1           0.4526     0.1796   2.521  0.01988 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43.73 on 21 degrees of freedom
## Multiple R-squared:  0.6893, Adjusted R-squared:  0.4674 
## F-statistic: 3.106 on 15 and 21 DF,  p-value: 0.008728

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

## 
## Time series regression with "ts" data:
## Start = 15, End = 50
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -87.64 -17.07  -0.15  26.23  52.23 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 179.0915    70.7902   2.530   0.0210 *
## X.t           4.1226     1.5013   2.746   0.0133 *
## X.1          -0.6491     1.8196  -0.357   0.7254  
## X.2          -4.9623     1.8180  -2.729   0.0138 *
## X.3           1.2462     1.5010   0.830   0.4173  
## X.4           0.4003     1.4750   0.271   0.7892  
## X.5          -2.5383     1.4198  -1.788   0.0907 .
## X.6          -0.6746     1.3393  -0.504   0.6206  
## X.7           0.8747     1.3226   0.661   0.5168  
## X.8           1.0106     1.3386   0.755   0.4600  
## X.9          -0.9203     1.4162  -0.650   0.5240  
## X.10          0.7917     1.3699   0.578   0.5705  
## X.11         -1.8974     1.3784  -1.377   0.1855  
## X.12         -1.8587     1.3906  -1.337   0.1980  
## X.13         -1.6866     1.3178  -1.280   0.2168  
## X.14          0.2098     1.0729   0.196   0.8472  
## Y.1           0.2744     0.2204   1.245   0.2290  
## Y.2           0.3912     0.2061   1.898   0.0739 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 42.69 on 18 degrees of freedom
## Multiple R-squared:  0.7456, Adjusted R-squared:  0.5054 
## F-statistic: 3.104 on 17 and 18 DF,  p-value: 0.01092

Peramalan dan Akurasi

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

## $forecasts
##  [1] 403.7630 313.0896 214.4709 283.5259 313.5155 190.9157 181.9814 234.3299
##  [9] 306.3070 239.6303 297.7806 168.7514 163.7663 171.0712
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 14)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

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

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

## [1] 0.1993977

#akurasi data training
GoF(model.ardl)

##             n     MAE       MPE     MAPE     sMAPE      MASE      MSE     MRAE
## model.ardl 36 24.2383 -2.431347 2.535382 0.1555511 0.6335628 911.4148 2.154835
##                GMRAE
## model.ardl 0.6506468

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. Fungsi dynlm secara umum adalah sebagai berikut.

dynlm(formula, data, subset, weights, na.action, method = "qr",
  model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE,
  contrasts = NULL, offset, start = NULL, end = NULL, ...)

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=3
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 = 50
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt), data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -148.885  -53.000   -0.522   37.637  179.156 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  111.937     27.943   4.006 0.000224 ***
## Xt             7.259      1.585   4.579 3.56e-05 ***
## L(Xt)          2.000      1.324   1.510 0.137860    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 73.12 on 46 degrees of freedom
## Multiple R-squared:  0.4314, Adjusted R-squared:  0.4067 
## F-statistic: 17.45 on 2 and 46 DF,  p-value: 2.294e-06

summary(cons_lm2)

## 
## Time series regression with "ts" data:
## Start = 2, End = 50
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -178.43  -28.82   11.38   28.39   88.94 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.99626   24.15139   1.118 0.269459    
## Xt           4.90093    1.17074   4.186 0.000127 ***
## L(Yt)        0.59311    0.08929   6.643 3.13e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 53.51 on 46 degrees of freedom
## Multiple R-squared:  0.6954, Adjusted R-squared:  0.6822 
## F-statistic: 52.51 on 2 and 46 DF,  p-value: 1.336e-12

summary(cons_lm3)

## 
## Time series regression with "ts" data:
## Start = 2, End = 50
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Yt), data = train.ts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -169.13  -32.02   11.48   30.98   93.68 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.1359    23.7405   1.185    0.242    
## Xt            5.3093     1.1774   4.509 4.62e-05 ***
## L(Xt)        -1.8098     1.1123  -1.627    0.111    
## L(Yt)         0.6792     0.1025   6.629 3.62e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 52.58 on 45 degrees of freedom
## Multiple R-squared:  0.7123, Adjusted R-squared:  0.6931 
## F-statistic: 37.14 on 3 and 45 DF,  p-value: 3.098e-12

summary(cons_lm4)

## 
## Time series regression with "ts" data:
## Start = 3, End = 50
## 
## Call:
## dynlm(formula = Yt ~ Xt + L(Xt) + L(Xt, 2), data = train.ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -148.741  -56.489    5.546   39.075  160.020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  96.6542    32.0225   3.018  0.00422 ** 
## Xt            7.6265     1.6093   4.739 2.27e-05 ***
## L(Xt)         3.4174     1.6187   2.111  0.04048 *  
## L(Xt, 2)     -0.6756     1.3417  -0.504  0.61709    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 72.7 on 44 degrees of freedom
## Multiple R-squared:  0.4529, Adjusted R-squared:  0.4156 
## F-statistic: 12.14 on 3 and 44 DF,  p-value: 6.427e-06

SSE

deviance(cons_lm1)

## [1] 245908.1

deviance(cons_lm2)

## [1] 131732.9

deviance(cons_lm3)

## [1] 124413.6

deviance(cons_lm4)

## [1] 232552.4

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     45 -1 43.9442 3.624e-08 ***
## M2 vs. ME     45 -1  2.6474    0.1107    
## ---
## 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.61011, p-value = 3.812e-09
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(cons_lm2)

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

dwtest(cons_lm3)

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

dwtest(cons_lm4)

## 
##  Durbin-Watson test
## 
## data:  cons_lm4
## DW = 0.62149, p-value = 1.093e-08
## alternative hypothesis: true autocorrelation is greater than 0

Heterogenitas

bptest(cons_lm1)

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

bptest(cons_lm2)

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

bptest(cons_lm3)

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

bptest(cons_lm4)

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

Kenormalan

shapiro.test(residuals(cons_lm1))

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

shapiro.test(residuals(cons_lm2))

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

shapiro.test(residuals(cons_lm3))

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

shapiro.test(residuals(cons_lm4))

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

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.4073759
## DLM 1          0.1596735
## DLM 2          0.2722425
## Autoregressive 0.1993977

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

Plot

par(mfrow=c(1,1))
plot(test$Xt, test$Yt, type="b", col="black")
points(test$Xt, fore.ardl$forecasts,col="green")
lines(test$Xt, fore.ardl$forecasts,col="green")
legend("topleft",c("aktual", "autoregressive"), lty=1, col=c("black","green"), cex=0.8)

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

Pertemuan 3 - Regresi dengan Peubah Lag

Kheni Hikmah Lestari

2023-08-28

Packages

Impor Data

Pembagian Data

Model Koyck

Pemodelan

Peramalan dan Akurasi Model Koyck

Regression with Distributed Lag

Pemodelan (Lag=2)```{r}

Peramalan dan Akurasi

Lag Optimum

Model Autoregressive Distributed Lag (ARDL)

Pemodelan

Lag Optimum

Peramalan dan Akurasi

Pemodelan DLM & ARDL dengan Library `dynlm`

Ringkasan Model

SSE

Uji Diagnostik

Autokorelasi

Heterogenitas

Kenormalan

Perbandingan Model

Plot

Pertemuan 3 - Regresi dengan Peubah Lag

Kheni Hikmah Lestari

2023-08-28

Packages

Impor Data

Pembagian Data

Model Koyck

Pemodelan

Peramalan dan Akurasi Model Koyck

Regression with Distributed Lag

Pemodelan (Lag=2)```{r}

Peramalan dan Akurasi

Lag Optimum

Model Autoregressive Distributed Lag (ARDL)

Pemodelan

Lag Optimum

Peramalan dan Akurasi

Pemodelan DLM & ARDL dengan Library dynlm

Ringkasan Model

SSE

Uji Diagnostik

Autokorelasi

Heterogenitas

Kenormalan

Perbandingan Model

Plot

Pemodelan DLM & ARDL dengan Library `dynlm`