2.1.1 Plot residual vs order

Uji eksploratif dapat dilakukan dengan melihat beberapa chart atau plot dari sisaan. Syntax yang digunakan untuk mendapat hal tersebut adalah sebagai berikut:

resi <- residuals(model)

par(mfrow = c(2,2))
plot(model, which=1:2)
hist(resi, col = "steelblue", main = "Histogram Sisaan")
plot(resi, type="o", 
     ylab = "Sisaan", xlab = "Order", main = "Sisaan vs Order")
abline(h = 0, col='red')

Terdapat 4 output yang dihasilkan dari syntax di atas, yaitu Residuals vs Fitted, Normal Q-Q Plot, Histogram sisaan, dan sisaan vs order. Keempat chart di atas diperlukan untuk memeriksa asumsi regresi. Dalam hal ini, fokus kita adalah pemeriksaan kebebasan galat yang terdapat pada plot keempat, yaitu plot Sisaan vs Order. Dari plot tersebut, terdapat pola naik turun yang mengindikasikan bahwa sisaan dari galatnya tidak saling bebas atau terdapat autokorelasi.

2.1.2 ACF dan PACF

Selain plot Sisaan vs Order, dapat dilihat Autocorrelation Function (ACF) dan plot Partial Autocorrelation Function (PACF). Pengujian kedua plot ini menggunakan lag sebagai pembandingnya. Lag adalah beda urutan suatu sisaan dengan sisaan sebelumnya. Misal, lag 1 berarti tiap sisaan dibandingkan dengan sisaan satu observasi sebelumnya. Jika garis vertikal di tiap lag lebih tinggi dari garis biru horizontal, maka dianggap terjadi autokorelasi di lag tersebut (kecuali lag 0 di plot ACF).

acf(resi)

pacf(resi)

Baik dari plot ACF maupun plot PACF, tidak ada garis vertikal di lag tertentu yang melebihi tinggi garis biru horizontal. Artinya, menurut kedua plot ini, tidak terdapat autokorelasi pada model.

2.2.1 Durbin Watson Test

Salah satu uji formal yang biasa digunakan untuk menguji kebebasan galat adalah Uji Durbin Watson. Uji ini hanya dapat memperhatikan ketidakbebasan pada lag 1.

Adapun asumsi-asumsi yang perlu dipenuhi adalah:

1. Terdapat intersep pada model regresi
2. Seluruh peubah penjelas bersifat tetap (fixed) pada penarikan contoh berulang
3. Galat mengikuti skema Autoregressive (AR) ordo ke-1 : \[u_t=\rho u_{t-1}+v_t\]dengan \(\rho\) adalah koefisien autokorelasi yang bernilai -1 s.d 1 dan \(v_t\) adalah komponen sisaan yang memenuhi asumsi kebebasan dan kehomogenan.
4. Galat menyebar normal
5. Lag dari peubah respon tidak disertakan sebagai peubah penjelas dalam model

Jika \(e_k\) adalah residual yang diurutkan berdasarkan waktu, statistik Durbin Watson adalah:

\[ d=\frac{\sum_{t=1}^{n}(e_t-e_{t-1})^2}{\sum_{t=1}^{n}e_t^2}; 0\le d\le4 \]

Interpretasi dari nilai statistik di atas jika dilakukan secara manual yaitu:

• Jika \(d\) mendekati 0, maka semakin besar kemungkinan adanya autokorelasi positif

• Jika \(d\) mendekati 2, maka belum cukup bukti adanya autokorelasi

• Jika \(d\) mendekati 0, maka semakin besar kemungkinan adanya autokorelasi negatif

lmtest::dwtest(model, alternative = 'two.sided')

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.87546, p-value = 0.007054
## alternative hypothesis: true autocorrelation is not 0

2.2.2 Breusch-Godfrey Test

Seperti dijelaskan sebelumnya, uji Durbin Watson hanya dapat menguji autokorelasi lag 1. Padahal, terlihat di plot ACF dan PACF bahwa autokorelasi dapat muncul di lag lainnya. Salah satu alternatif yang dapat dilakukan adalah dengan Breusch-Gofrey Test.

bgtest(IPM_Bali ~ Tahun, data=df, order=1)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  IPM_Bali ~ Tahun
## LM test = 1.8283, df = 1, p-value = 0.1763

2.2.3 Runs Test

Runs didefinisikan sebagai serangkaian nilai yang meningkat atau serangkaian nilai yang menurun. Misalnya ada 20 data, jika nilai data naik dari nilai sebelumnya dinotasikan positif (+), dan sebaliknya dinotasikan negatif (-). Jika proses tersebut menghasilkan + + + + − − − + + + − − + + + + - - - -, ada 6 run, 3 positif dan 3 negatif. Sebuah run merupakan kumpulan tanda positif atau negatif, dinotasikan sebagai \(R\). Dengan statistik Run berupa nilai Z sebagai berikut:

\[ Z=\frac{R-\bar{R}}{s_R} \]

Run Test akan menolak \(H_0\) jika \[|Z|>Z_{1-\alpha/2}\]

lawstat::runs.test(resid(model), alternative = 'two.sided')

## 
##  Runs Test - Two sided
## 
## data:  resid(model)
## Standardized Runs Statistic = 0, p-value = 1

Penanganan Cochrane-Orcutt:

