Library yang digunakan

## 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: AQI Peubah penjelas: o3

data <- rio::import("https://raw.githubusercontent.com/mrnabilnaufal07/mpdw/main/Pertemuan%202/NewDelhi_Air_quality.csv")
Yt = data$AQI
Xt = data$o3
data = data.frame(Yt,Xt)
str(data)
## 'data.frame':    72 obs. of  2 variables:
##  $ Yt: num  30.2 28.2 26.6 25 26 26 26 24 23 24 ...
##  $ Xt: num  55.8 54.9 54.6 55.1 55.8 ...
data
##      Yt       Xt
## 1  30.2 55.78995
## 2  28.2 54.93164
## 3  26.6 54.64554
## 4  25.0 55.07469
## 5  26.0 55.78995
## 6  26.0 56.50520
## 7  26.0 55.78995
## 8  24.0 52.92892
## 9  23.0 50.06790
## 10 24.0 51.49841
## 11 25.0 53.64418
## 12 26.0 55.78995
## 13 27.0 59.36623
## 14 29.0 62.22725
## 15 29.0 62.94250
## 16 29.0 62.94250
## 17 29.0 62.22725
## 18 28.0 60.79674
## 19 27.0 59.36623
## 20 27.0 58.65097
## 21 27.0 58.65097
## 22 27.0 57.93572
## 23 27.0 57.93572
## 24 27.0 58.65097
## 25 27.0 59.36623
## 26 28.0 61.51199
## 27 29.0 63.65776
## 28 31.0 66.51878
## 29 31.0 67.94930
## 30 32.0 68.66455
## 31 32.0 68.66455
## 32 31.0 67.23404
## 33 31.0 66.51878
## 34 30.0 65.80353
## 35 30.0 64.37302
## 36 29.0 62.22725
## 37 27.0 59.36623
## 38 26.0 56.50520
## 39 25.0 55.07469
## 40 25.0 53.64418
## 41 24.0 52.92892
## 42 24.0 52.92892
## 43 25.0 53.64418
## 44 25.0 53.64418
## 45 25.0 54.35944
## 46 26.0 55.78995
## 47 26.0 57.22046
## 48 27.0 57.93572
## 49 27.0 58.65097
## 50 27.0 58.65097
## 51 27.0 59.36623
## 52 28.0 60.08148
## 53 28.0 60.08148
## 54 27.0 59.36623
## 55 27.0 58.65097
## 56 26.0 57.22046
## 57 25.0 54.35944
## 58 23.0 50.78316
## 59 22.0 47.92213
## 60 21.0 45.41874
## 61 20.0 43.98823
## 62 19.0 41.84246
## 63 19.0 41.48483
## 64 20.0 43.27297
## 65 21.0 45.41874
## 66 21.0 45.41874
## 67 21.0 45.41874
## 68 22.0 47.20688
## 69 24.0 52.21367
## 70 26.0 56.50520
## 71 27.0 59.36623
## 72 28.0 60.79674

Pembagian Data

#Split data dengan proporsi: 80% training dan 20% testing
train<-data[1:57,] # 80% data pertama
test<-data[58:72,] # 20% data sisanya

Mengubah format data menjadi time series

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

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 
## -1.127670 -0.301407  0.001002  0.262601  1.103471 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.54798    0.92878   0.590    0.558    
## Y.1          0.41955    0.08223   5.102 4.63e-06 ***
## X.t          0.25830    0.04107   6.289 6.25e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5053 on 53 degrees of freedom
## Multiple R-Squared: 0.9463,  Adjusted R-squared: 0.9443 
## Wald test: 423.7 on 2 and 53 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                             alpha      beta       phi
## Geometric coefficients:  0.944071 0.2582957 0.4195528

Dari hasil tersebut, didapat bahwa peubah \(x_t\) dan \(y_{t-1}\) memiliki nilai \(P-Value<0.05\). Hal ini menunjukkan bahwa peubah \(x_t\) dan \(y_{t-1}\) berpengaruh signifikan terhadap \(y\). Artinya, menurut model Koyck, nilai AQI saat ini dipengaruhi oleh kandungan \(O_3\) pada saat ini, serta nilai AQI satu hari sebelumnya. Adapun model keseluruhannya adalah sebagai berikut

\[ \hat{Y_t}=0.54798+0.25830X_t+0.41955Y_{t-1} \]

