Package dan Data yang Digunakan

Package

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

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 zoo

library(TTR)

## Warning: package 'TTR' was built under R version 4.1.3

library(ggplot2)
library(tseries)

## Warning: package 'tseries' was built under R version 4.1.3

library(lmtest) #uji-Durbin Watson

## Warning: package 'lmtest' was built under R version 4.1.3

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(orcutt) #Cochrane-Orcutt

## Warning: package 'orcutt' was built under R version 4.1.3

library(HoRM)#Hildreth Lu

## Warning: package 'HoRM' was built under R version 4.1.3

Data yang Digunakan

ipmjabar <- read_excel("D://Kuliah//Semester 5//MPDW//Data IPM Jabar.xlsx")

Eksplorasi Data

Time Series Plot

z <- ggplot(ipmjabar, aes(x=Tahun, y=IPM)) +
  geom_line(lwd=1.2,col="red3")
z +labs(x="Year",y = "Indeks Pembangunan",
         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)

Korelasi x dan y

x <- ipmjabar$Tahun
y <- ipmjabar$IPM
cor(x,y)

## [1] 0.9939809

Data yang digunakan sebanyak 12 pengamatan. Secara visual, terlihat bahwa IPM dan tahun memiliki hubungan positif dengan korelasi sebesar 0.99. Artinya, IPM dan tahun memiliki hubungan yang sangat kuat.

Model Regresi Deret Waktu

x <- ipmjabar$Tahun
y <- ipmjabar$IPM
model <- lm(IPM~Tahun, data = ipmjabar)
summary(model)

## 
## Call:
## lm(formula = IPM ~ Tahun, data = ipmjabar)
## 
## 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 ***
## Tahun        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-11

Berdasarkan output diatas, diperoleh model regresi linier data deret waktu yaitu: IPM^ = -1.146 + 0.6032 (Tahun)

Deteksi Autokorelasi

Pendeteksian autokorelasi pada model regresi dapat dilakukan dengan dua cara, yaitu secara eksploratif dan secara formal. Secara eksploratif dapat menggunakan plot sisaan dan plot ACF/PACF. Secara formal dapat menggunakan uji Durbin-Watson dan Run’s Test.

1. Plot Sisaan

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

Berdasarkan plot sisaan vs order, terlihat amatan memiliki pola yang tidak acak. Dapat dilihat karena empat sisaan pertama membentuk garis naik, selanjutnya sisaan keempat sampai kesepuluh cenderung konstan, kemudian tiga sisaan terakhir membentuk garis turun. Oleh karena itu, diduga terdapat autokorelasi berdasarkan plot sisaan dan urutan.

2. ACF dan PACF Plot

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

Berdasarkan plot ACF dan PACF, terlihat ada garis vertikal yang melewati garis horizontal yang berwarna biru sehingga dapat diduga ada autokorelasi.

Uji Formal

1. Durbin Watson Test

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 0

Berdasarkan uji Durbin Watson didapatkan p-value < 0.05 yang artinya Tolak H0 (cukup bukti untuk menyatakan terdapat autokorelasi pada taraf 5%)

2. Run’s Test

H0: tidak ada autokorelasi H1: ada autokorelasi

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

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

Berdasarkan Run’s test didapatkan p-value > 0.05 yang artinya tak tolak H0 (cukup bukti untuk menyatakan terdapat autokorelasi pada taraf 5%)

Penanganan Autokorelasi

1. Cochrane-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 = ipmjabar)
## 
##  number of interaction: 748
##  rho 0.879734
## 
## Durbin-Watson statistic 
## (original):    0.79521 , p-value: 1.768e-03
## (transformed): 2.08182 , p-value: 4.03e-01
##  
##  coefficients: 
## (Intercept)       Tahun 
##   -633.0361      0.3495

#rho optimum
rho <- modelco$rho
rho

## [1] 0.8797335

y[-1]

##  [1] 66.67 67.32 68.25 68.80 69.50 70.05 70.69 71.30 72.03 72.09 72.45

y[-12]

##  [1] 66.15 66.67 67.32 68.25 68.80 69.50 70.05 70.69 71.30 72.03 72.09

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

##  [1] 8.475627 8.668166 9.026339 8.758187 8.974333 8.908520 9.064666 9.111637
##  [9] 9.304999 8.722794 9.030010

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

##  [1] 242.7356 242.8559 242.9761 243.0964 243.2167 243.3369 243.4572 243.5775
##  [9] 243.6977 243.8180 243.9383

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

## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35855 -0.08514 -0.00469  0.11189  0.26569 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -76.1330    38.7692  -1.964   0.0812 .
## x.trans       0.3495     0.1593   2.194   0.0559 .
## ---
## 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.3484, Adjusted R-squared:  0.276 
## F-statistic: 4.812 on 1 and 9 DF,  p-value: 0.05592

Pengecekan Autokorelasi

Uji Durbin Watson

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

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

Run’s Test

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

## 
##  Runs Test - Two sided
## 
## data:  resid(modelcorho)
## Standardized Runs Statistic = -0.93314, p-value = 0.3507

Berdasarkan uji formal Durbin Watson dan Run`s test, diperoleh p-value > 0.05 sehingga tak tolak H0, artinya tidak cukup bukti untuk menyatakan bahwa sisaan tidak ada autokorelasi pada taraf 5%.

Transformasi Balik

cat("IPM = ", coef(modelcorho)[1]/(1-rho), "+", coef(modelcorho)[2]," Tahun", sep = "")

## IPM = -633.0361+0.3494998 Tahun

2. Hidreth-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)[optrho,]

##   rho    SSE
## 9 0.9 0.3636

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

r <- seq(0.2,1.0, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r,function(i){deviance(hildreth.lu.func(i, model))}))
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

modelhl <- hildreth.lu.func(0.9, model)
summary(modelhl)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36185 -0.07965 -0.00773  0.11035  0.26443 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -51.7880    38.8082  -1.334    0.215
## x             0.2928     0.1916   1.528    0.161
## 
## Residual standard error: 0.201 on 9 degrees of freedom
## Multiple R-squared:  0.206,  Adjusted R-squared:  0.1177 
## F-statistic: 2.335 on 1 and 9 DF,  p-value: 0.1609

Pengecekan Autokorelasi

Uji Durbin Watson

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

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

Run’s Test

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

## 
##  Runs Test - Two sided
## 
## data:  resid(modelhl)
## Standardized Runs Statistic = -0.93314, p-value = 0.3507

Berdasarkan uji formal Durbin Watson dan Run`s test, diperoleh p-value > 0.05 sehingga tak tolak H0, artinya tidak cukup bukti untuk menyatakan bahwa sisaan tidak ada autokorelasi pada taraf 5%.

Transformasi Balik

cat("y = ", coef(modelhl)[1]/(1-0.9), "+", coef(modelhl)[2],"Tahun", sep = "")

## y = -517.8795+0.2928182Tahun