Package

library(readxl)
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(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
## Warning: package 'zoo' was built under R version 4.1.3
## 
## 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

Data yang digunakan merupakan Data IPM Provinsi Aceh tahun 2010-2021

data <- read_xlsx("C:/Users/Alam/Downloads/Tugas Individu MPDW (1).xlsx", sheet = 3)
head(data)
## # A tibble: 6 x 2
##   tahun  ACEH
##   <dbl> <dbl>
## 1     1  67.1
## 2     2  67.4
## 3     3  67.8
## 4     4  68.3
## 5     5  68.8
## 6     6  69.4

Eksplorasi Data

x <- data$tahun
y <- data$ACEH
tahungraph <- c(2010:2021)

#diagram pencar identifikasi model
plot(tahungraph,y,pch = 20, col = "blue", main = "Scatter Plot X vs Y",
     ylab = "IPM", xlab = "Tahun (Tahun 2010-2021)")

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

Berdasarkan eksplorasi data tahun 2010-2021, didapatkan korelasi positif antara kedua peubah dengan nilai korelasi sebesar 0.9944.

Model Regresi Deret Waktu

#model regresi
model <- lm(y~x)
summary(model)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35564 -0.14094 -0.00592  0.11916  0.38429 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 66.41606    0.12542  529.54  < 2e-16 ***
## x            0.50997    0.01704   29.93 4.06e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2038 on 10 degrees of freedom
## Multiple R-squared:  0.989,  Adjusted R-squared:  0.9879 
## F-statistic: 895.5 on 1 and 10 DF,  p-value: 4.061e-11

Model yang regresi yang didapatkan yakni b0 = 66.416 dan b1 = 0.510 dengan nilai adjusted R-squared 0.9879, yakni sebesar 98.79% nilai y dapat dijelaskan oleh peubah x.

Deteksi Autokorelasi

Grafik Autocol

1. Residual Plot

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

#Diagnostik dengan eksploratif
par(mfrow = c(2,2))

qqnorm(residual1)
qqline(residual1, col = "steelblue", lwd = 2)

plot(fit, residual1, col = "steelblue", pch = 20, xlab = "Sisaan", 
     ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)

hist(residual1, col = "steelblue")

plot(seq(1,12,1), residual1, col = "steelblue", pch = 20, 
     xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,12,1), residual1, col = "red")
abline(a = 0, b = 0, lwd = 2)

2.Plot ACF dan PACF

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

Uji Statistik

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.95635, p-value = 0.01294
## alternative hypothesis: true autocorrelation is not 0

Berdasarkan uji Durbin-Watson, didapatkan bahwa karena p-value > 0.05 H0 ditolak, sehingga terdapat cukup bukti untuk menyatakan bahwa tidak ada autokorelasi pada model tersebut.

2. Breusch-Godfrey Test

H0: tidak ada autokorelasi

H1: ada autokorelasi

