Tugas Penanganan Autokorelasi
Fathoni Sabri
2025-09-04
Instalasi Package
library(dplyr)## Warning: package 'dplyr' was built under R version 4.5.1##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(TTR)## Warning: package 'TTR' was built under R version 4.5.1library(forecast)## Warning: package 'forecast' was built under R version 4.5.1## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(lmtest)## Warning: package 'lmtest' was built under R version 4.5.1## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(orcutt)
library(HoRM)## Warning: package 'HoRM' was built under R version 4.5.1Data
Data yang digunakan adalah Indeks Pembangunan Manusia di Provinsi Sumatera Barat pada Tahun 2010-2024
# Membuat vektor untuk tahun
tahun <- c(2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024)
# Membuat vektor IPM
IPM <- c(67.25, 67.81, 68.36, 68.91, 69.36, 69.98, 70.73, 71.24, 71.73, 72.39, 72.38, 72.65, 73.26, 73.75, 74.49)
# Menggabungkan kedua vektor untuk menjadi Matriks
data <- cbind(tahun, IPM)
# Mengubah format data menjadi data frame
data <- as.data.frame(data)
data## tahun IPM
## 1 2010 67.25
## 2 2011 67.81
## 3 2012 68.36
## 4 2013 68.91
## 5 2014 69.36
## 6 2015 69.98
## 7 2016 70.73
## 8 2017 71.24
## 9 2018 71.73
## 10 2019 72.39
## 11 2020 72.38
## 12 2021 72.65
## 13 2022 73.26
## 14 2023 73.75
## 15 2024 74.49Eksplorasi Data
Eksplorasi data dilakukan untuk melihat sebaran data dan menentukan langkah lanjutan yang tepat.
data.ts<-ts(data$IPM)
data.ts## Time Series:
## Start = 1
## End = 15
## Frequency = 1
## [1] 67.25 67.81 68.36 68.91 69.36 69.98 70.73 71.24 71.73 72.39 72.38 72.65
## [13] 73.26 73.75 74.49#Membuat plot time series
ts.plot(data.ts, xlab="Periode Waktu(Tahun)", ylab="IPM", main= "Plot Deret Waktu IPM")
points(data.ts)
Berdasarkan plot di atas, dapat dilihat bahwa data memiliki tren. sehingga, selanjutnya kita akan melakukan peramalan dan pemulusan dengan metode DMA dan DES. ## DMA
# Metode DMA (Double Moving Average)
dt.sma <- SMA(data.ts, n = 3) # setiap nilai = rata-rata 3 observasi terakhir
dma <- SMA(dt.sma, n = 3)
At <- 2 * dt.sma - dma # At = 2*SMA - DMA
Bt <- 2/(3 - 1) * (dt.sma - dma) # Bt = (2/(n-1)) * (SMA - DMA), dengan n=3
# Hitung fitted values (ramalan historis)
dt.dma <- At + Bt # Yhat_t = At + Bt
dt.ramal <- c(NA, dt.dma)
t <- 1:5
f <- c()
# Hitung ramalan ke depan (out-of-sample forecast)
# Rumus: F_{t+m} = A_t + B_t * m
for (i in t) {
f[i] <- At[length(At)] + Bt[length(Bt)] * i
}dt.gab <- cbind(
aktual = c(data.ts, rep(NA, 5)), # data asli (time series) + NA agar panjang sama
pemulusan1 = c(dt.sma, rep(NA, 5)), # hasil SMA (pemulusan pertama) + NA
pemulusan2 = c(dt.dma, rep(NA, 5)), # hasil DMA (pemulusan kedua) + NA
At = c(At, rep(NA, 5)), # komponen level + NA
Bt = c(Bt, rep(NA, 5)), # komponen trend + NA
ramalan = c(dt.ramal, f[-1]) # fitted values + ramalan ke depan
)
dt.gab## aktual pemulusan1 pemulusan2 At Bt ramalan
## [1,] 67.25 NA NA NA NA NA
## [2,] 67.81 NA NA NA NA NA
## [3,] 68.36 67.80667 NA NA NA NA
## [4,] 68.91 68.36000 NA NA NA NA
## [5,] 69.36 68.87667 69.93444 69.40556 0.5288889 NA
## [6,] 69.98 69.41667 70.48111 69.94889 0.5322222 69.93444
## [7,] 70.73 70.02333 71.19222 70.60778 0.5844444 70.48111
## [8,] 71.24 70.65000 71.89000 71.27000 0.6200000 71.19222
## [9,] 71.73 71.23333 72.42889 71.83111 0.5977778 71.89000
## [10,] 72.39 71.78667 72.91333 72.35000 0.5633333 72.42889
## [11,] 72.38 72.16667 73.04222 72.60444 0.4377778 72.91333
## [12,] 72.65 72.47333 73.13556 72.80444 0.3311111 73.04222
## [13,] 73.26 72.76333 73.35444 73.05889 0.2955556 73.13556
## [14,] 73.75 73.22000 74.02222 73.62111 0.4011111 73.35444
## [15,] 74.49 73.83333 74.95556 74.39444 0.5611111 74.02222
## [16,] NA NA NA NA NA 74.95556
## [17,] NA NA NA NA NA 75.51667
## [18,] NA NA NA NA NA 76.07778
## [19,] NA NA NA NA NA 76.63889
## [20,] NA NA NA NA NA 77.20000# Plot deret waktu ari data aktual
ts.plot(dt.gab[,1], # kolom 1 = data aktual
xlab = "Periode Waktu", # label sumbu x
ylab = "IPM", # label sumbu y
main = "DMA Data IPM N=3", # judul grafik
ylim = c(66, 78)) # batas sumbu y (dari 62 sampai 75)
# Tambahkan titik-titik (points) untuk data aktual
points(dt.gab[,1])
# Tambahkan titik-titik untuk data pemulusan 2 (kolom 3 = dt.dma)
points(dt.gab[,3])
# Tambahkan titik-titik untuk data ramalan (kolom 6)
points(dt.gab[,6])
# Tambahkan garis untuk data pemulusan (warna hijau, tebal garis = 2)
lines(dt.gab[,3], col="green", lwd=2)
# Tambahkan garis untuk data ramalan (warna merah, tebal garis = 2)
lines(dt.gab[,6], col="red", lwd=2)
# Tambahkan legenda di pojok kiri atas
legend("topleft",
c("data aktual", "data pemulusan", "data peramalan"), # keterangan
lty = 8, # jenis garis
col = c("black", "green", "red"), # warna sesuai garis
cex = 0.8) # ukuran teks
Selanjutnya lakukan perhitungan untuk mengetahui keakuratan dari metode
DMA
# 1. Hitung error (selisih data aktual dengan ramalan historis)
error.dma <- data.ts - dt.ramal[1:length(data.ts)]
# 2. Hitung Sum of Squared Errors (SSE)
SSE.dma <- sum(error.dma[6:length(data.ts)]^2)
# 3. Hitung Mean Squared Error (MSE)
MSE.dma <- mean(error.dma[6:length(data.ts)]^2)
# 4. Hitung Mean Absolute Percentage Error (MAPE)
MAPE.dma <- mean(abs((error.dma[6:length(data.ts)] / data.ts[6:length(data.ts)]) * 100))
# 5. Gabungkan hasil ke dalam matriks ringkasan
akurasi.dma <- matrix(c(SSE.dma, MSE.dma, MAPE.dma))
# 6. Tambahkan nama baris dan kolom
row.names(akurasi.dma) <- c("SSE", "MSE", "MAPE")
colnames(akurasi.dma) <- c("Akurasi m = 3")
# 7. Tampilkan hasil akurasi
akurasi.dma## Akurasi m = 3
## SSE 0.92246543
## MSE 0.09224654
## MAPE 0.33717481Didapati MAPE sebesar 0.337, artinya ramalan DMA dengan menggunakan m=3 meleset sekitar 0.337% dari nilai aslinya. ## DES Selanjutnya kita akan mencoba menggunakan metode DES (Double Exponential Smoothing)
# 1. Membagi data menjadi training dan testing
training <- data[1:11, 2] # ambil 11 observasi pertama (kolom ke-2 = variabel target)
testing <- data[12:15, 2] # ambil observasi ke-12 dan ke-15 sebagai data uji
# 2. Membuat objek time series
training.ts <- ts(training) # data latih dalam bentuk time series
testing.ts <- ts(testing, start=12) # data uji sebagai time series, dimulai dari periode ke-12Plot Data Asli
# Eksplorasi data (visualisasi seluruh data asli)
plot(data.ts,
col = "red", # warna merah untuk garis data
main = "Plot semua data", # judul grafik
xlab = "Periode", # label sumbu x
ylab = "Nilai") # label sumbu y
points(data.ts) # tambahkan titik di setiap observasi
### Plot Data Training
plot(training.ts, col="blue",main="Plot data training")
points(training.ts)
Selanjutnya melakukan pemulusan dengan DES sekaligus mencari nilai
lambda dan gamma optimum. Nilai lambda dan gamma optimum dapat dilihat
pada smoothing parameters alpha untuk nilai lambda, dan beta untuk nilai
gamma.
des.opt<- HoltWinters(training.ts, gamma = FALSE)
des.opt## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = training.ts, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0
## beta : 0
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 72.85
## b 0.56# Membuat plot hasil pemodelan Holt-Winters
# Garis hitam = data aktual (training), garis merah = hasil pemulusan model
plot(des.opt)
# Menambahkan legenda di pojok kiri atas
# Hitam = data aktual, Merah = hasil peramalan (fit values)
legend("topleft", c("Data Aktual", "Peramalan"),
col = c("black", "red"),
lty = c(1, 1))
# Membuat ramalan menggunakan model Holt-Winters hasil training
# h = 5 → meramalkan 5 periode ke depan
rdesopt <- forecast(des.opt, h = 5)
# Menampilkan hasil ramalan
# Output biasanya berisi: nilai ramalan (Point Forecast),
# serta interval prediksi (Lo 80/95 dan Hi 80/95)
rdesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 12 73.41 73.18173 73.63827 73.06089 73.75911
## 13 73.97 73.74173 74.19827 73.62089 74.31911
## 14 74.53 74.30173 74.75827 74.18089 74.87911
## 15 75.09 74.86173 75.31827 74.74089 75.43911
## 16 75.65 75.42173 75.87827 75.30089 75.99911Selanjutnya, kita menghitung akurasi dari metode DES.
# Mengambil nilai SSE (Sum of Squared Errors) dari model Holt-Winters pada data training
ssedes.train <- des.opt$SSE
# Menghitung MSE (Mean Squared Error) pada data training
# Caranya: SSE dibagi jumlah data training
msedes.train <- ssedes.train / length(training.ts)
# Mengambil nilai residual (selisih antara data aktual dan hasil model)
# Residual ini bisa dipakai untuk cek pola kesalahan model
sisaandes <- rdesopt$residuals
# Menampilkan beberapa nilai awal residual untuk dicek
head(sisaandes)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] NA NA -0.01 -0.02 -0.13 -0.07dan Menghitung MAPE nya.
# Rumus: rata-rata dari |(aktual - ramalan)/aktual| * 100
# Indeks mulai dari 3 → menghindari error pada awal data yang kurang stabil
mapedes.train <- sum(abs(sisaandes[3:length(training.ts)] / training.ts[3:length(training.ts)]) * 100) / length(training.ts)
# Membuat matriks berisi nilai akurasi (SSE, MSE, MAPE)
akurasides.opt <- matrix(c(ssedes.train, msedes.train, mapedes.train))
# Memberi nama baris matriks sesuai metrik akurasi
row.names(akurasides.opt) <- c("SSE", "MSE", "MAPE")
# Memberi nama kolom matriks sebagai keterangan model yang digunakan
colnames(akurasides.opt) <- c("Akurasi lamda dan gamma optimum")
# Menampilkan matriks akurasi
akurasides.opt## Akurasi lamda dan gamma optimum
## SSE 0.27250000
## MSE 0.02477273
## MAPE 0.12604711#Akurasi data testing
selisihdesopt<- rdesopt$mean-testing.ts
selisihdesopt## Time Series:
## Start = 12
## End = 15
## Frequency = 1
## [1] 0.76 0.71 0.78 0.60Menghitung SSE nya.
# selisihdesopt = perbedaan antara data aktual (testing) dan hasil ramalan
SSEtestingdesopt <- sum(selisihdesopt^2)
# Menghitung MSE (Mean Squared Error) pada data testing
# Caranya SSE dibagi jumlah data testing
SSEtestingdesopt <- SSEtestingdesopt / length(testing.ts)
# Menghitung MAPE (Mean Absolute Percentage Error) pada data testing
# Rumus: rata-rata dari |(aktual - ramalan)/aktual| * 100
MAPEtestingdesopt <- sum(abs(selisihdesopt / testing.ts) * 100) / length(testing.ts)
# Membuat matriks berisi nilai akurasi (SSE, MSE, MAPE) untuk data testing
akurasiDesTesting <- matrix(c(SSEtestingdesopt, SSEtestingdesopt, MAPEtestingdesopt))
# Memberi nama baris matriks sesuai metrik akurasi
row.names(akurasiDesTesting) <- c("SSE", "MSE", "MAPE")
# Memberi nama kolom matriks sesuai model yang digunakan
colnames(akurasiDesTesting) <- c("Akurasi lamda dan gamma optimum")
# Menampilkan matriks akurasi pada data testing
akurasiDesTesting## Akurasi lamda dan gamma optimum
## SSE 0.5125250
## MSE 0.5125250
## MAPE 0.9695917Perbandingan metode DMA dan DES
# di kiri adalah nilai untuk DMA dan di kanan adalah nilai untuk DES
cbind(akurasi.dma, akurasides.opt)## Akurasi m = 3 Akurasi lamda dan gamma optimum
## SSE 0.92246543 0.27250000
## MSE 0.09224654 0.02477273
## MAPE 0.33717481 0.12604711Berdasarkan nilai akurasi yang dihitung, terlihat bahwa nilai SSE, MSE, dan MAPE metode DES lebih kecil dibanding DMA. Maka dari itu, metode yang paling baik untuk melakukan peramalan pada data ini adalah metode DES. ## Nilai Korelasi
cor(tahun,IPM)## [1] 0.9951967Berdasarkan nilai korelasinya, periode waktu dan IPM memiliki hubungan positif yang sangat kuat
Regresi
#Pembuatan Model Regresi
#model regresi
model<- lm(IPM~tahun, data = data)
summary(model)##
## Call:
## lm(formula = IPM ~ tahun, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31210 -0.15735 -0.07974 0.15351 0.43262
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -942.30169 27.64400 -34.09 4.21e-14 ***
## tahun 0.50236 0.01371 36.65 1.65e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2293 on 13 degrees of freedom
## Multiple R-squared: 0.9904, Adjusted R-squared: 0.9897
## F-statistic: 1343 on 1 and 13 DF, p-value: 1.653e-14\[ \hat{y} = -942 + 0.5024 \times \text{tahun} \] Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki \(p\text{-value} < \alpha \ (5\%)\). Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model.
Selanjutnya dapat dilihat juga nilai
\[ R^2 = 0.9904 \] Artinya, sebesar 99.04% keragaman nilai IPM dapat dijelaskan oleh peubah tahun. Hasil ini menunjukkan hasil yang bagus, seolah mendapatkan hasil terbaik. Namun, kita perlu melakukan uji terhadap sisaannya seperti berikut ini.
#sisaan dan fitted value
sisaan<- residuals(model)
fitValue<- predict(model)#Diagnostik dengan eksploratif
par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "steelblue", lwd = 2)
plot(fitValue, sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(sisaan, col = "steelblue")
plot(seq(1,15,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,15,1), sisaan, col = "red")
abline(a = 0, b = 0, lwd = 2)
Dua plot di samping kiri digunakan untuk melihat apakah sisaan menyebar
normal. Normal Q-Q Plot di atas menunjukkan bahwa sisaan cenderung
menyebar normal, tetapi histogram dari sisaan tidak menunjukkan
demikian. Selanjutnya, dua plot di samping kanan digunakan untuk melihat
autokorelasi. Plot Sisaan vs Fitted Value dan Plot Sisaan vs Order
menunjukkan adanya pola pada sisaan. Untuk lebih lanjut akan digunakan
uji formal melihat normalitas sisaan dan plot ACF dan PACF untuk melihat
apakah ada autokorelasi atau tidak.
#Melihat Sisaan Menyebar Normal/Tidak
#H0: sisaan mengikuti sebaran normal
#H1: sisaan tidak mengikuti sebaran normal
shapiro.test(sisaan)##
## Shapiro-Wilk normality test
##
## data: sisaan
## W = 0.91441, p-value = 0.1582ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.17568, p-value = 0.6805
## alternative hypothesis: two-sidedBerdasarkan uji formal Saphiro-Wilk dan Kolmogorov-Smirnov didapatkan nilai p-value > α (5%). Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal.
Selanjutnya, menggunakan ACF dan PACF untuk melihat ada atau tidaknya autokorelasi
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)
Dari grafik, dapat dilihat ada autokorelasi yang cukup signifikan pada
Lag 1 ACF dan PACF. Untuk memastikan keberadaan autokorelasi. Akan
dilakukan uji Durbin-Watson.
Durbin Watson
#Deteksi autokorelasi dengan uji-Durbin Watson
#H0: tidak ada autokorelasi
#H1: ada autokorelasi
dwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 0.72476, p-value = 0.0004945
## alternative hypothesis: true autocorrelation is greater than 0Nilai p-value < 0.05 dapat disimpulkan bahwa tolak H0, cukup bukti mengatakan adanya autokorelasi. Oleh karena itu, diperlukan penangan autokorelasi. Pada data ini kita akan menggunakan dua metode, yaitu Cochrane-Orcutt dan Hildret-Lu.
Penanganan
Metode Cochrane-Orcutt
Penanganan metode Cochrane-Orcutt dapat dilakukan dengan bantuan packages Orcutt pada aplikasi R
modelCO<-cochrane.orcutt(model)
modelCO## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = IPM ~ tahun, data = data)
##
## number of interaction: 8
## rho 0.611472
##
## Durbin-Watson statistic
## (original): 0.72476 , p-value: 4.945e-04
## (transformed): 1.61195 , p-value: 1.366e-01
##
## coefficients:
## (Intercept) tahun
## -912.553376 0.487642Hasilnya adalah sebagai berikut \[ \hat{y} = -912.5534 + 0.4876 \times \text{tahun} \] Nilai p-value < 0.05, artinya cukup bukti untuk menyatakan bahwa sisaan terdapat autokorelasi pada taraf nyata 5%. Untuk nilai ρ̂ optimum yang digunakan adalah 0.558556 . Nilai tersebut dapat diketahui dengan syntax berikut.
#Rho optimum
rho<- modelCO$rho
rho## [1] 0.6114721#Transformasi Manual
IPM.trans<- IPM[-1]-IPM[-15]*rho
tahun.trans<- tahun[-1]-tahun[-15]*rho
modelCOmanual<- lm(IPM.trans~tahun.trans)
summary(modelCOmanual)##
## Call:
## lm(formula = IPM.trans ~ tahun.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34419 -0.07243 -0.01291 0.09642 0.25885
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -354.55248 23.90257 -14.83 4.42e-09 ***
## tahun.trans 0.48764 0.03047 16.00 1.85e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1786 on 12 degrees of freedom
## Multiple R-squared: 0.9552, Adjusted R-squared: 0.9515
## F-statistic: 256.1 on 1 and 12 DF, p-value: 1.848e-09Hasil model transformasi bukan merupakan model sesungguhnya. Koefisien regresi masih perlu dicari kembali mengikuti β∗0=β0+ρ̂ β0 dan β∗1=β1 .
b0bintang <- modelCOmanual$coefficients[-2]
b0 <- b0bintang/(1-rho)
b1 <- modelCOmanual$coefficients[-1]
b0## (Intercept)
## -912.5534b1## tahun.trans
## 0.4876416Hasil perhitungan menghasilkan angka yang hampir sama dengan menggunakan package.
Metode Hildreth-Lu
Metode ini akan mencari nilai SSE terkecil dan dapat dicari secara manual maupun menggunakan packages. Jika menggunakan packages, gunakan library packages HoRM.
#Penanganan Autokorelasi 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))
}
#Pencarian rho yang meminimumkan SSE
r <- c(seq(0.1,0.9, by= 0.1))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4)## rho SSE
## 1 0.1 0.5613
## 2 0.2 0.4983
## 3 0.3 0.4489
## 4 0.4 0.4132
## 5 0.5 0.3911
## 6 0.6 0.3827
## 7 0.7 0.3880
## 8 0.8 0.4069
## 9 0.9 0.4395Setelah nilai SSE diperolah menggunakan ρ dari 0.1 hingga 0.9. Selanjutnya kita mencari ρ dengan nilai SSE minimum, yaitu ρ=0.6. Namun, hasilnya masih kurang teliti. Untuk itu kita akan mencari nilai ρ yang minimum.
#Rho optimal di sekitar 0.4
rOpt <- seq(0.5,0.7, by= 0.0001)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model))}))
head(tabOpt[order(tabOpt$SSE),])## rho SSE
## 1116 0.6115 0.3826005
## 1115 0.6114 0.3826005
## 1117 0.6116 0.3826005
## 1114 0.6113 0.3826006
## 1118 0.6117 0.3826006
## 1113 0.6112 0.3826006Nilai SSE minimum adalah 0.3826005 dengan ρ=0.6115, selanjutnya kita akan membuat plotnya.
#Grafik SSE optimum
par(mfrow = c(1,1))
plot(tab$SSE ~ tab$rho , type = "l", xlab = "Rho", ylab = "SSE")
abline(v = tabOpt[tabOpt$SSE==min(tabOpt$SSE),"rho"], lty = 2, col="red",lwd=2)
text(x=0.6115 , y=0.3826005, labels = "rho=0.3826005 ", cex = 0.8)
Selanjutnya, model yang sudah didapatkan dievaluasi nilai ρ ke dalam
fungsi hildreth.lu.func, serta dilanjutkan dengan menguji autokorelasi
dengan uji Durbin-Watson. Setelah dilakukan pengecekan, persamaan
tersebut kemudian ditransformasikan kembali dan penjadi persamaan yang
sesungguhnya.
#Model terbaik
modelHL <- hildreth.lu.func(0.6115, model)
summary(modelHL)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34420 -0.07243 -0.01290 0.09641 0.25884
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -354.52619 23.90257 -14.83 4.42e-09 ***
## x 0.48764 0.03047 16.00 1.85e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1786 on 12 degrees of freedom
## Multiple R-squared: 0.9552, Adjusted R-squared: 0.9515
## F-statistic: 256.1 on 1 and 12 DF, p-value: 1.85e-09#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.559), "+", coef(modelHL)[2],"x", sep = "")## y = -803.9143+0.4876406x\[ \hat{y} = -803.9143 + 0.4876\times \text{tahun} \]
#Deteksi autokorelasi
dwtest(modelHL)##
## Durbin-Watson test
##
## data: modelHL
## DW = 1.612, p-value = 0.1367
## alternative hypothesis: true autocorrelation is greater than 0Hasil uji Durbin-Watson menunjukkan bawah nilai p-value sebesar 0.1367 , di mana p-value > α =5%. Artinya tak tolak H0 atau belum cukup bukti untuk menyatakan bahwa ada autokorelasi dalam data nilai IPM dengan metode Hildreth-Lu pada taraf nyata 5%.
Membandingkan Ketiga Metode
Selanjutnya kita akan membandingkan ketiga metode yang telah digunakan (metode awal, metode Cohcran-Orcutt, dan Hildreth-Lu)
#Perbandingan
sseModelawal <- anova(model)$`Sum Sq`[-1]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-1]
sseModelHL <- anova(modelHL)$`Sum Sq`[-1]
mseModelawal <- sseModelawal/length(IPM)
mseModelCO <- sseModelCO/length(IPM)
mseModelHL <- sseModelHL/length(IPM)
akurasi <- matrix(c(sseModelawal,sseModelCO,sseModelHL,
mseModelawal,mseModelCO,mseModelHL),nrow=2,ncol=3,byrow = T)
colnames(akurasi) <- c("Model Awal", "Model Cochrane-Orcutt", "Model Hildreth-Lu")
row.names(akurasi) <- c("SSE","MSE")
akurasi## Model Awal Model Cochrane-Orcutt Model Hildreth-Lu
## SSE 0.68373762 0.3826005 0.3826005
## MSE 0.04558251 0.0255067 0.0255067Berdasarkan hasil tersebut dapat diketahui bahwa penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu menghasilkan SSE dan MSE yang sama, yaitu 0.3826005 dan 0.0255067. Nilai ini lebih baik dibandingkan model awal sebelum penanganan autokorelasi, yang memiliki SSE 0.68373762 dan MSE 0.04558251.
Kesimpulan
Dari hasil analisis diketahui bahwa data IPM masih mengandung autokorelasi. Hal ini wajar, karena indikator penyusun IPM, khususnya yang berkaitan erat dengan aspek ekonomi, cenderung saling memengaruhi dari waktu ke waktu. Akibat adanya autokorelasi, model regresi yang dibangun menjadi kurang optimal, karena galat atau error model ikut membesar.
Pendeteksian autokorelasi dapat dilakukan dengan melihat pola residual melalui grafik (residual plot, ACF, dan PACF), maupun dengan uji formal seperti Durbin-Watson. Upaya perbaikan menggunakan metode Cochrane-Orcutt dan Hildreth-Lu berhasil menurunkan nilai SSE dan MSE dibandingkan model awal. Hasil uji Durbin–Watson pada model yang diperbaiki dengan metode Cochrane-Orcutt dan Hildreth-Lu menunjukkan bahwa autokorelasi juga sudah berhasil diatasi. Dengan demikian, asumsi bebas autokorelasi pada residual terpenuhi. Hal ini menjadikan model yang diperoleh dapat diandalkan untuk interpretasi dan prediksi, karena error tidak lagi saling berkorelasi. Selain itu, penurunan nilai SSE dan MSE dibandingkan model awal memperkuat bukti bahwa kedua metode perbaikan ini mampu meningkatkan kualitas model regresi.