1. Mencari Korelasi

round(cor(x,y), 2)

## [1] 0.99

2. Membuat Scatterplot

#diagram pencar identifikasi model
plot(x,y,pch = 20, col = "blue", main = "Scatter Plot Tahun vs Nilai IPM", ylab = "Nilai IPM", xlab = "Tahun")

Dapat dilihat dengan jelas dari scatterplot di atas bahwa terdapat hubungan positif antara variabel tahun dan IPM. Nilai korelasi antara kedua variabel ini juga hampir sempurna, yaitu sebesar \(0.99\). Jika dianalisis lebih lanjut, plot di atas juga menunjukkan kenaikan nilai IPM yang cenderung lebih rendah pada tahun 2020 dan 2021, yaitu ketika COVID-19 menyerang Indonesia.

3. Membuat Model Regresi Deret Waktu

Model regresi yang digunakan adalah \(Y=β0+β1x\)

model <- lm(y~x, data = dataregresi)
summary(model)

## 
## Call:
## lm(formula = y ~ x, data = dataregresi)
## 
## 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

Berdasarkan hasil di atas, didapat model regresi linear adalah: \[ Y = -974.75 + 0.52x \] Artinya nilai dugaan rataan Y bertambah sebesar 0.52 jika X bertambah satu satuan atau bertambah satu tahun lamanya. Berdasarkan uji-t parsial untuk intersep dan koefisien X, terlihat bahwa kedua parameter berpengaruh nyata pada taraf 5%. Hal ini juga sejalan dengan hasil uji F yang dilakukan didapat bahwa \(P−Value<0.05\), sehingga dapat disimpulkan bahwa terdapat minimal 1 variabel yang berpengaruh nyata secara signifikan pada taraf 5%. Nilai \(R^2=0.9848\) menujukkan bahwa variabel x (tahun) dapat menjelaskan keragaman variabel y (IPM) sebesar 98.48%.

1. Plot Residual vs Order

#Residual 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(model, which=1)

hist(resi1, col = "steelblue")

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

1. Normal Q-Q Plot

   Berdasarkan hasil eksplorasi di atas, terlihat bahwa banyak amatan yang telah mendekati garis qq-plot distribusi normal. Sehingga secara eksploratif, dapat disimpulkan bahwa sisaan cukup menyebar normal.

2. Residual vs Fitted Values

   Berdasarkan eksplorasi plot Residuals vs Fitted di atas, data sisaan cukup menyebar di sekitar titik 0, sehingga asumsi kehomogenan ragam masih terpenuhi.

3. Histogram

   Berdasarkan histogram di atas, data terlihat cukup menjulur ke kiri.

4. Residual vs Order

   Berdasarkan hasil eksplorasi di atas, titik pada plot kebebasan sisaan mayoritas bergerak di sekitar titik nol. Namun, terdapat beberapa titik amatan yang terletak cukup jauh dari titik nol. Sehingga, belum dapat disimpulkan apakah terdapat autokorelasi atau tidak.

2. Plot ACF PACF

Autocorrelation Function (ACF) merupakan korelasi antar data di dalam time series yang dipisahkan dalam lag. Lag adalah beda urutan suatu sisaan dengan sisaan sebelumnya. Artinya, jika lag =2, maka nilai t akan dikorelasikan dengan t+2 (artinya data pertama akan dikorelasikan dengan data ketiga), jika menggunakan lag=4, maka kemudian t akan dikorelasikan dengan t+4, begitu seterusnya. Hal ini biasanya ditulis dengan rumus korelasi (X_t, X_k+t) dengan k adalah jumlah lag dan t adalah nilainya.

Sama halnya autocorrelation, Partial Autocorrelation Function (PACF) juga menghitung korelasi antar data di time series. Perbedaannya adalah jika ACF mencari korelasi antara X_t dan X_t+k ( k = lag), maka PACF juga menghitung korelasi antara X_t dan X_t+k namun sebelumnya menghilangkan data antara X_t dan X_t+k terlebih dahulu. Berarti data yang dihilangkan adalah data X_t+1 sampai dengan X_t+k-1

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

Berdasarkan kedua plot di atas, terlihat bahwa garis pada plot tidak melebihi garis batas biru. Hal ini mengindikasikan tidak adanya autokorelasi.

3. Uji Durbin-Watson

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

1. Cochrane-Orcutt

Metode Cochrane-Orcutt dilakukan dengan menghitung nilai koefisien autokorelasi (\(\rho\)) menggunakan nilai error pada model regresi. Pada metode ini autokorelasi dihilangkan secara bertahap dari bentuk yang paling sederhana sehingga autokorelasi dapat diatasi. Setelah ditemukan nilai duga bagi \(\rho\), maka langkah selanjutnya adalah melakukan transformasi sebagai berikut:

\[ y_t^*=y_t-\hat{p}y_{t-1};x_t^*=z_t-\hat{p}x_{t-1} \]

