Penanganan Autokorelasi pada Regresi Linear Sederhana untuk Data Deret Waktu

Salah satu asumsi pada regresi linear sederhana yaitu kebebasan antar sisaan. Pada analisis regresi sangat mungkin untuk mengalami autokorelasi, yaitu ketika amatan tertentu dipengaruhi oleh amatan sebelumnya yang menyebabkan tidak saling bebasnya sisaan antar amatan.

Angka IPM adalah salah satu contoh data deret waktu, IPM di Indonesia sejak tahun 2014 disajikan secara tahunan. Pengertian Indeks Pembangunan Manusia (IPM) sendiri sebagaimana yang dirilis oleh United Nations Development Programme (UNDP) yaitu “merupakan salah satu pendekatan untuk mengukur tingkat keberhasilan pembangunan manusia”. IPM di Indonesia juga turut disajikan untuk setiap Provinsi secara tahunan oleh BPS.

Selanjutnya dengan menggunakan IPM sebagai peubah respon, dan tahun rilis IPM sebagai peubah bebas, akan dibangun model regresi linear sederhana. Model regresi ini selanjutnya akan dilakukan deteksi autokorelasi serta dilakukan penanganannya.

package

library(dplyr)
library(readxl)
library(forecast)
library(TTR)
library(tseries)
library(lmtest) #uji-Durbin Watson
library(orcutt) #Cochrane-Orcutt
library(HoRM)#Hildreth Lu

DATA

Data yang digunakan pada pemodelan regresi linear sederhana untuk data deret waktu ini adalah Angka IPM Provinsi Kepulauan Bangka Belitung untuk tahun 2010-2021. Dengan jumlah amatan sebanyak 12 amatan, Angka IPM sebagai peubah respon (y) dan tahun sebagai peubah bebas (x).

#membaca data
data <- read_excel("C:/Users/HP/Downloads/Data Deret Waktu IPM untuk Regresi.xlsx")
head (data)
## # A tibble: 6 x 2
##   Tahun   IPM
##   <dbl> <dbl>
## 1  2010  66.0
## 2  2011  66.6
## 3  2012  67.2
## 4  2013  67.9
## 5  2014  68.3
## 6  2015  69.0
#menginisialisasi peubah
x <- data$Tahun
y <- data$IPM

EKSPLORASI DATA

Eksplorasi data dilakukan melalui grafik dan juga korelasi antara x dengan y. Eksplorasi ini dilakukan untuk mengetahui hubungan antara x (Tahun) dengan y (IPM)

#Scatter plot X vs Y
plot(x,y,pch = 20, col = "blue", main = "Scatter Plot Tahun vs IPM",
     ylab = "Tahun", xlab = "IPM")

#korelasi x dan y
cor(x,y)
## [1] 0.995183

Dapat terlihat bahwa sejak tahun 2010 hingga 2019, IPM Provinsi Kepulauan Bangka Belitung menunjukkan nilai yang terus meningkat seiring bertambahnya tahun. Namun di tahun 2020, kenaikan IPM ini mengalami sedikit penurunan (yang ditunjukkan oleh titik amatan yang cenderung membuat pola garis sebelumnya menjadi lebih landai). Secara keseluruhan dapat disimpulkan bahwa semakin bertambahnya tahun IPM Provinsi Kepulauan Bangka Belitung juga meningkat. Hal ini selaras dengan nilai korelasinya yang tinggi yaitu 0.995.

MODEL REGRESI DERET WAKTU

model <- lm(y~x, data)
summary(model)
## 
## Call:
## lm(formula = y ~ x, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42205 -0.10830 -0.00455  0.14593  0.26718 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.018e+03  3.388e+01  -30.06 3.88e-11 ***
## x            5.396e-01  1.681e-02   32.10 2.03e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.201 on 10 degrees of freedom
## Multiple R-squared:  0.9904, Adjusted R-squared:  0.9894 
## F-statistic:  1030 on 1 and 10 DF,  p-value: 2.026e-11

