Package dan Data

library(lmtest) 
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(orcutt)
library(HoRM)   
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
# Data IPM Kabupaten Kuningan (2004-2013)
t <- 2004:2013
ipm <- c(67.70, 68.53, 69.21, 69.70, 70.12, 70.42, 70.89, 71.55, 71.99, 72.47)
d <- data.frame(t, ipm)

# Lihat datanya
print(d)
##       t   ipm
## 1  2004 67.70
## 2  2005 68.53
## 3  2006 69.21
## 4  2007 69.70
## 5  2008 70.12
## 6  2009 70.42
## 7  2010 70.89
## 8  2011 71.55
## 9  2012 71.99
## 10 2013 72.47
plot(ts(ipm))
points(ipm, col="green")

Model Awal

model_awal <- lm(ipm ~ t, data = d)
summary(model_awal)
## 
## Call:
## lm(formula = ipm ~ t, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30200 -0.07750  0.00200  0.09417  0.20533 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -936.67000   35.87924  -26.11 4.98e-09 ***
## t              0.50133    0.01786   28.06 2.81e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1623 on 8 degrees of freedom
## Multiple R-squared:  0.9899, Adjusted R-squared:  0.9887 
## F-statistic: 787.6 on 1 and 8 DF,  p-value: 2.807e-09

Investigasi Sisaan

.

# Ambil sisaan dari model
sisaan_awal <- residuals(model_awal)

# Plot sisaan terhadap waktu
plot(d$t, sisaan_awal, type="o", pch=20, col="red",
     main="Plot Sisaan vs Waktu", xlab="Tahun", ylab="Sisaan")
abline(h=0, lty=2)

Plot di atas menunjukkan pola yang. Perhatikan bagaimana sisaan cenderung berkelompok: titik pertama di bawah nol, lalu naik ke atas nol, dan turun lagi. Ini adalah gejala visual yang kuat dari autokorelasi positif. Sisaan positif cenderung diikuti sisaan positif, dan sebaliknya.

Uji Durbin-Watson

Untuk uji formal, kita gunakan Uji Durbin-Watson.

dwtest(model_awal)
## 
##  Durbin-Watson test
## 
## data:  model_awal
## DW = 1.0332, p-value = 0.01238
## alternative hypothesis: true autocorrelation is greater than 0

P-value yang dihasilkan 0.01238, yaitu di bawah tingkat signifikansi 0.05 dan Tolak H0 yang berarti terbukti ada autokorelasi

Penanganan

Cochrane-Orcutt*

model_co <- cochrane.orcutt(model_awal)
summary(model_co)
## Call:
## lm(formula = ipm ~ t, data = d)
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -862.893466   31.992176 -26.972 2.469e-08 ***
## t              0.464640    0.015921  29.184 1.428e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0881 on 7 degrees of freedom
## Multiple R-squared:  0.9918 ,  Adjusted R-squared:  0.9907
## F-statistic: 851.7 on 1 and 7 DF,  p-value: < 1.428e-08
## 
## Durbin-Watson statistic 
## (original):    1.03315 , p-value: 1.238e-02
## (transformed): 1.59521 , p-value: 1.273e-01
# rho paling optimum
rho <- model_co$rho
cat("Rho optimum:", rho)
## Rho optimum: 0.2855958

Nilai rho yang diestimasi adalah yang digunakan untuk memperbaiki model. Lihat statistik Durbin-Watson baru. Nilainya sekarang jauh lebih dekat ke 2, dan p-value-nya (> 0.05) menunjukkan bahwa autokorelasi sudah berhasil diatasi. Perhatikan koefisien (Intercept) dan tahun yang baru. Ini adalah estimasi yang lebih valid.

ipm.trans <- d$ipm[-1] - d$ipm[-10]*rho

tahun.trans <- d$t[-1] - d$t[-10]*rho
model_co_manual <- lm(ipm.trans~tahun.trans)

b0_co_manual <- coef(model_co_manual)[1]/(1-rho)
b1_co_manual <- coef(model_co_manual)[2]
cat("Koefisien manual Cochrane-Orcutt:\n")
## Koefisien manual Cochrane-Orcutt:
cat("b0:", b0_co_manual, "\n")
## b0: -862.8935
cat("b1:", b1_co_manual, "\n")
## b1: 0.46464