Peramalan dan Akurasi Model Koyck

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

#h =15, merupakan 15 periode yang akan diprediksi selanjutnya
fore.koyck <- forecast(model = model.koyck, x=test$Xt, h=15)
fore.koyck
## $forecasts
##  [1] 24.15387 23.05989 21.95429 21.12094 20.21706 19.74546 20.00947 20.67448
##  [9] 20.95348 21.07054 21.58152 23.08914 24.83015 26.29958 27.28558
## 
## $call
## forecast.koyckDlm(model = model.koyck, x = test$Xt, h = 15)
## 
## 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.03187326
## 
## $`MAPE Training.MAPE`
## [1] 0.01435646

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

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.

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.35051 22.97129 47.02321 14499.82 0.49589  0.98302  0.354745

Diperoleh lag optimum untuk peubah \(O_3\) adalah 10 hari sebelumnya. Selanjutnya dilakukan pemodelan kembali dengan \(q=10\)

Peramalan dan Akurasi Model DLM (Distributed Lag Model)

#model dlm dengan lag optimum
model.dlm <- dlm(x = train$Xt,y = train$Yt , q = 10) #nilai q diganti dengan lag optimum yang diperoleh pada langkah sebelumnya
summary(model.dlm)
## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5204 -0.1663  0.0356  0.1624  0.5130 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.294780   0.931786  -0.316  0.75361    
## x.t          0.405638   0.063702   6.368 2.54e-07 ***
## x.1          0.159000   0.133129   1.194  0.24038    
## x.2         -0.052741   0.126788  -0.416  0.67997    
## x.3         -0.148621   0.121282  -1.225  0.22860    
## x.4          0.143544   0.124404   1.154  0.25638    
## x.5         -0.007047   0.123348  -0.057  0.95476    
## x.6         -0.127683   0.123191  -1.036  0.30709    
## x.7          0.287823   0.121937   2.360  0.02396 *  
## x.8         -0.399437   0.121443  -3.289  0.00230 ** 
## x.9          0.304147   0.103813   2.930  0.00594 ** 
## x.10        -0.098958   0.048025  -2.061  0.04684 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2715 on 35 degrees of freedom
## Multiple R-squared:  0.9871, Adjusted R-squared:  0.983 
## F-statistic: 243.1 on 11 and 35 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##        AIC      BIC
## 1 22.97129 47.02321

Dari hasil tersebut terdapat beberapa peubah yang berpengaruh signifikan terhadap taraf nyata 10% yaitu \(x_t\) , \(x_{t-7}\) , \(x_{t-8}\) , \(x_{t-9}\), dan \(x_{t-10}\). Artinya, menurut model DLM dengan \(q=10\), nilai AQI saat ini dipengaruhi oleh kandungan \(O_3\) pada saat ini, 7 hari sebelumnya, 8 hari sebelumnya, 9 hari sebelumnya, dan 10 hari sebelumnya. Adapun keseluruhan model yang terbentuk adalah

\[ \hat{Y_t}=-0.294780 + 0.405638X_t+...-0.098958X_{t-10} \]

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

#peramalan dan akurasi
fore.dlm <- forecast(model = model.dlm, x=test$Xt, h=15)

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

#akurasi data training
mape.dlm.train = GoF(model.dlm)["MAPE"]

c("MAPE Testing"=mape.dlm,"MAPE Training"=mape.dlm.train)
## $`MAPE Testing`
## [1] 0.01274404
## 
## $`MAPE Training.MAPE`
## [1] 0.007238929

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

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 ARDL

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.

Lag Optimum untuk ARDL

#penentuan lag optimum
model.ardl.opt <- ardlBoundOrders(data = data.frame(data), ic = "AIC", 
                                  formula = Yt ~ Xt )