dengan model ini, selanjutnya akan dilakukan deteksi autokorelasi dalam model dengan metode grafik dan pengujian statistik

DETEKSI AUTOKORELASI

GRAFIK

Pendeteksian autokorlasi dalam model regresi dapat dilakukan melalui residual plot dan ACF dan PACF

1. Residual Plot

Residual Plot ini meliputi normal Q-Q plot dan Histogram dari Sisaan untuk mengidentifikasi kenormalan data, serta Sisaan vs Fitted Values dan Sisaan vs Order untuk mengidentifikasi adanya autokorelasi.

#sisaan dan fitted value
res <- residuals(model)
fit <- predict(model)

#Diagnostik dengan eksploratif
par(mfrow = c(2,2))

qqnorm(res)
qqline(res, col = "steelblue", lwd = 2)

plot(fit, res, col = "steelblue", pch = 20, xlab = "Fitted Values", 
     ylab = "Sisaan", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)

hist(res, col = "steelblue", main = "Histogram dari Sisaan")

plot(seq(1,12,1), res, col = "steelblue", pch = 20, 
     xlab = "Order", ylab = "Sisaan", main = "Sisaan vs Order")
lines(seq(1,12,1), res, col = "red")
abline(a = 0, b = 0, lwd = 2)

  • Melalui Normal Q-Q Plot dan juga HIstogram dari Sisaan diindikasikan bahwa sisaan tidaklah menyebar normal, hal ini terlihat melalui Normal Q-Q Plot dimana terdapat beberapa amatan yang berada lumayan jauh dari garis Normal, melalui histogram dari sisaan yang cenderung tidak menunjukkan bentuk simetris serta terlihat cenderung menjulur ke kiri.

  • Melalui Plot Sisaan vs Fitted Values, terlihat bahwa sisaan tidaklah menyebar secara acak yang ditunjukkan oleh adanya pola kenaikan dan penurunan. Hal ini ditunjukkan pada Plot Sisaan vs Order, maka terlihat terdapat pola naik untuk sisaan amatan 1 hingga amatan 4, pola turun untuk sisaan amatan 10 hinggan amatan 12. Hal ini menunjukkan bahwa terindikasi terdapat autokorelasi pada model regresi ini.

2. ACF dan PCF

PLot ACF (autocorrelation function) dan Plot PACF (partial autocorrelation function) dapat menunjukkan indikasi autokorelasi dalam model jika terdapat lag yang nilai ACF ataupun PACF (yang ditunjukkan garis hitam) yang melewati batas signifikansi (garis putus-putus biru)

par(mfrow = c(2,1))
acf(res)
pacf(res)

Baik Plot ACF ataupun Plot PACF, terlihat tidak terdapat lag yang berada diluar garis putus-putus biru (tidak ada nilai ACF ataupun PACF yang signifikan). Hal ini mengindikasikan bahwa tidak terdapat autokorelasi (sisaan saling bebas) pada model regresi ini.

UJI STATISTIK

Selain melalui grafik, pendeteksian autokorelasi dalam model regresi dapat dilakukan melalu pengujian statistik yaitu:

  • Durbin Watson Test
  • Breusch-Godfrey Test
  • Run’s Test.

Pengujian dengan uji statistik ini dilakukan dengan hipotesis sebagai berikut:
H0 : tidak terdapat autokorelasi
H1 : terdapat autokorelasi