Hasil koefisien manual kita sangat mirip dengan hasil dari fungsi otomatis. Ini membuktikan bahwa kita memahami mekanisme di baliknya. (Catatan: Perbedaan kecil mungkin terjadi karena fungsi otomatis melakukan beberapa iterasi).

Hildreth-Lu*

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

#Pencariab 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_awal))}))
round(tab, 4)
##   rho    SSE
## 1 0.1 0.0615
## 2 0.2 0.0559
## 3 0.3 0.0544
## 4 0.4 0.0571
## 5 0.5 0.0639
## 6 0.6 0.0749
## 7 0.7 0.0900
## 8 0.8 0.1093
## 9 0.9 0.1327

Pertama-tama akan dicari di mana kira-kira ρ yang menghasilkan SSE minimum. Pada hasil di atas terlihat ρ minimum ketika 0.3. Namun, hasil tersebut masih kurang teliti sehingga akan dicari kembali ρ yang lebih optimum dengan ketelitian yang lebih. Jika sebelumnya jarak antar ρ yang dicari adalah 0.1, kali ini jarak antar ρ adalah 0.001 dan dilakukan pada selang 0.1 sampai dengan 0.5.

#Rho optimal di sekitar 0.5
rOpt <- seq(0.1,0.5, by= 0.001)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model_awal))}))
head(tabOpt[order(tabOpt$SSE),])
##       rho        SSE
## 187 0.286 0.05433650
## 186 0.285 0.05433654
## 188 0.287 0.05433688
## 185 0.284 0.05433700
## 189 0.288 0.05433767
## 184 0.283 0.05433787
#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.286, y=0.05433650, labels = "rho=0.286", cex = 0.8)

model_hl <- hildreth.lu(d$ipm, d$t, rho = 0.286)
summary(model_hl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15305 -0.02398  0.01276  0.05069  0.09326 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -616.07312   22.85536  -26.95 2.48e-08 ***
## x              0.46462    0.01593   29.17 1.43e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0881 on 7 degrees of freedom
## Multiple R-squared:  0.9918, Adjusted R-squared:  0.9907 
## F-statistic: 850.6 on 1 and 7 DF,  p-value: 1.434e-08
#Transformasi Balik
cat("y = ", coef(model_hl)[1]/(1-0.286), "+", coef(model_hl)[2],"x", sep = "")
## y = -862.8475+0.4646172x
#Deteksi autokorelasi
dwtest(model_hl)
## 
##  Durbin-Watson test
## 
## data:  model_hl
## DW = 1.595, p-value = 0.1273
## alternative hypothesis: true autocorrelation is greater than 0

Plot yang dihasilkan secara visual menunjukkan nilai ρ mana yang memiliki SSE terendah. Hasil summary memberikan model akhir yang juga sudah “sehat”, dengan statistik D-W yang sudah membaik.

Evaluasi

# Menghitung SSE untuk setiap model secara manual
sse_awal <- sum(residuals(model_awal)^2)
sse_co <- sum(residuals(model_co)^2)
sse_hl <- sum(residuals(model_hl)^2)

# Membuat tabel perbandingan
data.frame(
  Metode = c("Model Awal (Sakit)", "Cochrane-Orcutt (Sehat)", "Hildreth-Lu (Sehat)"),
  SSE = c(sse_awal, sse_co, sse_hl),
  DW_Statistic = c(dwtest(model_awal)$statistic, model_co$DW[3], dwtest(model_hl)$statistic)
)
##                    Metode       SSE DW_Statistic
## 1      Model Awal (Sakit) 0.2106133     1.033154
## 2 Cochrane-Orcutt (Sehat) 0.3825796     1.595213
## 3     Hildreth-Lu (Sehat) 0.0543365     1.595034

Kesimpulan:

Hasil analisis menunjukkan bahwa model awal memiliki masalah autokorelasi yang cukup serius (DW ≈ 1.03). Setelah dilakukan koreksi dengan metode Cochrane-Orcutt dan Hildreth-Lu, nilai Durbin–Watson meningkat menjadi sekitar 1.59, mendekati 2, sehingga masalah autokorelasi berhasil diatasi.

Dari sisi ketepatan model, metode Hildreth-Lu memberikan hasil terbaik karena mampu menurunkan SSE secara signifikan dibandingkan metode lain, sekaligus menjaga kestabilan residual. Dengan demikian, model Hildreth-Lu dapat dipilih sebagai model akhir yang paling sehat dan dapat dipercaya untuk analisis selanjutnya.