Handling Autocorrelation

by Muhammad Nachnoer Novatron Fitra Arss

A. Input Data

Memanggil Packages

lapply(c("readxl","dplyr","ggplot2","hrbrthemes","car","orcutt","lawstat","knitr"),library,character.only=T)[[1]]

## [1] "readxl"    "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"

Memasukkan data awal (IPM 34 Provinsi di Indonesia)

ipm<-read_excel("C:/Users/falco/Documents/IPMPDW.xlsx")
colnames(ipm)[-1]<-as.numeric(colnames(ipm)[-1])
kable(ipm[1:5,])

Provinsi	2010	2011	2012	2013	2014	2015	2016	2017	2018	2019	2020	2021
ACEH	67.09	67.45	67.81	68.30	68.81	69.45	70.00	70.60	71.19	71.90	71.99	72.18
SUMATERA UTARA	67.09	67.34	67.74	68.36	68.87	69.51	70.00	70.57	71.18	71.74	71.77	72.00
SUMATERA BARAT	67.25	67.81	68.36	68.91	69.36	69.98	70.73	71.24	71.73	72.39	72.38	72.65
RIAU	68.65	68.9	69.15	69.91	70.33	70.84	71.20	71.79	72.44	73.00	72.71	72.94
JAMBI	65.39	66.14	66.94	67.76	68.24	68.89	69.62	69.99	70.65	71.26	71.29	71.63

Mentranspose data dengan menghilangkan baris pertama dan menyesuaikan nama kolom

IPM<-as.data.frame(t(ipm))[-1,];colnames(IPM)<-ipm$Provinsi
kable(IPM[,1:4])

	ACEH	SUMATERA UTARA	SUMATERA BARAT	RIAU
2010	67.09	67.09	67.25	68.65
2011	67.45	67.34	67.81	68.9
2012	67.81	67.74	68.36	69.15
2013	68.30	68.36	68.91	69.91
2014	68.81	68.87	69.36	70.33
2015	69.45	69.51	69.98	70.84
2016	70.00	70.00	70.73	71.20
2017	70.60	70.57	71.24	71.79
2018	71.19	71.18	71.73	72.44
2019	71.90	71.74	72.39	73.00
2020	71.99	71.77	72.38	72.71
2021	72.18	72.00	72.65	72.94

Membuat subset data baru untuk Provinsi DIY dengan peubah Tahun dan IPM yang masing-masing telah di reformat menjadi peubah numerik

diy<-data.frame(Tahun=as.numeric(rownames(IPM)),IPM=as.numeric(IPM$`DI YOGYAKARTA`));diy

##    Tahun   IPM
## 1   2010 75.37
## 2   2011 75.93
## 3   2012 76.15
## 4   2013 76.44
## 5   2014 76.81
## 6   2015 77.59
## 7   2016 78.38
## 8   2017 78.89
## 9   2018 79.53
## 10  2019 79.99
## 11  2020 79.97
## 12  2021 80.22

B.Eksplorasi

Scatterplot

ggplot(diy,aes(x=Tahun,y=IPM))+geom_point(cex=5,col="coral")+geom_line(lwd=2,col="green")+geom_smooth(method = "lm",lwd=2,col="lightblue")+theme_tinyhand(axis_title_just = "center")+theme(plot.background = element_rect(fill="#000044",color="#000044"),plot.title = element_text(color = "white"),axis.title = element_text(color="white"),axis.text = element_text(color="white"))+labs(title="Scatterplot Tahun VS IPM", x="\nTahun", y="IPM\n")

cor(diy)

##           Tahun       IPM
## Tahun 1.0000000 0.9879122
## IPM   0.9879122 1.0000000

Interpretasi: IPM cenderung meningkat tiap tahunnya dan keduanya memiliki korelasi yang kuat. Artinya, terdapat hubungan positif kuat antara urutan waktu dengan IPM di DI Yogyakarta.

Deteksi Autokol

model<-lm(IPM~Tahun,diy);par(mfrow=c(1,2))
plot(1:12,model$residuals, col="blue", main="Sisaan vs Order",xlab="Sisaan", ylab="Order");lines(1:12,model$residuals,col="coral",lwd=2);acf(model$residuals,main="ACF Plot")