min_p=c()
for(i in 1:15){
  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         2        13 16.55044

Dari tabel di atas, dapat terlihat bahwa nilai AIC terendah didapat ketika \(p=13\) dan \(q=2\), yaitu sebesar 16.55044. Artinya, lag optimum untuk peubah \(X_t\) \((O_3)\) adalah 13 hari sebelumnua sedangkan lag optimum untuk peubah \(Y_t\) \((AQI)\) adalah 2 hari. Selanjutnya nilai ini akan dimasukkan ke dalam proses pembentukan model ardl.

model.ardl <- ardlDlm(x = train$Xt, y = train$Yt, p = 13 , q = 2)
summary(model.ardl)
## 
## Time series regression with "ts" data:
## Start = 14, End = 57
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45939 -0.19244  0.04237  0.18424  0.45823 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.77136    1.12356  -0.687   0.4982    
## X.t          0.39407    0.07142   5.517 7.61e-06 ***
## X.1          0.20786    0.18280   1.137   0.2655    
## X.2         -0.02966    0.21795  -0.136   0.8928    
## X.3         -0.13619    0.16006  -0.851   0.4023    
## X.4          0.05226    0.16212   0.322   0.7497    
## X.5          0.04957    0.14512   0.342   0.7353    
## X.6         -0.13909    0.13618  -1.021   0.3162    
## X.7          0.27019    0.13651   1.979   0.0581 .  
## X.8         -0.31492    0.14669  -2.147   0.0409 *  
## X.9          0.15294    0.16112   0.949   0.3509    
## X.10         0.09506    0.16928   0.562   0.5790    
## X.11        -0.16098    0.15567  -1.034   0.3103    
## X.12         0.04585    0.12175   0.377   0.7094    
## X.13         0.01492    0.05471   0.273   0.7871    
## Y.1         -0.12772    0.19313  -0.661   0.5140    
## Y.2          0.06633    0.21010   0.316   0.7546    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2868 on 27 degrees of freedom
## Multiple R-squared:  0.9884, Adjusted R-squared:  0.9814 
## F-statistic: 143.2 on 16 and 27 DF,  p-value: < 2.2e-16

Terdapat 2 peubah yang berpengaruh signifikan terhadap nilai \(AQI\) pada selang kepercayaan 95% yaitu \(X_t\) dan \(X_{t-8}\). Artinya, menurut model ARDL dengan \(p=13\) dan \(q=2\), nilai AQI saat ini dipengaruhi oleh kandungan \(O_3\) pada saat ini, serta 8 hari sebelumnya. Model keseluruhannya adalah sebagai berikut:

\[ \hat{Y}=-0.77136+ 0.39407X_t+0.20786 X_{t-1}+...-0.12772Y_{t-1}+0.06633Y_{t-2} \]

Peramalan dan Akurasi Model ARDL

fore.ardl <- forecast(model = model.ardl, x=test$Xt, h=15)
fore.ardl
## $forecasts
##  [1] 23.31655 21.80904 20.61682 20.06357 19.24186 18.82466 19.48054 20.86204
##  [9] 21.19795 20.60591 21.49901 23.57998 26.28620 27.75765 28.02269
## 
## $call
## forecast.ardlDlm(model = model.ardl, x = test$Xt, h = 15)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

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

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

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

c("MAPE Testing"=mape.ardl,"MAPE Training"=mape.ardl.train)
## $`MAPE Testing`
## [1] 0.01378089
## 
## $`MAPE Training.MAPE`
## [1] 0.007069452

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

Perbandingan Model

Akurasi

akurasi <- matrix(c(mape.koyck, mape.dlm, mape.ardl))
row.names(akurasi)<- c("Koyck","DLM","Autoregressive")
colnames(akurasi) <- c("MAPE")
akurasi
##                      MAPE
## Koyck          0.03187326
## DLM            0.01274404
## Autoregressive 0.01378089

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

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.ardl$forecasts,col="green")
lines(test$Xt, fore.ardl$forecasts,col="green")
legend("topleft",c("aktual", "koyck","DLM", "autoregressive"), lty=1, col=c("black","red","blue","green"), cex=0.8)

Berdasarkan plot tersebut, terlihat bahwa plot yang paling mendekati data aktualnya adalah Distributed Lag Model (DLM), sehingga dapat disimpulkan model terbaik dalam hal ini adalah model DLM

Kesimpulan

Dari ketiga model yang dicobakan terhadap pengaruh kadar \(O_3\) terhadap \(AQI\) di kota New Delhi, diperoleh kesimpulan bahwa Distributed Lag Model (DLM) adalah yang paling baik dalam peramalan data tersebut.

Daftar Pustaka

Gujarati, D. N. (2004). Basic Econometrics.Fourth Edition. New York: The McGraw Hill

https://www.kaggle.com/datasets/anuragbantu/new-delhi-air-quality