#Durbin Watson Test
lmtest::dwtest(model, alternative = 'two.sided')
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.0105, p-value = 0.01864
## alternative hypothesis: true autocorrelation is not 0
#Breusch-Godfrey Test
lmtest::bgtest(y ~ x, data=data, order=1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  y ~ x
## LM test = 1.8237, df = 1, p-value = 0.1769
#Run's Test
lawstat::runs.test(resid(model), alternative = 'two.sided')
## 
##  Runs Test - Two sided
## 
## data:  resid(model)
## Standardized Runs Statistic = -1.2111, p-value = 0.2259

Diantara ketiga pengujian statistik untuk mengidentifikasi autokorelasi dalam model, hasil pengujian Durbin-Watson menunjukkan p-value < 5%, tolak H0, sehingga cukup bukti untuk menyatakan terdapat autokorelasi dalam model pada taraf 5%. Sedangkan untuk Breush-Godfrey Test dan Run’s Test, hasil pengujian keduanya menunjukkan p-value > 5%, gagal tolak H0, sehingga tidak terdapat cukup bukti untuk menyatakan terdapat autokorelasi dalam model pada taraf 5%

PENANGANAN AUTOKORELASI

Penanganan autokorelasi dalam model regresi dapat dilakukan dengan metode:

  • Cochrane-Orcutt
  • Hildreth-lu

1. Cochrane-Orcutt

#Interactive method using to solve first order autocorrelation problems. 

modelco <- cochrane.orcutt(model, convergence = 7, max.iter = 1000)
modelco
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x, data = data)
## 
##  number of interaction: 213
##  rho 0.754277
## 
## Durbin-Watson statistic 
## (original):    1.01045 , p-value: 9.321e-03
## (transformed): 1.92933 , p-value: 3.032e-01
##  
##  coefficients: 
## (Intercept)           x 
## -805.787582    0.434258
#rho optimum
rho <- modelco$rho

rho-optimum hasil iterasi adalah 0.754277, nilai ini selanjutnya akan digunakan dalam transformasi balik kedalam bentuk model linier hasil penanganan autokorelasi (yaitu modelcorho)

#transformasi terhadap y dan x
y.trans <- y[-1]-y[-12]*rho
x.trans <- x[-1]-x[-12]*rho
#model hasil penanganan dengan metode cochrane-orcutt
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho) 
## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.20061 -0.15420 -0.02557  0.15479  0.22623 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -198.00082   34.15317  -5.797  0.00026 ***
## x.trans        0.43426    0.06884   6.308  0.00014 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1774 on 9 degrees of freedom
## Multiple R-squared:  0.8156, Adjusted R-squared:  0.7951 
## F-statistic:  39.8 on 1 and 9 DF,  p-value: 0.0001396
#pengujian kembali
lmtest::dwtest(modelcorho,alternative = 'two.sided')
## 
##  Durbin-Watson test
## 
## data:  modelcorho
## DW = 1.9293, p-value = 0.6064
## alternative hypothesis: true autocorrelation is not 0

Untuk memeriksa apakah autokorelasi yang ada pada model sebelumnya sudah tertangani, maka dilakukan Durbin-Watson test (uji statistik yang sebelumnya menunjukkan autokorelasi pada model awal) pada modelcorho. Dengan Hipotesis yaitu:
H0 : tidak terdapat autokorelasi
H1 : terdapat autokorelasi
didapatkan p-value sebesar 0.6064>5%, gagal tolak H0, sehingga tidak terdapat cukup bukti untuk menyatakan terdapat autokorelasi pada modelcorho.

Penanganan autokorelasi dengan metode cochrane-orcutt berhasil dilakukan pada model awal, sehingga tidak terdapat autokorelasi pada modelorcho.

# Transformasi Balik
cat("y = ", modelcorho$coefficients[1]/(1-rho), "+", modelcorho$coefficients[2],"x", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal
## y = -805.7876+0.4342584x

Persamaan regresi setelah dilakukan penanganan autokorelasi pada model yaitu :
y = -805.7876+0.4342584x

2. 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))
}

diperlukan fungsi khusus untuk melakukan iterasi dalam rangka mencari rho optimum menggunakan metode hildreth-lu seperti fungsi diatas.

