Tugas Individu 1 Metode Peramalan Deret Waktu

Data

Data yang digunakan adalah data Indeks Pembangunan Manusia (IPM) di Bali dari tahun 2010-2021. Data tersebut bersumber dari Badan Pusat Statistik (BPS).

Package yang perlu di run :

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

Data yang dipilih adalah data Indeks Pembangunan Manusia (IPM) di Bali dengan tahun sebagai x dan IPM sebagai y.

#membuka file data
data_ipm <- read_excel("E:/KULIAH/SEMESTER 5/STA1341 Metode Peramalan Deret Waktu/TUGAS/Tugas Individu MPDW.xlsx", sheet = 2)
knitr::kable(data_ipm,align = "c")
Tahun IPM
2010 70.10
2011 70.87
2012 71.62
2013 72.09
2014 72.48
2015 73.27
2016 73.65
2017 74.30
2018 74.77
2019 75.38
2020 75.50
2021 75.69

Eksplorasi Data

x <- data_ipm$Tahun
y <- data_ipm$IPM
#diagram pencar identifikasi model
plot(x,y,pch = 20, col = "blue", main = "Scatter Plot IPM Di Bali",
     ylab = "IPM", xlab = "Tahun")

Berdasarkan scatter plot di atas, IPM di Bali mengalami kenaikan secara konstan dari tahun ke tahun.

#korelasi x dan y
cor(x,y)
## [1] 0.9923534

Model Regresi Deret Waktu

#model regresi
model <- lm(y~x, data = data_ipm)
summary(model)
## 
## Call:
## lm(formula = y ~ x, data = data_ipm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4800 -0.1125  0.0800  0.1725  0.2500 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -974.75000   41.22291  -23.65 4.15e-10 ***
## x              0.52000    0.02045   25.42 2.03e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2446 on 10 degrees of freedom
## Multiple R-squared:  0.9848, Adjusted R-squared:  0.9832 
## F-statistic: 646.4 on 1 and 10 DF,  p-value: 2.033e-10

Deteksi Autokorelasi

1. 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)

Secara eksploratif berdasarkan Plot Sissan vs Order terlihat bahwa plot membentuk pola tertentu. Maka asumsi kebebasan tidak terpenuhi, sehingga terdeteksi adanya gejala autokorelasi.

2. ACF dan PACF Plot

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

Berdasarkna Plot ACF dan PACF tidak ada yang signifikan. Sehingga dapat disimpulkan sisaan saling bebas.

Uji Statistik

1. Durbin Watson Test

H0: tidak ada autokorelasi

H1: ada autokorelasi

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

p-value < 5% : Tolak H0 sehingga cukup bukti untuk menyatakan bahwa terdapat autokorelasi (positif) pada taraf 5%.

2. Breusch-Godfrey Test

H0: tidak ada autokorelasi

H1: ada autokorelasi

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

p-value > 5% : Tak Tolak H0 sehingga tidak cukup bukti untuk menyatakan bahwa ada autokorelasi pada taraf 5% (tidak terdapat autokorelasi).

3. 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 = 0, p-value = 1

p-value > 5% : Tak Tolak H0 sehingga tidak cukup bukti untuk menyatakan bahwa ada autokorelasi pada taraf 5% (tidak terdapat autokorelasi).

Berdasarkan uji statistik, pada uji Durbin Watson didapatkan p-value < 0.05 yang artinya Tolak H0 (terdapat autokorelasi). Sementara itu, pada uji Breusch-Godfrey dan Runs test p-value > 0.05 yang artinya Terima H0 (tidak terdapat autokorelasi). Karenakan pada uji Durbin-Watson terdeteksi adanya autokorelasi, maka diperlukan adanya penanganan menggunakan Cochran Orcutt dan Hildreth-lu.

Penanganan Autokorelasi

1. Cochrane-Orcutt