Interpretasi: Plot sisaan terhadap order membentuk pola tertentu, artinya sisaan tidak saling bebas atau terdapat korelasi serial. Namun, pada ACF Plot, tidak terdapat kolerasi serial yang signifikan karena nilai ACF tiap lag selain lag 0, tidak ada yang melebihi batas.

C.Pengujian

Durbin Watson Test

dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.77348, p-value = 0.00144
## alternative hypothesis: true autocorrelation is greater than 0

Interpretasi: P-value < 0,05, artinya cukup bukti untuk menyatakan bahwa terdapat autokorelasi positif pada model.

Breusch-Godfrey Test

bgtest(model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 4.2573, df = 1, p-value = 0.03908

Interpretasi: P-value < 0.05, artinya cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada model.

D. Penanganan

Fungsi Hildreth lu

myhildr<-function(x,y,metrics=c("RSE","MSE","RMSE")){
  library(HoRM)
  dat<-data.frame()
  for(rho in seq(0.01,0.99,by=0.01)){
    h<-hildreth.lu(y,x,rho)
    rse<-summary(h)$sigma
    a<-anova(h)
    mse<-tail(a$`Mean Sq`,1);rmse<-sqrt(mse)
    da<-cbind(rho,"RSE"=rse,"MSE"=mse,"RMSE"=rmse)
    dat<-rbind(dat,da)
  }
  d<-dat[which.min(dat[,which(colnames(dat)==metrics)]),c("rho",metrics)]
  h<-hildreth.lu(y,x,d$rho);s<-summary(h)
  return(list(overall.rho=dat,rho.optimum=d,model=h,summary=s))
}

Output Hildreth lu

m<-myhildr(diy$Tahun,diy$IPM,"RSE");m

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

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## $overall.rho
##     rho       RSE        MSE      RMSE
## 1  0.01 0.2997192 0.08983162 0.2997192
## 2  0.02 0.2982008 0.08892372 0.2982008
## 3  0.03 0.2966976 0.08802944 0.2966976
## 4  0.04 0.2952097 0.08714878 0.2952097
## 5  0.05 0.2937375 0.08628174 0.2937375
## 6  0.06 0.2922812 0.08542831 0.2922812
## 7  0.07 0.2908410 0.08458850 0.2908410
## 8  0.08 0.2894172 0.08376231 0.2894172
## 9  0.09 0.2880100 0.08294974 0.2880100
## 10 0.10 0.2866196 0.08215078 0.2866196
## 11 0.11 0.2852463 0.08136545 0.2852463
## 12 0.12 0.2838903 0.08059373 0.2838903
## 13 0.13 0.2825520 0.07983562 0.2825520
## 14 0.14 0.2812315 0.07909114 0.2812315
## 15 0.15 0.2799290 0.07836027 0.2799290
## 16 0.16 0.2786450 0.07764302 0.2786450
## 17 0.17 0.2773795 0.07693939 0.2773795
## 18 0.18 0.2761329 0.07624938 0.2761329
## 19 0.19 0.2749054 0.07557298 0.2749054
## 20 0.20 0.2736973 0.07491020 0.2736973
## 21 0.21 0.2725088 0.07426104 0.2725088
## 22 0.22 0.2713402 0.07362550 0.2713402
## 23 0.23 0.2701917 0.07300357 0.2701917
## 24 0.24 0.2690637 0.07239526 0.2690637
## 25 0.25 0.2679563 0.07180057 0.2679563
## 26 0.26 0.2668698 0.07121950 0.2668698
## 27 0.27 0.2658045 0.07065205 0.2658045
## 28 0.28 0.2647607 0.07009821 0.2647607
## 29 0.29 0.2637385 0.06955799 0.2637385
## 30 0.30 0.2627383 0.06903139 0.2627383
## 31 0.31 0.2617602 0.06851840 0.2617602
## 32 0.32 0.2608046 0.06801904 0.2608046
## 33 0.33 0.2598717 0.06753329 0.2598717
## 34 0.34 0.2589617 0.06706116 0.2589617
## 35 0.35 0.2580749 0.06660265 0.2580749
## 36 0.36 0.2572115 0.06615775 0.2572115
## 37 0.37 0.2563717 0.06572647 0.2563717
## 38 0.38 0.2555559 0.06530881 0.2555559
## 39 0.39 0.2547641 0.06490477 0.2547641
## 40 0.40 0.2539967 0.06451434 0.2539967
## 41 0.41 0.2532539 0.06413754 0.2532539
## 42 0.42 0.2525358 0.06377435 0.2525358
## 43 0.43 0.2518428 0.06342477 0.2518428
## 44 0.44 0.2511749 0.06308882 0.2511749
## 45 0.45 0.2505324 0.06276648 0.2505324
## 46 0.46 0.2499155 0.06245776 0.2499155
## 47 0.47 0.2493244 0.06216266 0.2493244
## 48 0.48 0.2487593 0.06188118 0.2487593
## 49 0.49 0.2482203 0.06161331 0.2482203
## 50 0.50 0.2477076 0.06135907 0.2477076
## 51 0.51 0.2472214 0.06111844 0.2472214
## 52 0.52 0.2467619 0.06089142 0.2467619
## 53 0.53 0.2463291 0.06067803 0.2463291
## 54 0.54 0.2459233 0.06047825 0.2459233
## 55 0.55 0.2455445 0.06029209 0.2455445
## 56 0.56 0.2451929 0.06011955 0.2451929
## 57 0.57 0.2448686 0.05996062 0.2448686
## 58 0.58 0.2445717 0.05981532 0.2445717
## 59 0.59 0.2443023 0.05968363 0.2443023
## 60 0.60 0.2440606 0.05956556 0.2440606
## 61 0.61 0.2438465 0.05946110 0.2438465
## 62 0.62 0.2436601 0.05937027 0.2436601
## 63 0.63 0.2435016 0.05929305 0.2435016
## 64 0.64 0.2433710 0.05922945 0.2433710
## 65 0.65 0.2432683 0.05917946 0.2432683
## 66 0.66 0.2431935 0.05914310 0.2431935
## 67 0.67 0.2431468 0.05912035 0.2431468
## 68 0.68 0.2431280 0.05911122 0.2431280
## 69 0.69 0.2431372 0.05911571 0.2431372
## 70 0.70 0.2431744 0.05913381 0.2431744
## 71 0.71 0.2432397 0.05916554 0.2432397
## 72 0.72 0.2433329 0.05921088 0.2433329
## 73 0.73 0.2434540 0.05926984 0.2434540
## 74 0.74 0.2436030 0.05934241 0.2436030
## 75 0.75 0.2437798 0.05942860 0.2437798
## 76 0.76 0.2439845 0.05952842 0.2439845
## 77 0.77 0.2442168 0.05964185 0.2442168
## 78 0.78 0.2444768 0.05976889 0.2444768
## 79 0.79 0.2447643 0.05990956 0.2447643
## 80 0.80 0.2450792 0.06006384 0.2450792
## 81 0.81 0.2454216 0.06023174 0.2454216
## 82 0.82 0.2457911 0.06041326 0.2457911
## 83 0.83 0.2461877 0.06060839 0.2461877
## 84 0.84 0.2466113 0.06081714 0.2466113
## 85 0.85 0.2470618 0.06103951 0.2470618
## 86 0.86 0.2475389 0.06127550 0.2475389
## 87 0.87 0.2480426 0.06152511 0.2480426
## 88 0.88 0.2485726 0.06178833 0.2485726
## 89 0.89 0.2491288 0.06206517 0.2491288
## 90 0.90 0.2497111 0.06235563 0.2497111
## 91 0.91 0.2503192 0.06265971 0.2503192
## 92 0.92 0.2509530 0.06297740 0.2509530
## 93 0.93 0.2516122 0.06330871 0.2516122
## 94 0.94 0.2522967 0.06365364 0.2522967
## 95 0.95 0.2530063 0.06401219 0.2530063
## 96 0.96 0.2537407 0.06438436 0.2537407
## 97 0.97 0.2544998 0.06477014 0.2544998
## 98 0.98 0.2552833 0.06516954 0.2552833
## 99 0.99 0.2560909 0.06558256 0.2560909
## 
## $rho.optimum
##     rho      RSE
## 68 0.68 0.243128
## 
## $model
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##   -270.8356       0.4586  
## 
## 
## $summary
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32527 -0.20946  0.09704  0.17256  0.30371 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -270.83558   46.78291  -5.789 0.000263 ***
## x              0.45858    0.07244   6.330 0.000136 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2431 on 9 degrees of freedom
## Multiple R-squared:  0.8166, Adjusted R-squared:  0.7962 
## F-statistic: 40.07 on 1 and 9 DF,  p-value: 0.000136

Interpretasi: Didapatkan model dengan penduga, b0 = -270.84 dan b1 = 0.46, dari rho optimum = 0,68 beserta besaran metrics, RSE = 0,24 dan R-Square = 0,82. Dalam hal ini, model juga memiliki peforma yang tinggi, dan merupakan model yang layak karena p-value uji F < 0,05.

Pengujian pasca-Hildreth lu

bgtest(m$model)#tidak ada autokol

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  m$model
## LM test = 2.1982, df = 1, p-value = 0.1382

dwtest(m$model)#ada autokol

## 
##  Durbin-Watson test
## 
## data:  m$model
## DW = 1.06, p-value = 0.01373
## alternative hypothesis: true autocorrelation is greater than 0

runs.test(m$model$residuals)#tidak ada autokol

## 
##  Runs Test - Two sided
## 
## data:  m$model$residuals
## Standardized Runs Statistic = -1.3416, p-value = 0.1797

Interpretasi: Pengujian menggunakan Breusch-Godfrey p-value > 0,05, sehingga menghasilkan keputusan, bahwa terdapat cukup bukti untuk menyatakan bahwa tidak terdapat autokorelasi pada model. Namun, ketika diuji menggunakan Durbin-Watson, p-value < 0,05, sehingga menghasilkan keputusan yang bertolak belakang dari sebelumnya, yaitu terdapat cukup bukti bahwa terdeteksi ada autokorelasi positif pada model. Maka, dilakukan pengujian sekali lagi dengan metode yang berbeda, yaitu Runs-test. P-value Runs-test > 0,05, sehingga menghasilkan keputusan yang sama dengan Breusch Godfrey, yaitu tidak terdapat autokorelasi. Dengan demikian, dapat disimpulkan bahwa autokorelasi sudah teratasi.

Model Hildreth lu Pasca-Transformasi

boh<-m$model$coefficients[1]
b1h<-m$model$coefficients[2]
cat("IPM =",boh/(1-m$rho.optimum$rho),"Tahun","+",b1h)

## IPM = -846.3612 Tahun + 0.4585795

Cochrane orcutt

c<-cochrane.orcutt(model);c

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = IPM ~ Tahun, data = diy)
## 
##  number of interaction: 16
##  rho 0.681705
## 
## Durbin-Watson statistic 
## (original):    0.77348 , p-value: 1.44e-03
## (transformed): 1.06089 , p-value: 1.381e-02
##  
##  coefficients: 
## (Intercept)       Tahun 
## -845.830309    0.458316