lmtest::bgtest(y ~ x, order=1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  y ~ x
## LM test = 2.3518, df = 1, p-value = 0.1251

Berdasarkan uji Breusch-Godfrey, didapatkan bahwa karena p-value > 0.05 H0 diterima, sehingga terdapat cukup bukti untuk menyatakan bahwa tidak ada autokorelasi pada model tersebut.

3. Runs 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.8166, p-value = 0.06928

Berdasarkan uji Runs, didapatkan bahwa karena p-value > 0.05 H0 diterima, sehingga terdapat cukup bukti untuk menyatakan bahwa tidak ada autokorelasi pada model tersebut.

Penanganan Autokorelasi

1. Cochrane-Orcutt

modelco <- orcutt::cochrane.orcutt(model)
modelco 
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x)
## 
##  number of interaction: 18
##  rho 0.535106
## 
## Durbin-Watson statistic 
## (original):    0.95635 , p-value: 6.471e-03
## (transformed): 1.10751 , p-value: 1.797e-02
##  
##  coefficients: 
## (Intercept)           x 
##   66.371964    0.506876
rho <- modelco$rho
y
##  [1] 67.09 67.45 67.81 68.30 68.81 69.45 70.00 70.60 71.19 71.90 71.99 72.18
#transformasi terhadap y dan x
(y.trans <- y[-1]-y[-12]*rho)
##  [1] 31.54971 31.71707 32.01443 32.26223 32.62933 32.83686 33.14255 33.41149
##  [9] 33.80578 33.51585 33.65769
(x.trans <- x[-1]-x[-12]*rho)
##  [1] 1.464894 1.929787 2.394681 2.859574 3.324468 3.789362 4.254255 4.719149
##  [9] 5.184042 5.648936 6.113830
modelcorho <- lm(y.trans~x.trans)
lmtest::dwtest(modelcorho,alternative = 'two.sided') 
## 
##  Durbin-Watson test
## 
## data:  modelcorho
## DW = 1.1075, p-value = 0.03595
## alternative hypothesis: true autocorrelation is not 0
lmtest::bgtest(modelcorho, order=1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelcorho
## LM test = 1.6603, df = 1, p-value = 0.1976

Setelah dilakukan penanganan autokorelasi Cochrane Orcutt, didapatkan bahwa H0 ditolak pada Durbin-Watson test, sementara H0 diterima pada Breusch-Godfrey test. Sehingga, didapatkan bahwa setelah dilaksanakan penanganan Cochrane-Orcutt, terdapat autokorelasi pada Durbin-Watson test, sementara tidak ada autokorelasi pada Breusch-Godfrey test.

#model baru
summary(modelcorho) 
## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29716 -0.08613 -0.04312  0.10930  0.32221 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 30.85590    0.15461  199.57  < 2e-16 ***
## x.trans      0.50688    0.03804   13.32 3.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1855 on 9 degrees of freedom
## Multiple R-squared:  0.9518, Adjusted R-squared:  0.9464 
## F-statistic: 177.6 on 1 and 9 DF,  p-value: 3.14e-07

Setelah dilakukan penanganan Cochran-Orcutt, didapatkan model regresi yakni b0 = 30.856 dan b1 = 0.507 dengan nilai adjusted R-squared 0.9879, yakni sebesar 94.64% nilai y dapat dijelaskan oleh peubah x.

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.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))}))
round(tab, 4) #0,40 memiliki SSE Terkecil
##     rho    SSE
## 1  0.10 0.3543
## 2  0.20 0.3361
## 3  0.30 0.3226
## 4  0.40 0.3139
## 5  0.50 0.3099
## 6  0.60 0.3106
## 7  0.70 0.3160
## 8  0.80 0.3262
## 9  0.90 0.3410
## 10 0.91 0.3428
## 11 0.92 0.3446
## 12 0.93 0.3464
## 13 0.94 0.3483
## 14 0.95 0.3502
## 15 0.96 0.3522
## 16 0.97 0.3542
## 17 0.98 0.3563
## 18 0.99 0.3584
#rho optimal di sekitar 0.5
r <- seq(0.4, 0.6, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0,53 memiliki SSE Terkecil
##     rho    SSE
## 1  0.40 0.3139
## 2  0.41 0.3133
## 3  0.42 0.3127
## 4  0.43 0.3122
## 5  0.44 0.3117
## 6  0.45 0.3113
## 7  0.46 0.3109
## 8  0.47 0.3106
## 9  0.48 0.3103
## 10 0.49 0.3101
## 11 0.50 0.3099
## 12 0.51 0.3098
## 13 0.52 0.3097
## 14 0.53 0.3096
## 15 0.54 0.3096
## 16 0.55 0.3097
## 17 0.56 0.3098
## 18 0.57 0.3099
## 19 0.58 0.3101
## 20 0.59 0.3103
## 21 0.60 0.3106
#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.53, model)
summary(modelhl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29800 -0.08626 -0.04388  0.10878  0.32269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 31.19437    0.15423  202.26  < 2e-16 ***
## x            0.50712    0.03763   13.48 2.85e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1855 on 9 degrees of freedom
## Multiple R-squared:  0.9528, Adjusted R-squared:  0.9475 
## F-statistic: 181.6 on 1 and 9 DF,  p-value: 2.847e-07
# Deteksi autokorelasi
lmtest::dwtest(modelhl)
## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 1.106, p-value = 0.01783
## alternative hypothesis: true autocorrelation is greater than 0
lmtest::bgtest(modelhl, order=1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelhl
## LM test = 1.668, df = 1, p-value = 0.1965

Setelah dilakukan penanganan autokorelasi Hildreth-Lu, didapatkan bahwa H0 ditolak pada Durbin-Watson test, sementara H0 diterima pada Breusch-Godfrey test. Sehingga, didapatkan bahwa setelah dilaksanakan penanganan Hildreth-Lu, terdapat autokorelasi pada Durbin-Watson test, sementara tidak ada autokorelasi pada Breusch-Godfrey test.

# Transformasi Balik
cat("y = ", coef(modelhl)[1]/(1-0.44), "+", coef(modelhl)[2],"x", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal
## y = 55.70424+0.5071199x
#fungsi Hildreth Lu dengan library HoRM
modelhl2 <- HoRM::hildreth.lu(y, x, 0.53) #fungsi hildreth lu dengan rho tertentu
summary(modelhl2)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29800 -0.08626 -0.04388  0.10878  0.32269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 31.19437    0.15423  202.26  < 2e-16 ***
## x            0.50712    0.03763   13.48 2.85e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1855 on 9 degrees of freedom
## Multiple R-squared:  0.9528, Adjusted R-squared:  0.9475 
## F-statistic: 181.6 on 1 and 9 DF,  p-value: 2.847e-07

Setelah dilakukan penanganan Hildreth-Lu, didapatkan model regresi yakni b0 = 31.194 dan b1 = 0.507 dengan nilai adjusted R-squared 0.9475, yakni sebesar 98.79% nilai y dapat dijelaskan oleh peubah x.