Deteksi dan Penanganan Autokorelasi pada Data IPM Provinsi Jawa Barat 2010-2021
M. Ficky Haris Ardiansyah
9/7/2022
#Library yang digunakan
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.1.2## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.3 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.1.2## Warning: package 'tidyr' was built under R version 4.1.1## Warning: package 'readr' was built under R version 4.1.1## Warning: package 'purrr' was built under R version 4.1.1## Warning: package 'dplyr' was built under R version 4.1.1## Warning: package 'stringr' was built under R version 4.1.1## Warning: package 'forcats' was built under R version 4.1.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(ggplot2)
library(dplyr)
library(readxl)
library(forecast)## Warning: package 'forecast' was built under R version 4.1.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(readxl)
library(TTR)## Warning: package 'TTR' was built under R version 4.1.3library(tseries)## Warning: package 'tseries' was built under R version 4.1.3library(lmtest) #uji-Durbin Watson## Warning: package 'lmtest' was built under R version 4.1.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.1.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(orcutt) #Cochrane-Orcutt## Warning: package 'orcutt' was built under R version 4.1.3library(HoRM)#Hildreth Lu## Warning: package 'HoRM' was built under R version 4.1.3library(knitr)## Warning: package 'knitr' was built under R version 4.1.3Input Data IPM Provinsi Jawa Barat 2010-2021
dataregresi <- read_excel("Tugas mpdw jabar.xlsx")
head(dataregresi)## # A tibble: 6 x 2
## Tahun IPM
## <dbl> <dbl>
## 1 2010 66.2
## 2 2011 66.7
## 3 2012 67.3
## 4 2013 68.2
## 5 2014 68.8
## 6 2015 69.5View(dataregresi)Eksplorasi Data
p <- ggplot(dataregresi, aes(x=Tahun, y=IPM)) +
geom_line(lwd=1.3,col="blue")
p +labs(x="Tahun",y = "IPM",
title="Time Series Plot Indeks Pembangunan Manusia Jawa Barat",
subtitle = "Periode 2010 - 2021")+
theme_bw()+
theme(
plot.title = element_text(size = 14L,
face = "bold",
hjust = 0.5),
plot.subtitle = element_text(size = 11L,
face = "plain",
hjust = 0.5)
)+geom_point(size=2) +geom_text(aes(label=paste(IPM,"%")),vjust=-0.8,size=3)
Indeks Pembangunan Manusia di Provinsi Jawa Barat tahun 2010-2021 mengalami kenaikan yang konstan. Pada tahun 2019, IPM Provinsi Jawa Barat mengalami kenaikan paling kecil yaitu senilai 0.06%.
Pemodelan Regresi
x <- dataregresi$Tahun
y <- dataregresi$IPM
model<- lm(y~x, data = dataregresi)
summary(model)##
## Call:
## lm(formula = y ~ x, data = dataregresi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4760 -0.1888 0.1183 0.1785 0.3104
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.146e+03 4.237e+01 -27.05 1.10e-10 ***
## x 6.032e-01 2.102e-02 28.69 6.16e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2514 on 10 degrees of freedom
## Multiple R-squared: 0.988, Adjusted R-squared: 0.9868
## F-statistic: 823.2 on 1 and 10 DF, p-value: 6.159e-11cor(x, y)## [1] 0.9939809Berdasarkan output di atas diperoleh model regresi linier deret waktu yaitu: IPM^ = -1146 + 0.6032 (Tahun)
Deteksi Auotokorelasi
1. Residual Plot
#sisaan dan fitted value
resi1 <- residuals(model)
fit <- predict(model)
fit## 1 2 3 4 5 6 7 8
## 66.29064 66.89386 67.49707 68.10029 68.70351 69.30672 69.90994 70.51316
## 9 10 11 12
## 71.11638 71.71959 72.32281 72.92603#Diagnostik dengan eksploratif
par(mfrow = c(2,2))
qqnorm(resi1)
qqline(resi1, col = "steelblue", lwd = 2)
plot(fit, resi1, col = "steelblue", pch = 20, xlab = "Sisaan",
ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(resi1, col = "steelblue")
plot(seq(1,12,1), resi1, col = "steelblue", pch = 20,
xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,12,1), resi1, col = "red")
abline(a = 0, b = 0, lwd = 2)
1. Normal QQ-Plot Sisaan berada di sepanjang garis normal sehingga mengindikasikan bahwa sisaan menyebar normal
Residual vs Fitted Values Lebar pita relatif sama sehingga mengindikasikan ragam sisaan homogen
Histogram Sisaan menjulur ke kiri yang berarti banyak sisaan yang bernilai besar
Residual vs Urutan Tebaran amatan memiliki pola yang tidak acak. Diduga terdapat autokorelasi
2. ACF dan PACF Plot
#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(1,2))
acf(resi1)
pacf(resi1)
Berdasarkan plot nilai ACF dan PACF, terlihat terdapat garis vertikal yang melewati garis horizontal berwarna biru sehingga diduga terdapat autokorelasi.
Uji Statistik
1. Durbin Watson Test
Asumsi : H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::dwtest(model, alternative = 'two.sided') #ada autokorelasi##
## Durbin-Watson test
##
## data: model
## DW = 0.79521, p-value = 0.003536
## alternative hypothesis: true autocorrelation is not 0Berdasarkan uji Durbin Watson Test, didapatkan nilai p-value (<0.05) sehingga tolak H0. Cukup bukti untuk menyatakan bahwa ada autokorelasi antar amatan pada taraf 5%.
2. Breusch-Godfrey Test
H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::bgtest(y ~ x, data=dataregresi, order=1) #Perform Breusch-Godfrey test for first-order serial correlation##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: y ~ x
## LM test = 3.9403, df = 1, p-value = 0.04714Berdasarkan uji Breusch-Godfrey Test, didapatkan nilai p-value (<0.05) sehingga tolak H0 atau cukup bukti untuk menyatakan bahwa ada autokorelasi antar amatan pada taraf 5%.
3. 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.2259Berdasarkan uji Run’Test didapatkan nilai p-value > 5% sehingga tak tolak H0 sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%
Berdasarkan uji statistik pada uji Durbin Watson dan uji Breusch-Godfrey didapatkan nilai p-value < 0.05 yang artinya Tolak H0 (Terdapat autokorelasi). Oleh karena itu, diperlukan adanya penanganan menggunakan Cochran Orcutt dan Hildreth-lu.
Penanganan Autokorelasi
1. Cochrane-Orcutt
#Interactive method using to solve first order autocorrelation problems.
modelco <- orcutt::cochrane.orcutt(model, convergence = 5, max.iter=1000)
modelco ## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = y ~ x, data = dataregresi)
##
## number of interaction: 412
## rho 0.878545
##
## Durbin-Watson statistic
## (original): 0.79521 , p-value: 1.768e-03
## (transformed): 2.07953 , p-value: 4.014e-01
##
## coefficients:
## (Intercept) x
## -638.590309 0.352236#rho optimum
rho <- modelco$rho
y## [1] 66.15 66.67 67.32 68.25 68.80 69.50 70.05 70.69 71.30 72.03 72.09 72.45y[-1]## [1] 66.67 67.32 68.25 68.80 69.50 70.05 70.69 71.30 72.03 72.09 72.45y[-12]## [1] 66.15 66.67 67.32 68.25 68.80 69.50 70.05 70.69 71.30 72.03 72.09View(rho)#transformasi terhadap y dan x
(y.trans <- y[-1]-y[-12]*rho)## [1] 8.554229 8.747386 9.106331 8.839284 9.056084 8.991103 9.147903 9.195634
## [9] 9.389721 8.808383 9.115670(x.trans <- x[-1]-x[-12]*rho)## [1] 245.1240 245.2454 245.3669 245.4883 245.6098 245.7312 245.8527 245.9742
## [9] 246.0956 246.2171 246.3385#### Model Baru
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho)##
## Call:
## lm(formula = y.trans ~ x.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35835 -0.08547 -0.00451 0.11199 0.26577
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -77.5598 38.7676 -2.001 0.0765 .
## x.trans 0.3522 0.1578 2.233 0.0525 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.201 on 9 degrees of freedom
## Multiple R-squared: 0.3564, Adjusted R-squared: 0.2849
## F-statistic: 4.985 on 1 and 9 DF, p-value: 0.05247Pengecekan Autokorelasi
lmtest::dwtest(modelcorho,alternative = 'two.sided') ##
## Durbin-Watson test
##
## data: modelcorho
## DW = 2.0795, p-value = 0.8028
## alternative hypothesis: true autocorrelation is not 0Berdasarkan uji formal Durbin Watson dan Run`s test p-value > 0.05 Tak Tolak H0, artinya tidak cukup bukti untuk menyatakan bahwa ada autokorelasi pada taraf 5%.
Transformasi Balik
cat("IPM = ", coef(modelcorho)[1]/(1-rho), "+", coef(modelcorho)[2]," Tahun", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal## IPM = -638.5903+0.3522361 TahunMaka model regresi setelah dilakukan transformasi balik adalah IPM^ = -638.5903 + 0.3522361 (Tahun)
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))
}
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))}))
optrho <- which.min(round(tab,4)[,2])
round(tab, 4) #0,90 memiliki SSE Terkecil## rho SSE
## 1 0.10 0.5525
## 2 0.20 0.5071
## 3 0.30 0.4680
## 4 0.40 0.4350
## 5 0.50 0.4083
## 6 0.60 0.3878
## 7 0.70 0.3735
## 8 0.80 0.3655
## 9 0.90 0.3636
## 10 0.91 0.3638
## 11 0.92 0.3640
## 12 0.93 0.3643
## 13 0.94 0.3646
## 14 0.95 0.3650
## 15 0.96 0.3655
## 16 0.97 0.3660
## 17 0.98 0.3666
## 18 0.99 0.3672Dari hasil diatas, diketahui bahwa nilai ρ optimum yang meminimumkan SSE adalah 0.9
#grafik rho dan SSE
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)
#rho optimal di sekitar 0.95
r <- seq(0.8, 0.9, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0,90 memiliki SSE Terkecil## rho SSE
## 1 0.80 0.3655
## 2 0.81 0.3650
## 3 0.82 0.3646
## 4 0.83 0.3643
## 5 0.84 0.3640
## 6 0.85 0.3638
## 7 0.86 0.3636
## 8 0.87 0.3635
## 9 0.88 0.3635
## 10 0.89 0.3635
## 11 0.90 0.3636#grafik SSE optimum
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)
Model Baru
Terlihat bahwa nilai ρ optimum yang menghasilkan SES terkecil adalah sebesar 0.88. Dengan demikian, model optimum yang dibuat dengan metode Hidreth-Lu yaitu:
modelhl <- hildreth.lu.func(0.88, model)
summary(modelhl)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35859 -0.08506 -0.00473 0.11187 0.26568
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -75.8130 38.7695 -1.955 0.0822 .
## x 0.3489 0.1597 2.185 0.0567 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.201 on 9 degrees of freedom
## Multiple R-squared: 0.3466, Adjusted R-squared: 0.274
## F-statistic: 4.774 on 1 and 9 DF, p-value: 0.05672Deteksi Autokorelasi
lmtest::dwtest(modelhl,alternative = 'two.sided')##
## Durbin-Watson test
##
## data: modelhl
## DW = 2.0823, p-value = 0.8066
## alternative hypothesis: true autocorrelation is not 0Pengecekan Asumsi Autokorelasi
#fungsi Hildreth Lu dengan library HoRM
modelhl2 <- HoRM::hildreth.lu(y, x, 0.88) #fungsi hildreth lu dengan rho tertentu
summary(modelhl2)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35859 -0.08506 -0.00473 0.11187 0.26568
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -75.8130 38.7695 -1.955 0.0822 .
## x 0.3489 0.1597 2.185 0.0567 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.201 on 9 degrees of freedom
## Multiple R-squared: 0.3466, Adjusted R-squared: 0.274
## F-statistic: 4.774 on 1 and 9 DF, p-value: 0.05672Model Hildreth-Lu dengan memanfaatkan nilai rho yang didapatkan melalui metode Cochrane-Orcutt juga menghasilkan model yang sudah bebas dari autokorelasi
Transformasi Balik
cat("y = ", coef(modelhl)[1]/(1-0.44), "+", coef(modelhl)[2],"x", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal## y = -135.3804+0.3488788xModel regresi yang terbentuk setelah dilakukan penanganan autokorelasi yaitu IPM^ = -135.3804+0.3488788 (Tahun)