Memanggil Library yang digunakan

Diperlukan beberapa package dalam eksekusi regresi time series dan juga mengecek gejala autokorelasi dan penanganannya. Berikut adalah list packagenya :

• tidyverse
• TTR
• forecast
• readxl
• tseries
• lmtest
• orcutt
• HoRM

lapply(c("tidyverse","forecast","TTR","readxl",
         "tseries","lmtest","orcutt","HoRM","rmarkdown"),
       library,character.only=T)

Deklarasi Data Indeks Pembangunan Bali 2010-2021

IPM <- read_excel("D:\\Kuliah Semester 5\\MPDW\\Tugas Individu MPDW.xlsx",
                  sheet = 3)
knitr::kable(IPM,align = "c")

Indeks Pembangunan Manusia

Indeks Pembangunan Manusia merupakan indikator penting untuk mengukur keberhasilan dalam upaya membangun kualitas hidup masyarakat. Kecepatan IPM menggambarkan upaya yang dilakukan untuk meningkatkan pembangunan manusia dalam suatu periode.

Time Series Plot

p <- ggplot(IPM, aes(x=Tahun, y=IPM)) +
  geom_line(lwd=1.2,col="red3")
p +labs(x="Year",y = "Indeks Pembangunan",
         title="Time Series Plot Indeks Pembangunan Manusia Bali ",
        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 Bali mengalami kenaikan secara konstan dari tahun ke tahun. Pada tahun 2020 ,IPM mengalami peningkatan sebesar 0,12 persen. Pada tahun 2019, IPM Bali mencapai 75,38 dan selanjutnya pada tahun 2020 naik mencapai 75,5 .Walaupun masih tergolong tinggi, peningkatan IPM di Bali mengalami perlambatan jika dibandingkan dengan tahun-tahun sebelumnya. Pada tahun sebelumnya rata-rata peningkatan IPM mencapai 0,86 persen.

Deteksi Autokorelasi

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

#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",main = "Histogram Residual")

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)

#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(1,2))
acf(resi1)
pacf(resi1)

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

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

bgtest(IPM ~ Tahun, data=IPM, order=1)

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

Cochran Orcutt

modelco <- orcutt::cochrane.orcutt(model,convergence=6,max.iter = 1000)
modelco 

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = IPM ~ Tahun, data = IPM)
## 
##  number of interaction: 177
##  rho 0.735319
## 
## Durbin-Watson statistic 
## (original):    0.87546 , p-value: 3.527e-03
## (transformed): 2.13015 , p-value: 4.363e-01
##  
##  coefficients: 
## (Intercept)       Tahun 
## -682.302244    0.375136

rho <- modelco$rho
rho

## [1] 0.7353396

y[-1]

##  [1] 70.87 71.62 72.09 72.48 73.27 73.65 74.30 74.77 75.38 75.50 75.69

y[-12]

##  [1] 70.10 70.87 71.62 72.09 72.48 73.27 73.65 74.30 74.77 75.38 75.50

(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 dengan peubah transformasi

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

Pengecekan Autokorelasi

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

bgtest(modelcorho)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelcorho
## LM test = 0.20182, df = 1, p-value = 0.6533

Transformasi Balik

cat("IPM = ", coef(modelcorho)[1]/(1-rho), "+", coef(modelcorho)[2]," Tahun", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal

## IPM = -682.2747+0.3751221 Tahun

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.999, by= 0.001))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
tab$rho[which.min(tab$SSE)]#rho optimal

## [1] 0.735

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

Model dengan peubah transformasi

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

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18521 -0.15406 -0.03417  0.13618  0.24766 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -180.9202    33.8105  -5.351 0.000462 ***
## x              0.3753     0.0632   5.939 0.000218 ***
## ---
## 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.7967, Adjusted R-squared:  0.7741 
## F-statistic: 35.27 on 1 and 9 DF,  p-value: 0.0002183

Pengecekan Asumsi Autokorelasi

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

lmtest::dwtest(modelhl,alternative = "two.sided")

## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 2.1296, p-value = 0.8718
## alternative hypothesis: true autocorrelation is not 0

bgtest(modelcorho)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelcorho
## LM test = 0.20182, df = 1, p-value = 0.6533

bgtest(modelhl)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelhl
## LM test = 0.20115, df = 1, p-value = 0.6538

Berdasarkan Hasil di atas dapat disimpulkan hasil transformasi memenuhi asumsi bebas autokorelasi.

Transformasi Balik

cat("IPM = ", coef(modelhl)[1]/(1-0.735), "+", coef(modelhl)[2]," Tahun", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal

## IPM = -682.7179+0.3753413 Tahun