Tahap Pra-proses data

Sebelum melakukan eksplorasi, maka kita perlu memanggil library packages yang akan dipakai.

library(tseries)
library(forecast)
library(TTR)
library(readxl)
library(orcutt)
library(HoRM)
library(rmarkdown)
library(lmtest)

Pengenalan Data

IPM atau Indeks Pembangunan Manusia merupakan indikator penting untuk mengukur keberhasilan dalam upaya membangun kualitas hidup manusia (masyarakat/penduduk). Bagi Indonesia, IPM merupakan data strategis karena selain sebagai ukuran kinerja Pemerintah, IPM juga digunakan sebagai salah satu alokator penentuan Dana Alokasi Umum (BPS, 2021).

Olah data

Kemudian, input data yang digunakan, yakni IPM Indonesia periode tahun 2010 - 2021 (n=12)

data_p3 <- read_excel("C:/Users/Dicky Girsang/Desktop/Semester 5/STA1341 Metode Peramalan Deret Waktu/Pertemuan 3/Tugas Individu MPDW 3.xlsx")
data_ts <- ts(data_p3$y)

Eksplorasi Data

par(mfrow=c(1,1))
plot(data_ts, col="red",main="IPM Indonesia Periode tahun 2010 - 2021",
     xlab = "Periode Tahun", ylab= "IPM Indonesia")
points(data_ts)

x_data <- data_p3$tahun
model <- lm(data_ts~x_data, data_ts)
summary(model)
## 
## Call:
## lm(formula = data_ts ~ x_data, data = data_ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47192 -0.13823 -0.00544  0.20485  0.28867 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.046e+03  4.162e+01  -25.13 2.28e-10 ***
## x_data       5.535e-01  2.065e-02   26.80 1.21e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2469 on 10 degrees of freedom
## Multiple R-squared:  0.9863, Adjusted R-squared:  0.9849 
## F-statistic: 718.5 on 1 and 10 DF,  p-value: 1.207e-10

Variabel data

x <- data_p3$tahun
y <- data_ts

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 = "darkorange",main = "Histogram Residual")

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)

Dari masing-masing plot di atas, diperoleh hasil bahwa terdapat kecenderungan data menyebar tidak normal, secara observatif sebagian data menyebar tidak di sekitar garis qqplot. Terdapat juga indikasi bahwa data yang digunakan mengikuti pola tertentu. Dengan kata lain, terjadi heteroskedasitas pada data tersebut. Dengan tujuan memenuhi asumsi kebebasan klasik, maka perlu dilakukan langkah untuk mendeteksi ada atau tidaknya autokorelasi pada data yang digunakan.

Identifikasi Autokorelasi

ACF dan PACF

par(mfrow = c(1,2))
acf(resi1)
pacf(resi1)

Durbin-Watson dan Breusch-Godfrey Test

H0: Tidak terdapat autokorelasi pada data

H1: Terdapat autokorelasi pada data

lmtest::dwtest(model, alternative = 'two.sided') 
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.5839, p-value = 0.0003055
## alternative hypothesis: true autocorrelation is not 0
lmtest::bgtest(y ~ x, order=1) 
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  y ~ x
## LM test = 6.254, df = 1, p-value = 0.01239

Dari kedua uji di atas, diperoleh hasil bahwa P-value < 0.5, sehingga disimpulkan bahwa terdapat autokorelasi pada data yang digunakan

Penanganan Autokorelasi

Cochrane Orcutt

modelco <- orcutt::cochrane.orcutt(model,convergence=6,max.iter = 10000)
rho <- modelco$rho
rho
## [1] 0.9767815
y[c(-1, -12)]
##  [1] 67.09 67.70 68.31 68.90 69.55 70.18 70.81 71.39 71.92 71.94
(y.trans <- y[-1]-y[-12]*rho)
##  [1] 2.104726 2.167729 2.181892 2.176055 2.249754 2.244846 2.259474 2.224102
##  [9] 2.187568 1.689874 2.020338
(x.trans <- x[-1]-x[-12]*rho)
##  [1] 47.66917 47.69239 47.71561 47.73883 47.76205 47.78526 47.80848 47.83170
##  [9] 47.85492 47.87814 47.90136

Model hasil transformasi

model_trans <- lm(y.trans~x.trans)
summary(model_trans) 
## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36669 -0.03286 -0.00108  0.10941  0.14263 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  43.4951    30.9662   1.405    0.194
## x.trans      -0.8655     0.6480  -1.336    0.214
## 
## Residual standard error: 0.1578 on 9 degrees of freedom
## Multiple R-squared:  0.1654, Adjusted R-squared:  0.07268 
## F-statistic: 1.784 on 1 and 9 DF,  p-value: 0.2145
lmtest::dwtest(model_trans,alternative = 'two.sided') 
## 
##  Durbin-Watson test
## 
## data:  model_trans
## DW = 1.6509, p-value = 0.3095
## alternative hypothesis: true autocorrelation is not 0

Setelah dilakukan penanganan autokorelasi menggunakan dw-test (seperti pada contoh di atas), p-value (0.3095) > 0.05, artinya autokorelasi pada data sudah berhasil diatasi.

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.5185
## 2  0.20 0.4621
## 3  0.30 0.4116
## 4  0.40 0.3669
## 5  0.50 0.3281
## 6  0.60 0.2952
## 7  0.70 0.2682
## 8  0.80 0.2471
## 9  0.90 0.2318
## 10 0.91 0.2306
## 11 0.92 0.2295
## 12 0.93 0.2284
## 13 0.94 0.2274
## 14 0.95 0.2264
## 15 0.96 0.2255
## 16 0.97 0.2247
## 17 0.98 0.2239
## 18 0.99 0.2231
#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)