#Interactive method using to solve first order autocorrelation problems.
modelco <- orcutt::cochrane.orcutt(model,convergence = 7,max.iter = 10^6)
modelco

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = IPM_Bali ~ Tahun, data = df)
## 
##  number of interaction: 234
##  rho 0.73534
## 
## Durbin-Watson statistic 
## (original):    0.87546 , p-value: 3.527e-03
## (transformed): 2.13019 , p-value: 4.363e-01
##  
##  coefficients: 
## (Intercept)       Tahun 
## -682.274737    0.375122

# Mencari nilai optimum bagi rho
rho <- modelco$rho
rho

## [1] 0.7353396

Telah didapat nilai rho (\(\rho\)) optimum sebesar \(0.7353\). Selanjutnya akan dilakukan transformasi peubah X dan Y menggunakan nilai \(\rho\).

# Membuat variabel baru
y <- df$IPM_Bali
x <- df$Tahun

#transformasi terhadap y dan x
(y.trans <- y[-1]-y[-12]*rho)

##  [1] 19.32270 19.50649 19.42498 19.46937 19.97259 19.77167 20.14224 20.13427
##  [9] 20.39866 20.07010 20.17186

(x.trans <- x[-1]-x[-12]*rho)

##  [1] 532.9675 533.2321 533.4968 533.7615 534.0261 534.2908 534.5554 534.8201
##  [9] 535.0848 535.3494 535.6141

#model baru
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho) 

## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18524 -0.15403 -0.03408  0.13616  0.24765 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -180.57113   33.81062  -5.341 0.000468 ***
## x.trans        0.37512    0.06328   5.928 0.000221 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1757 on 9 degrees of freedom
## Multiple R-squared:  0.7961, Adjusted R-squared:  0.7734 
## F-statistic: 35.14 on 1 and 9 DF,  p-value: 0.0002213

Didapat model transformasi dengan intersep sebesar \(-180.571\) dan koefisien x.trans sebesar \(0.375\).

lmtest::dwtest(modelcorho,alternative = 'two.sided') 

## 
##  Durbin-Watson test
## 
## data:  modelcorho
## DW = 2.1302, p-value = 0.8726
## alternative hypothesis: true autocorrelation is not 0

Berdasarkan hasil di atas, didapat nilai \(P-Value=0.8726>0.05\), sehingga \(H_0\) tidak ditolak. Artinya, tidak ada autokorelasi atau terbukti bahwa autokorelasi telah berhasil ditangani dengan metode Cochrane-Orcutt. Selanjutnya akan dilakukan transformasi balik untuk mendapatkan model regresi deret waktu terbaru.

(b0 <- modelcorho$coefficients[1]/(1-rho))

## (Intercept) 
##   -682.2747

(b1 <- modelcorho$coefficients[2])

##   x.trans 
## 0.3751221

Setelah dilakukan transformasi balik, didapat model regresi yang baru yaitu:

\[ \hat{y_t}=-682.275+0.375x_t \]

Penanganan Hildreth-lu:

# Hildreth-Lu
hildreth.lu.func<- function(r, model){
  x <- model.matrix(model)[,-1]
  y <- model.response(model.frame(model))
  n <- length(y)
  t <- 2:n
  y <- y[t]-r*y[t-1]
  x <- x[t]-r*x[t-1]
  
  return(lm(y~x))
}

#mencari rho yang meminimumkan SSE (iteratif)
r <- c(seq(0.1,0.8, by= 0.1), seq(0.9,0.99, by= 0.01))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0,40 memiliki SSE Terkecil

#grafik rho dan SSE
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

Pada iterasi pertama dapat terlihat bahwa \(\rho\) optimum yang menghasilkan SES terkecil berada di sekitar \(0.7\). Untuk hasil yang lebih akurat, selang dapat diperkecil.

#rho optimal di sekitar 0.7
r <- seq(0.68, 0.75, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4)

#grafik SSE optimum
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

Dari grafik di atas, dapat terlihat bahwa nilai \(\rho\) optimum berada di \(\rho=0.74\). Sehingga nilai ini akan digunakan untuk merumuskan model terbaik pada tahap selanjutnya:

# Model terbaik
modelhl <- hildreth.lu.func(0.74, model)
summary(modelhl)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18566 -0.15366 -0.03286  0.13587  0.24747 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -175.77964   33.81212  -5.199 0.000565 ***
## x              0.37206    0.06442   5.776 0.000267 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1757 on 9 degrees of freedom
## Multiple R-squared:  0.7875, Adjusted R-squared:  0.7639 
## F-statistic: 33.36 on 1 and 9 DF,  p-value: 0.0002675

# Deteksi autokorelasi
lmtest::dwtest(modelhl)

## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 2.1385, p-value = 0.4421
## alternative hypothesis: true autocorrelation is greater than 0

Pada uji Durbin-Watson di atas, dapat terlihat bahwa nilai \(P-Value=0.4421>0.05\)sehingga \(H_0\) tidak ditolak. Artinya, tidak ada autokorelasi atau terbukti bahwa autokorelasi telah berhasil ditangani dengan metode Hidreth-lu. Selanjutnya akan dilakukan transformasi balik untuk mendapatkan model regresi deret waktu terbaru.

# Transformasi Balik
cat("IPM = ", coef(modelhl)[1]/(1-0.74), " + ", coef(modelhl)[2]," tahun", sep = "")

## IPM = -676.0755 + 0.3720559 tahun

Setelah dilakukan transformasi balik, didapat model regresi yang baru yaitu:

\[ \hat{y_t}=-676.0755+0.372x_t \]