Interpretasi:Didapatkan model dengan penduga, b0 = -845.83 dan b1 = 0.46, dari rho optimum = 0,68. Selain itu, p-value hasil uji Durbin Watson setelah transformasi juga > 0,05 yaitu sebesar 0,138. Artinya, terdapat cukup bukti bahwa autokorelasi sudah teratasi.

summary(c)

## Call:
## lm(formula = IPM ~ Tahun, data = diy)
## 
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -845.83031  146.98048  -5.755 0.0002747 ***
## Tahun          0.45832    0.07283   6.293 0.0001422 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2431 on 9 degrees of freedom
## Multiple R-squared:  0.8148 ,  Adjusted R-squared:  0.7942
## F-statistic: 39.6 on 1 and 9 DF,  p-value: < 1.422e-04
## 
## Durbin-Watson statistic 
## (original):    0.77348 , p-value: 1.44e-03
## (transformed): 1.06089 , p-value: 1.381e-02

##R-square 0.82
##RSE 0.24

Interpretasi: Besaran metrics yang didapatkan sama dengan hasil dari hildreth lu, yaitu RSE = 0,24 dan R-Square = 0,82. Dalam hal ini, model juga memiliki peforma yang tinggi, dan merupakan model yang layak karena p-value uji F < 0,05.

Model Orcutt Pasca-Transformasi

boc<-c$coefficients[1]
b1c<-c$coefficients[2]
cat("IPM =",boc/(1-c$rho),"Tahun","+",b1c)

## IPM = -2657.375 Tahun + 0.4583163