Diperoleh hasil bahwa rho yang memiliki nilai SSE terkecil adalah 0.99 (ketelitian 0.01)

#rho optimal di sekitar 0.99
r <- seq(0.9, 1, by= 0.001)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0.999 memiliki SSE Terkecil
##       rho    SSE
## 1   0.900 0.2318
## 2   0.901 0.2317
## 3   0.902 0.2316
## 4   0.903 0.2315
## 5   0.904 0.2314
## 6   0.905 0.2312
## 7   0.906 0.2311
## 8   0.907 0.2310
## 9   0.908 0.2309
## 10  0.909 0.2308
## 11  0.910 0.2306
## 12  0.911 0.2305
## 13  0.912 0.2304
## 14  0.913 0.2303
## 15  0.914 0.2302
## 16  0.915 0.2301
## 17  0.916 0.2300
## 18  0.917 0.2298
## 19  0.918 0.2297
## 20  0.919 0.2296
## 21  0.920 0.2295
## 22  0.921 0.2294
## 23  0.922 0.2293
## 24  0.923 0.2292
## 25  0.924 0.2291
## 26  0.925 0.2290
## 27  0.926 0.2288
## 28  0.927 0.2287
## 29  0.928 0.2286
## 30  0.929 0.2285
## 31  0.930 0.2284
## 32  0.931 0.2283
## 33  0.932 0.2282
## 34  0.933 0.2281
## 35  0.934 0.2280
## 36  0.935 0.2279
## 37  0.936 0.2278
## 38  0.937 0.2277
## 39  0.938 0.2276
## 40  0.939 0.2275
## 41  0.940 0.2274
## 42  0.941 0.2273
## 43  0.942 0.2272
## 44  0.943 0.2271
## 45  0.944 0.2270
## 46  0.945 0.2269
## 47  0.946 0.2268
## 48  0.947 0.2267
## 49  0.948 0.2266
## 50  0.949 0.2265
## 51  0.950 0.2264
## 52  0.951 0.2263
## 53  0.952 0.2262
## 54  0.953 0.2261
## 55  0.954 0.2261
## 56  0.955 0.2260
## 57  0.956 0.2259
## 58  0.957 0.2258
## 59  0.958 0.2257
## 60  0.959 0.2256
## 61  0.960 0.2255
## 62  0.961 0.2254
## 63  0.962 0.2253
## 64  0.963 0.2253
## 65  0.964 0.2252
## 66  0.965 0.2251
## 67  0.966 0.2250
## 68  0.967 0.2249
## 69  0.968 0.2248
## 70  0.969 0.2247
## 71  0.970 0.2247
## 72  0.971 0.2246
## 73  0.972 0.2245
## 74  0.973 0.2244
## 75  0.974 0.2243
## 76  0.975 0.2243
## 77  0.976 0.2242
## 78  0.977 0.2241
## 79  0.978 0.2240
## 80  0.979 0.2240
## 81  0.980 0.2239
## 82  0.981 0.2238
## 83  0.982 0.2237
## 84  0.983 0.2237
## 85  0.984 0.2236
## 86  0.985 0.2235
## 87  0.986 0.2234
## 88  0.987 0.2234
## 89  0.988 0.2233
## 90  0.989 0.2232
## 91  0.990 0.2231
## 92  0.991 0.2231
## 93  0.992 0.2230
## 94  0.993 0.2229
## 95  0.994 0.2229
## 96  0.995 0.2228
## 97  0.996 0.2227
## 98  0.997 0.2227
## 99  0.998 0.2226
## 100 0.999 0.2225
## 101 1.000 0.3463

Setelah dilakukan iterasi pada sekitar nilai 0.99 untuk mendapatkan hasil dengan ketelitian yang lebih tinggi, bahwa rho sebesar 0.999 memiliki nilai SSE terkecil.

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

modelhl <- hildreth.lu.func(0.999, model)
summary(modelhl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36934 -0.03111 -0.00074  0.10680  0.14003 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    99.99      45.20   2.212   0.0543 .
## x             -32.97      14.99  -2.199   0.0554 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1572 on 9 degrees of freedom
## Multiple R-squared:  0.3495, Adjusted R-squared:  0.2772 
## F-statistic: 4.835 on 1 and 9 DF,  p-value: 0.05545
lmtest::dwtest(modelhl)
## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 1.6973, p-value = 0.1758
## alternative hypothesis: true autocorrelation is greater than 0

Sebagaii hasil akhir dari penanganan autokorelasi di atas, diperoleh p-value (0.1758) > 0.05, artinya terima H0, yakni tidak terdapat autokorelasi pada model data yang baru.

Model Akhir

Setelah dilakukan penanganan, diperoleh model seperti di bawah ini:

cat("Indeks Pembangunan Manusia = ", coef(model_trans)[1]/(1-0.999), "+", coef(model_trans)[2]," Tahun", sep = "") 
## Indeks Pembangunan Manusia = 43495.11+-0.8655004 Tahun