#mencari rho yang meminimumkan SSE (iteratif)
r <- c(seq(0.1,0.99, by= 0.01))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
head(round(tab, 4))
##    rho    SSE
## 1 0.10 0.3481
## 2 0.11 0.3461
## 3 0.12 0.3442
## 4 0.13 0.3423
## 5 0.14 0.3404
## 6 0.15 0.3386
tail(round(tab, 4))
##     rho    SSE
## 85 0.94 0.2885
## 86 0.95 0.2891
## 87 0.96 0.2897
## 88 0.97 0.2903
## 89 0.98 0.2910
## 90 0.99 0.2917

Pencarian rho optimum dilakukan melalui iterasi fungsi pada semua nilai rho pada rentang 0.1 hingga 0.99, dengan meminimumkan nilai Sum Squares of Error (SSE) dari model yang mungkin terbuat.

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

Melalui grafik diatas, terlihat bahwa nilai SSE optimum yaitu ketika rho berada pada 0.7 < rho < 0.8. Setelah dilakukan pencarian nilai rho dengan SSE minimum, didapatkan nilai rho-optimum untuk model hildreth-lu yaitu sebesar 0.75.

#rho optimum
rho_optimum<-tab[which.min(tab$SSE),1]
rho_optimum
## [1] 0.75

Selanjutnya dengan rho optimum yang sudah ditemukan pada iterasi sebelumnya, dilakukan pemodelan kembali yaitu modelhl sebagai model regresi hasil penanganan autokorelasi dengan metode hildreth-lu.

# Model terbaik
modelhl <- hildreth.lu.func(rho_optimum, model)
summary(modelhl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.20170 -0.15408 -0.02525  0.15459  0.22652 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -202.55680   34.15215  -5.931 0.000220 ***
## x              0.43645    0.06766   6.451 0.000118 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1774 on 9 degrees of freedom
## Multiple R-squared:  0.8222, Adjusted R-squared:  0.8024 
## F-statistic: 41.61 on 1 and 9 DF,  p-value: 0.0001181
# pengujian kembali
lmtest::dwtest(modelhl)
## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 1.9229, p-value = 0.2993
## alternative hypothesis: true autocorrelation is greater than 0

Hasil pengujian kembali modelhl menggunakan uji statistik Durbin-Watson test dengan hipotesis:
H0 : tidak terdapat autokorelasi
H1 : terdapat autokorelasi
didapatkan p-value sebesar 0.2993>5%, gagal tolak H0, sehingga tidak terdapat cukup bukti untuk menyatakan terdapat autokorelasi pada modelhl.

Penanganan autokorelasi dengan metode hildreth-lu berhasil dilakukan pada model awal, sehingga tidak terdapat autokorelasi pada modelhl.

# Transformasi Balik
cat("y = ", coef(modelhl)[1]/(1-rho_optimum), "+", coef(modelhl)[2],"x", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal
## y = -810.2272+0.4364545x

Persamaan regresi setelah dilakukan penanganan autokorelasi pada model yaitu :
y = -810.2272+0.4364545x

KESIMPULAN

  • Pemodelan IPM sebagai peubah respon dan waktu rilis IPM sebagai peubah bebas menunjukkan adanya autokorelasi baik dalam pendeteksian secara grafik ataupun uji statistik (Durbin Watson test).

  • Penanganan yang dilakukan pada model awal berhasil dilakukan baik dengan metode Cochrane-Orcutt dengan rho sebesar 0.7542767 ataupun Hildreth-lu dengan rho sebesar 0.75. Hal ini ditunjukkan dengan tidak signifikannya pengujian autokorelasi dengan Durbin Watson test.

  • Perbandingan model adalah sebagai berikut:
    MODEL AWAL (y = -1018.4506410+0.5396154x)
    MODEL COCHRANE-ORCUTT (y = -805.7876+0.4342584x)
    MODEL HILDRETH-LU (y = -810.2272+0.4364545x)
    yang mana terlihat bahwa perubahan terbesar terjadi pada intercept, sedangkan untuk beta1 tidak terjadi perbedaan yang begitu jauh.