Tugas Individu 1 - Metode Peramalan Deret Waktu

Pra-Processing Data

Untuk menangani kasus penanganan autokorelasi, diperlukan beberapa library untuk menunjang proses analisis data sebagai berikut.

library(readxl)
library(tidyverse)
library(TTR)
library(forecast)
library(tseries)
library(lmtest)
library(orcutt)
library(HoRM)
library(lawstat)

Pemanggilan data yang berformat .xlsx, import data menuju R dapat dilakukan menggunakan fungsi read_excel dengan rangkaian sintaks sebagai berikut.

setwd("D:/DATA AKMAL/IPB/Semester 5/Metode Peramalan Deret Waktu/Pertemuan 3")
IPM.papua <- as.data.frame(read_excel("Tugas Individu MPDW.xlsx", sheet = "Papua"))
IPM.papua$Tahun <- as.integer(IPM.papua$Tahun)
knitr::kable(IPM.papua, align = "c")
Tahun IPM
2010 54.45
2011 55.01
2012 55.55
2013 56.25
2014 56.75
2015 57.25
2016 58.05
2017 59.09
2018 60.06
2019 60.84
2020 60.44
2021 60.62

Eksplorasi Data

ggplot(IPM.papua, aes(x=Tahun, y=IPM)) + 
  ylim(54, 62) +
  geom_line(color="red", size=1) +
  geom_point(color="black", size=1.5) +
  geom_text(aes(label=IPM), size=3.2, hjust=0.7, vjust=-1.5) +
  labs(x="Tahun", y="Indeks Pembangunan Manusia",
       title="Indeks Pembangunan Manusia Provinsi Papua",
       subtitle=c("Periode 2010-2021"),
       caption="sumber: bps.go.id") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.subtitle = element_text(hjust = 0.5))

Berdasarkan plot di atas, dapat diamati bahwa secara umum IPM Provinsi Papua mengalami peningkatan yang stabil tiap tahunnya, menunjukkan bahwa kesejahteraan masyarakat Provinsi Papua terus meningkat setiap tahun. Penurunan IPM praktis hanya terjadi sekali, yaitu pada tahun 2020, akibat adanya pandemi COVID-19.

Pendeteksian Autokorelasi

Pemodelan Regresi Linear

Pendeteksian autokorelasi dilakukan pada model regresi linear antara peubah periode tahun dengan peubah IPM Provinsi Papua. Proses pemodelan regresi linear dapat diamati melalui rangkaian sintaks sebagai berikut.

model <- lm(IPM ~ Tahun, data = IPM.papua)
summary(model)
## 
## Call:
## lm(formula = IPM ~ Tahun, data = IPM.papua)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71718 -0.19083 -0.06851  0.11520  0.76604 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.215e+03  7.110e+01  -17.09 9.92e-09 ***
## Tahun        6.316e-01  3.528e-02   17.91 6.31e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4218 on 10 degrees of freedom
## Multiple R-squared:  0.9698, Adjusted R-squared:  0.9667 
## F-statistic: 320.6 on 1 and 10 DF,  p-value: 6.312e-09

Deteksi Autokorelasi dengan Metode Grafik

res <- residuals(model)
fitted <- predict(model)

par(mfrow = c(2,2))

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

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

hist(res, col = "steelblue",main = "Histogram Residual")

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

par(mfrow = c(1,2))

acf(res)
pacf(res)

Berdasarkan keenam plot residual di atas, secara eksploratif dapat diamati beberapa hal berikut.

  1. Melalui QQPlot, sisaan teramati menyebar menyerupai sebaran Normal.

  2. Teramati bahwa plot sisaan menyebar di sekitar nilai 0, mengindikasikan bahwa model telah memenuhi asumsi nilai harapan sisaan sama dengan 0.

  3. Plot sisaan teramati bergerak seperti dalam suatu tren yang mengindikasikan adanya ketidakbebasan sisaan.

Deteksi Autokorelasi dengan Statistik Uji

Selanjutnya, asumsi ini perlu diuji menggunakan statistik uji untuk menguji validitasnya dengan rangkaian hipotesis sebagai berikut.

H0 : Sisaan tidak mengalami autokorelasi

H1 : Sisaan mengalami autokorelasi

dwtest(model, alternative = 'two.sided') 
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.92878, p-value = 0.01062
## alternative hypothesis: true autocorrelation is not 0
bgtest(IPM ~ Tahun, data=IPM.papua, order=1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  IPM ~ Tahun
## LM test = 3.0355, df = 1, p-value = 0.08146
runs.test(resid(model), alternative = "two.sided")
## 
##  Runs Test - Two sided
## 
## data:  resid(model)
## Standardized Runs Statistic = -0.60553, p-value = 0.5448

Berdasarkan ketiga statistik uji di atas, Durbin-Watson test menghasilkan keputusan tolak H0 (p-value < 0.05). Dengan demikian, perlu dilakukan penanganan terhadap autokorelasi menggunakan dua metode, yaitu metode Chocrann-Ocutt dan Hildreth-Lu.

Penanganan Autokorelasi

Metode Cochrann-Ocutt

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.papua)
## 
##  number of interaction: 19
##  rho 0.651261
## 
## Durbin-Watson statistic 
## (original):    0.92878 , p-value: 5.308e-03
## (transformed): 1.14402 , p-value: 2.183e-02
##  
##  coefficients: 
## (Intercept)       Tahun 
## -1086.96736     0.56802
rho <- modelco$rho
rho
## [1] 0.6512606
y <- IPM.papua$IPM
y[-1]
##  [1] 55.01 55.55 56.25 56.75 57.25 58.05 59.09 60.06 60.84 60.44 60.62
y[-12]
##  [1] 54.45 55.01 55.55 56.25 56.75 57.25 58.05 59.09 60.06 60.84 60.44
y.tr <- y[-1]-y[-12]*rho
y.tr
##  [1] 19.54886 19.72416 20.07248 20.11659 20.29096 20.76533 21.28433 21.57701
##  [9] 21.72529 20.81731 21.25781
x <- IPM.papua$Tahun
x[-1]
##  [1] 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021
x[-12]
##  [1] 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020
x.tr <- x[-1]-x[-12]*rho
x.tr
##  [1] 701.9663 702.3150 702.6638 703.0125 703.3613 703.7100 704.0587 704.4075
##  [9] 704.7562 705.1049 705.4537

