Include Package

##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)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(TTR)
library(tseries)
library(lmtest) #uji-Durbin Watson
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(orcutt) #Cochrane-Orcutt
library(HoRM)#Hildreth Lu

Input Data

Data berupa data IPM tahun 2010-2021

dataregresi <- read_xlsx("D:/Kuliah/Sem7/MPDW/Minggu 3/TugasMinggu3/dataIPM.xlsx")
head(dataregresi)
## # A tibble: 6 × 2
##   tahun     y
##   <dbl> <dbl>
## 1  2010  66.1
## 2  2011  66.6
## 3  2012  67.2
## 4  2013  68.0
## 5  2014  68.8
## 6  2015  69.5

Eksplorasi Data

Diketahui data bersifat trend

x <- dataregresi$"tahun"
y <- dataregresi$"y"

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

cor(x,y)
## [1] 0.9925192

Model Regresi Deret Waktu

model <- lm(y~x)
summary(model)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.51436 -0.20550  0.09494  0.20190  0.31494 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.106e+03  4.572e+01  -24.19 3.32e-10 ***
## x            5.832e-01  2.269e-02   25.71 1.82e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2713 on 10 degrees of freedom
## Multiple R-squared:  0.9851, Adjusted R-squared:  0.9836 
## F-statistic: 660.9 on 1 and 10 DF,  p-value: 1.822e-10

Deteksi Autokorelasi

#Residual plot#

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

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(2,1))
acf(resi1)
pacf(resi1)

Uji Durbin-Watson

lmtest::dwtest(model,alternative = 'two.sided')
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.53458, p-value = 0.0001388
## alternative hypothesis: true autocorrelation is not 0

nilai p < 5%. Maka tolak H0, yang artinya cukup bukti bahwa terdapat autokorelasi pada taraf 5%

Penanganan Autokorelasi

  1. Metode Cochrane-Orcutt
#Cochrane-Orcutt
modelco <- orcutt::cochrane.orcutt(model)
## Warning in orcutt::cochrane.orcutt(model): Did not converge
modelco 
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x)
## 
##  number of interaction: 100
##  rho 0.923636
## 
## Durbin-Watson statistic 
## (original):    0.53458 , p-value: 6.938e-05
## (transformed): NA , p-value: NA
##  
##  coefficients: 
## [1] NA
# Hasil sebelum Cochrane-Orcutt
lmtest::dwtest(model,alternative = 'two.sided') 
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.53458, p-value = 0.0001388
## alternative hypothesis: true autocorrelation is not 0
# Hasil setelah Cochrane-Orcutt
lmtest::dwtest(modelco,alternative='two.sided')
## 
##  Durbin-Watson test
## 
## data:  modelco
## DW = 1.4187, p-value = 0.146
## alternative hypothesis: true autocorrelation is not 0

Diketahui setelah digunakan fungsi Cochrane-Orcutt, nilai statistik uji DW bertambah dari 0.534 menjadi 1.418 dan nilai p meningkat dari 0.0001388 ke 0.146

  1. Metode 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,1, by= 0.1))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0,40 memiliki SSE Terkecil
##    rho    SSE
## 1  0.1 0.6116
## 2  0.2 0.5399
## 3  0.3 0.4754
## 4  0.4 0.4181
## 5  0.5 0.3680
## 6  0.6 0.3252
## 7  0.7 0.2895
## 8  0.8 0.2611
## 9  0.9 0.2399
## 10 1.0 0.3832
#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.9
r <- seq(0.8, 1.0, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0,99 memiliki SSE Terkecil
##     rho    SSE
## 1  0.80 0.2611
## 2  0.81 0.2587
## 3  0.82 0.2563
## 4  0.83 0.2540
## 5  0.84 0.2518
## 6  0.85 0.2496
## 7  0.86 0.2475
## 8  0.87 0.2455
## 9  0.88 0.2436
## 10 0.89 0.2417
## 11 0.90 0.2399
## 12 0.91 0.2382
## 13 0.92 0.2365
## 14 0.93 0.2349
## 15 0.94 0.2334
## 16 0.95 0.2320
## 17 0.96 0.2306
## 18 0.97 0.2293
## 19 0.98 0.2281
## 20 0.99 0.2270
## 21 1.00 0.3832
#grafik SSE optimum
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

# Model terbaik
modelhl <- hildreth.lu.func(0.99, model)
summary(modelhl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26082 -0.10654  0.02656  0.12776  0.17137 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   68.302     32.024   2.133   0.0617 .
## x             -3.171      1.514  -2.094   0.0658 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1588 on 9 degrees of freedom
## Multiple R-squared:  0.3276, Adjusted R-squared:  0.2529 
## F-statistic: 4.385 on 1 and 9 DF,  p-value: 0.06576
# Hasil sebelum Hildreth-lu
lmtest::dwtest(model,alternative = 'two.sided') 
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.53458, p-value = 0.0001388
## alternative hypothesis: true autocorrelation is not 0
# Hasil setelah Hidlreth-lu
lmtest::dwtest(modelhl,alternative='two.sided')
## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 1.5498, p-value = 0.2287
## alternative hypothesis: true autocorrelation is not 0

Diketahui setelah digunakan fungsi Hildreth-lu, nilai statistik uji DW bertambah lebih dari 1 dan nilai p meningkat dari 0.0001388 ke 0.2287

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

Kesimpulan

Setelah digunakan fungsi Cochrane-Orcutt dan Hildreth-Lu, nilai p yang sebelumnya kurang dari 5% berubah menjadi lebih dari 5%. Artinya jika sebelumnya cukup bukti pada taraf 5% bahwa terdapat autokorelasi, setelah digunakan fungsi-fungsi tersebut, kesimpulannya berubah menjadi tidak cukup bukti bahwa terjadi autokorelasi.