#Interactive method using to solve first order autocorrelation problems.
modelco <- orcutt::cochrane.orcutt(model,convergence=5,max.iter = 1000)
modelco 
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x, data = data_ipm)
## 
##  number of interaction: 131
##  rho 0.735166
## 
## Durbin-Watson statistic 
## (original):    0.87546 , p-value: 3.527e-03
## (transformed): 2.12988 , p-value: 4.361e-01
##  
##  coefficients: 
## (Intercept)           x 
## -682.501632    0.375234
#rho optimum
rho <- modelco$rho
rho
## [1] 0.7351658
y
##  [1] 70.10 70.87 71.62 72.09 72.48 73.27 73.65 74.30 74.77 75.38 75.50 75.69
y[-1]
##  [1] 70.87 71.62 72.09 72.48 73.27 73.65 74.30 74.77 75.38 75.50 75.69
y[-12]
##  [1] 70.10 70.87 71.62 72.09 72.48 73.27 73.65 74.30 74.77 75.38 75.50
#transformasi terhadap y dan x
(y.trans <- y[-1]-y[-12]*rho)
##  [1] 19.33488 19.51880 19.43742 19.48190 19.98518 19.78440 20.15504 20.14718
##  [9] 20.41165 20.08320 20.18498
(x.trans <- x[-1]-x[-12]*rho)
##  [1] 533.3167 533.5816 533.8464 534.1112 534.3761 534.6409 534.9057 535.1706
##  [9] 535.4354 535.7002 535.9651
#model baru (model dengan peubah transformasi)
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho)
## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18523 -0.15405 -0.03413  0.13617  0.24765 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -180.74977   33.81058  -5.346 0.000465 ***
## x.trans        0.37523    0.06324   5.934 0.000220 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1757 on 9 degrees of freedom
## Multiple R-squared:  0.7964, Adjusted R-squared:  0.7738 
## F-statistic: 35.21 on 1 and 9 DF,  p-value: 0.0002198
#Melakukan pengecekan autokorelasi

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

p-value > 5% : Tak Tolak H0 sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5% (tidak terdapat autokorelasi).

Transformasi Balik

#Trasformasi Balik

cat("IPM = ", coef(modelcorho)[1]/(1-rho), "+", coef(modelcorho)[2],"Tahun", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal
## IPM = -682.5016+0.3752344Tahun

2. 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,99 memiliki SSE Terkecil
##     rho    SSE
## 1  0.10 0.3873
## 2  0.20 0.3555
## 3  0.30 0.3291
## 4  0.40 0.3082
## 5  0.50 0.2927
## 6  0.60 0.2827
## 7  0.70 0.2780
## 8  0.80 0.2788
## 9  0.90 0.2851
## 10 0.91 0.2860
## 11 0.92 0.2869
## 12 0.93 0.2880
## 13 0.94 0.2891
## 14 0.95 0.2902
## 15 0.96 0.2914
## 16 0.97 0.2926
## 17 0.98 0.2939
## 18 0.99 0.2953
tab$rho[which.min(tab$SSE)]#rho optimal
## [1] 0.7
#grafik rho dan SSE
par(mfrow = c(1,1))
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

#rho optimal di sekitar 0.7
r <- seq(0.6, 0.8, by= 0.005)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0,735 memiliki SSE Terkecil
##      rho    SSE
## 1  0.600 0.2827
## 2  0.605 0.2823
## 3  0.610 0.2820
## 4  0.615 0.2816
## 5  0.620 0.2813
## 6  0.625 0.2810
## 7  0.630 0.2807
## 8  0.635 0.2804
## 9  0.640 0.2802
## 10 0.645 0.2799
## 11 0.650 0.2797
## 12 0.655 0.2794
## 13 0.660 0.2792
## 14 0.665 0.2790
## 15 0.670 0.2789
## 16 0.675 0.2787
## 17 0.680 0.2785
## 18 0.685 0.2784
## 19 0.690 0.2782
## 20 0.695 0.2781
## 21 0.700 0.2780
## 22 0.705 0.2779
## 23 0.710 0.2779
## 24 0.715 0.2778
## 25 0.720 0.2778
## 26 0.725 0.2777
## 27 0.730 0.2777
## 28 0.735 0.2777
## 29 0.740 0.2777
## 30 0.745 0.2777
## 31 0.750 0.2778
## 32 0.755 0.2778
## 33 0.760 0.2779
## 34 0.765 0.2779
## 35 0.770 0.2780
## 36 0.775 0.2781
## 37 0.780 0.2782
## 38 0.785 0.2784
## 39 0.790 0.2785
## 40 0.795 0.2787
## 41 0.800 0.2788
tab$rho[which.min(tab$SSE)]#rho optimal
## [1] 0.735
#grafik SSE optimum
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

# Model terbaik (model dengan peubah transformasi)
modelhl <- hildreth.lu.func(0.735, model)
summary(modelhl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18521 -0.15406 -0.03417  0.13618  0.24766 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -180.9202    33.8105  -5.351 0.000462 ***
## x              0.3753     0.0632   5.939 0.000218 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1757 on 9 degrees of freedom
## Multiple R-squared:  0.7967, Adjusted R-squared:  0.7741 
## F-statistic: 35.27 on 1 and 9 DF,  p-value: 0.0002183
# Melakukan pengecekan autokorelasi

lmtest::dwtest(modelhl)
## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 2.1296, p-value = 0.4359
## alternative hypothesis: true autocorrelation is greater than 0

p-value > 5% : Tak Tolak H0 sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5% (tidak terdapat autokorelasi).

Transformasi Balik

#Transformasi Balik

cat("IPM = ", coef(modelhl)[1]/(1-0.735), "+", coef(modelhl)[2],"Tahun", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal
## IPM = -682.7179+0.3753413Tahun