#Interactive method using to solve first order autocorrelation problems. 
modelco <- cochrane.orcutt(model, convergence = 7, max.iter = 1000)
modelco 

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x, data = dataregresi)
## 
##  number of interaction: 234
##  rho 0.73534
## 
## Durbin-Watson statistic 
## (original):    0.87546 , p-value: 3.527e-03
## (transformed): 2.13019 , p-value: 4.363e-01
##  
##  coefficients: 
## (Intercept)           x 
## -682.274737    0.375122

#rho optimum
rho <- modelco$rho
rho

## [1] 0.7353396

Telah didapat nilai (\(\rho\)) optimum sebesar \(0.735\). Selanjutnya akan dilakukan transformasi peubah X dan Y menggunakan nilai (\(\rho\)).

#transformasi terhadap y dan x
(y.trans <- y[-1]-y[-12]*rho)

##  [1] 19.32270 19.50649 19.42498 19.46937 19.97259 19.77167 20.14224 20.13427
##  [9] 20.39866 20.07010 20.17186

(x.trans <- x[-1]-x[-12]*rho)

##  [1] 532.9675 533.2321 533.4968 533.7615 534.0261 534.2908 534.5554 534.8201
##  [9] 535.0848 535.3494 535.6141

#model baru
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho) 

## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18524 -0.15403 -0.03408  0.13616  0.24765 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -180.57113   33.81062  -5.341 0.000468 ***
## x.trans        0.37512    0.06328   5.928 0.000221 ***
## ---
## 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.7961, Adjusted R-squared:  0.7734 
## F-statistic: 35.14 on 1 and 9 DF,  p-value: 0.0002213

lmtest::dwtest(modelcorho,alternative = 'two.sided') 

## 
##  Durbin-Watson test
## 
## data:  modelcorho
## DW = 2.1302, p-value = 0.8726
## alternative hypothesis: true autocorrelation is not 0

Berdasarkan hasil di atas, didapat nilai \(P-Value = 0.8726\). Artinya pada taraf nyata 5%, ada bukti untuk menyatakan bahwa tidak ada autokorelasi atau terbukti bahwa autokorelasi telah berhasil ditangani dengan metode Cochrane-Orcutt. Selanjutnya akan dilakukan transformasi untuk mendapatkan model regresi deret waktu terbaru.

b0 <- modelcorho$coefficients[1]/(1-rho)
b0

## (Intercept) 
##   -682.2747

b1 <- modelcorho$coefficients[2]
b1

##   x.trans 
## 0.3751221

Setelah dilakukan transformasi balik, didapat model regresi deret waktu yaitu: \[ \hat{y_t} = -682,27 + 0.375x_t \]

2. Hildreth-lu

Prosedur Hildreth-Lu menggunakan pendekatan yang sama untuk memperkirakan parameter autokorelasi ρ untuk digunakan dalam transformasi. Nilai ρ dipilih dengan prosedur Hildreth-Lu, yang meminimalkan sum square error untuk model regresi yang ditransformasi.

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

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

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

Pada iterasi pertama dapat terlihat bahwa \(ρ\) optimum yang menghasilkan SSE terkecil berada di sekitar \(0.7\). Sehingga selang kita perkecil.

#rho optimal di sekitar 0.7
r <- seq(0.68, 0.75, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4)

##    rho    SSE
## 1 0.68 0.2785
## 2 0.69 0.2782
## 3 0.70 0.2780
## 4 0.71 0.2779
## 5 0.72 0.2778
## 6 0.73 0.2777
## 7 0.74 0.2777
## 8 0.75 0.2778

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

Dapat dilihat bahwa nilai \(ρ\) optimum yang menghasilkan SSE terkecil adalah sebesar \(0.74\). Sehingga model optimum yang dapat dibuat dengan metode Hildreth-Lu adalah sebagai berikut:

# Model terbaik
modelhl <- hildreth.lu.func(0.74, model)
summary(modelhl)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18566 -0.15366 -0.03286  0.13587  0.24747 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -175.77964   33.81212  -5.199 0.000565 ***
## x              0.37206    0.06442   5.776 0.000267 ***
## ---
## 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.7875, Adjusted R-squared:  0.7639 
## F-statistic: 33.36 on 1 and 9 DF,  p-value: 0.0002675

# Deteksi autokorelasi
dwtest(modelhl)

## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 2.1385, p-value = 0.4421
## alternative hypothesis: true autocorrelation is greater than 0

Pada uji Durbin-Watson, dapat dilihat bahwa nilai \(P−value = 0.4421 > 0.05\). Sehingga tak tolak H0, atau belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf nyata 5%.

# Transformasi Balik
cat("IPM = ", coef(modelhl)[1]/(1-0.74), "+", coef(modelhl)[2],"tahun", sep = "")

## IPM = -676.0755+0.3720559tahun

Model regresi setelah melakukan penanganan autokorelasi adalah:

\[ \hat{y_t}= = -676.075 + 0.372x_t \]