Selanjutnya, dilakukan pemodelan regresi linear dengan peubah x dan y hasil transformasi dengan rangkaian sintaks sebagai berikut.

model.coor <- lm(y.tr ~ x.tr)
summary(model.coor) 
## 
## Call:
## lm(formula = y.tr ~ x.tr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6278 -0.1518 -0.1134  0.2730  0.5281 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -379.0684    73.9294  -5.127 0.000622 ***
## x.tr           0.5680     0.1051   5.407 0.000429 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3843 on 9 degrees of freedom
## Multiple R-squared:  0.7646, Adjusted R-squared:  0.7384 
## F-statistic: 29.23 on 1 and 9 DF,  p-value: 0.0004292

Selanjutnya, dilakukan uji statistik untuk menguji validitasnya dengan rangkaian hipotesis sebagai berikut.

dwtest(model.coor, alternative = 'two.sided') 
## 
##  Durbin-Watson test
## 
## data:  model.coor
## DW = 1.144, p-value = 0.04367
## alternative hypothesis: true autocorrelation is not 0
bgtest(model.coor)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model.coor
## LM test = 1.7632, df = 1, p-value = 0.1842
runs.test(resid(model.coor), alternative = "two.sided")
## 
##  Runs Test - Two sided
## 
## data:  resid(model.coor)
## Standardized Runs Statistic = 0, p-value = 1

Dapat diamati bahwa uji asumsi menggunakan Breusch-Godfrey dan Runs test menghasilkan keputusan tak tolak H0 yang bermakna sisaan tidak mengalami autokorelasi.

Sementara itu, keputusan dari Durbin-Watson test masih menghasilkan tolak H0. Meskipun demikian, nilai p-value dari tes ini sudah lebih besar dibandingkan nilai p-value sebelum proses transformasi.

Selanjutnya, dilakukan proses transformasi balik pada persamaan regresi linear dengan rangkaian sintaks sebagai berikut.

cat("IPM = ", coef(model.coor)[1]/(1-rho), " + ", coef(model.coor)[2]," Tahun", sep = "")
## IPM = -1086.967 + 0.5680197 Tahun

Metode 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))
}

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)]
## [1] 0.651
par(mfrow = c(1,1))

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

Langkah selanjutnya adalah membuat model regresi linear dengan parameter rho 0.651 dengan rangkaian sintaks sebagai berikut.

model.hl <- hildreth.lu.func(0.651, model)
summary(model.hl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6277 -0.1519 -0.1134  0.2730  0.5281 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -379.4056    73.9293  -5.132 0.000618 ***
## x              0.5681     0.1050   5.412 0.000427 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3843 on 9 degrees of freedom
## Multiple R-squared:  0.7649, Adjusted R-squared:  0.7388 
## F-statistic: 29.29 on 1 and 9 DF,  p-value: 0.0004265

Selanjutnya, dilakukan uji statistik untuk menguji validitasnya dengan rangkaian hipotesis sebagai berikut.

dwtest(model.hl, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  model.hl
## DW = 1.1439, p-value = 0.04364
## alternative hypothesis: true autocorrelation is not 0
bgtest(model.hl)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model.hl
## LM test = 1.7638, df = 1, p-value = 0.1841
runs.test(resid(model.hl), alternative = "two.sided")
## 
##  Runs Test - Two sided
## 
## data:  resid(model.hl)
## Standardized Runs Statistic = 0, p-value = 1

Dapat diamati bahwa uji asumsi menggunakan Breusch-Godfrey dan Runs test menghasilkan keputusan tak tolak H0 yang bermakna sisaan tidak mengalami autokorelasi.

Sementara itu, keputusan dari Durbin-Watson test masih menghasilkan tolak H0. Meskipun demikian, nilai p-value dari tes ini sudah lebih besar dibandingkan nilai p-value sebelum proses transformasi.

Selanjutnya, dilakukan proses transformasi balik pada persamaan regresi linear dengan rangkaian sintaks sebagai berikut.

cat("IPM = ", coef(model.hl)[1]/(1-0.651), " + ", coef(model.hl)[2]," Tahun", sep = "")
## IPM = -1087.122 + 0.5680964 Tahun

Interpretasi Model Regresi Transformasi Balik

Berdasarkan dua metode yang telah digunakan, model regresi hasil transformasi balik cenderung sama (perbedaan tidak signifikan). Hal ini mengisyaratkan bahwa kedua metode menghasilkan dugaan model regresi linear yang serupa, yaitu b0 = -1087 dan b1 = 0.57. Dengan demikian, dapat disimpulkan bahwa kenaikan satu tahun meningkatkan dugaan rataan besar IPM Provinsi Papua sebesar 0.57.