Dalam prosesnya dibutuhkan beberapa package yang bisa digunakan antara lain :
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
library(lawstat)
##
## Attaching package: 'lawstat'
## The following object is masked from 'package:tseries':
##
## runs.test
IPM menjelaskan bagaimana penduduk dapat mengakses hasil pembangunan dalam memperoleh pendapatan, kesehatan, pendidikan, dan sebagainya. Adapun dimensi dasar IPM terdiri dari 3 point yaitu Umur panjang dan hidup sehat, Pengetahuan, Standar hidup layak. Salah satu manfaat IPM yaitu sebagai indikator penting untuk mengukur keberhasilan dalam upaya membangun kualitas hidup manusia (masyarakat/penduduk). (BPS 2020)
Berikut data yang digunakan merupakan data Indeks Pembangunan Sumatera Selatan pada Tahun 2010 - 2021
Tugas_Individu_MPDW <- read_excel("C:/Users/asus/OneDrive/Documents/MK Semester 5/Metode Peramalan Deret Waktu/Tugas Individu MPDW.xlsx",
sheet = "Sheet1")
Peubah x kita misalkan sebagai kolom dalam tahun dan Peubah y sebagai kolom IPM (Indeks pembangunan manusia)
x <- Tugas_Individu_MPDW$Tahun
y <- Tugas_Individu_MPDW$IPM
#korelasi x dan y
cor(x,y)
## [1] 0.9913599
#diagram pencar identifikasi model
plot(x,y,pch = 20, col = "blue", main = "Scatter Plot X vs Y",
ylab = "Nilai Peubah Y", xlab = "Nilai Peubah X")
Dapat dilihat pada scatter plot di atas bahwa terdapat korelasi yang kuat dan positif antara peubah x dan y, dengan nilai korelasinya yaitu 0.9932939.
Selanjutnya dilakukan pendugaan model regresi awal dan didapatkan persamaan sebagai berikut :
#model regresi
model <- lm(y~x, data = Tugas_Individu_MPDW)
summary(model)
##
## Call:
## lm(formula = y ~ x, data = Tugas_Individu_MPDW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.55667 -0.15076 -0.01212 0.25902 0.34697
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.065e+03 4.738e+01 -22.47 6.85e-10 ***
## x 5.618e-01 2.351e-02 23.90 3.74e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2811 on 10 degrees of freedom
## Multiple R-squared: 0.9828, Adjusted R-squared: 0.9811
## F-statistic: 571.2 on 1 and 10 DF, p-value: 3.738e-10
Model regresi awal yang diperoleh adalah y = -1065 + 0.5618x
#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(1:12, resi1, col = "steelblue", pch = 20,
xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(1:12, resi1, col = "red")
abline(a = 0, b = 0, lwd = 2)
Pengujian pertama dilakukan secara eksploratif dengan melihat Normal QQ plot, plot sisaan vs fitted values, histogram, dan plot sisaan vs order, dari plot diatas yaitu normal QQ plot dan histogram dapat kita lihat datanya berbentuk linear dan cenderung menjulur ke kiri. Dan pada plot sisaan vs fitted value dan plot sisaan vs order, dapat kita lihat bahwasannya data tersebut tidak terlihat pola khusus sehingga diduga tidak terdapat autokorelasi antar amatan, sehingga butuh dilakukan pengujian lebih lanjut yaitu dengan melihat nilai ACF dan PACF
#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(2,1))
acf(resi1)
pacf(resi1)
Dari plot diatas yaitu ACF dan PACF terlihat tidak ada garis yang melewati garis biru sehingga dari sini dapat diketahui tidak ada autokorelasi di kedua plot tersebut.
Selanjutnya dilakukan uji formal durbin watson untuk memastikan kembali
H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::dwtest(model, alternative = 'two.sided')
##
## Durbin-Watson test
##
## data: model
## DW = 0.73266, p-value = 0.001914
## alternative hypothesis: true autocorrelation is not 0
karena p-value < 5% maka Tolak H0 sehingga cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::bgtest(y ~ x, data=Tugas_Individu_MPDW, order=1)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: y ~ x
## LM test = 4.5615, df = 1, p-value = 0.0327
karena p-value < 5% maka Tolak H0 sehingga cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
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 = -1.2111, p-value = 0.2259
karena p-value > 5% maka Terima H0 sehingga tidak cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
#Interactive method using to solve first order autocorrelation problems.
modelco <- orcutt::cochrane.orcutt(model,convergence = 6,max.iter = 10000)
modelco
## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = y ~ x, data = Tugas_Individu_MPDW)
##
## number of interaction: 8355
## rho 0.967478
##
## Durbin-Watson statistic
## (original): 0.73266 , p-value: 9.57e-04
## (transformed): 1.61531 , p-value: 1.397e-01
##
## coefficients:
## (Intercept) x
## 1388.336118 -0.637736
#rho optimum
rho <- modelco$rho
y
## [1] 64.44 65.12 65.79 66.16 66.75 67.46 68.24 68.86 69.39 70.02 70.01 70.24
y[-1]
## [1] 65.12 65.79 66.16 66.75 67.46 68.24 68.86 69.39 70.02 70.01 70.24
y[-12]
## [1] 64.44 65.12 65.79 66.16 66.75 67.46 68.24 68.86 69.39 70.02 70.01
#transformasi terhadap y dan x
(y.trans <- y[-1]-y[-12]*rho)
## [1] 2.775750 2.787865 2.509655 2.741688 2.880877 2.973968 2.839335 2.769499
## [9] 2.886736 2.267225 2.506900
(x.trans <- x[-1]-x[-12]*rho)
## [1] 66.37022 66.40274 66.43526 66.46779 66.50031 66.53283 66.56535 66.59788
## [9] 66.63040 66.66292 66.69544
#model baru
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho)
##
## Call:
## lm(formula = y.trans ~ x.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37158 -0.08045 -0.01687 0.13833 0.25220
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.1522 40.6228 1.111 0.295
## x.trans -0.6377 0.6106 -1.044 0.323
##
## Residual standard error: 0.2083 on 9 degrees of freedom
## Multiple R-squared: 0.1081, Adjusted R-squared: 0.009016
## F-statistic: 1.091 on 1 and 9 DF, p-value: 0.3235
H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::dwtest(modelcorho,alternative = 'two.sided')
##
## Durbin-Watson test
##
## data: modelcorho
## DW = 1.6153, p-value = 0.2793
## alternative hypothesis: true autocorrelation is not 0
karena p-value > 5% maka Tak Tolak H0 sehingga tidak cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::bgtest(modelcorho, order=1)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelcorho
## LM test = 0.34608, df = 1, p-value = 0.5563
karena p-value > 5% maka Tak tolak H0 sehingga tidak cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
H0: tidak ada autokorelasi
H1: ada autokorelasi
lawstat::runs.test(resid(modelcorho), alternative = 'two.sided')
##
## Runs Test - Two sided
##
## data: resid(modelcorho)
## Standardized Runs Statistic = -2.2162, p-value = 0.02668
karena p-value < 5% maka Tolak H0 sehingga cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
# 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))}))
tab$rho[which.min(tab$SSE)]
## [1] 0.99
#grafik rho dan SSE
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)
#Model terbaik
modelhl <- hildreth.lu.func(0.99, model)
summary(modelhl)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37552 -0.07479 -0.01715 0.13791 0.25256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.277 41.973 1.746 0.115
## x -3.408 1.985 -1.717 0.120
##
## Residual standard error: 0.2081 on 9 degrees of freedom
## Multiple R-squared: 0.2468, Adjusted R-squared: 0.1631
## F-statistic: 2.949 on 1 and 9 DF, p-value: 0.1201
H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::dwtest(modelhl,alternative = "two.sided")
##
## Durbin-Watson test
##
## data: modelhl
## DW = 1.6491, p-value = 0.308
## alternative hypothesis: true autocorrelation is not 0
karena p-value > 5% maka Tak tolak H0 sehingga tidak cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
H0: tidak ada autokorelasi
H1: ada autokorelasi
bgtest(modelhl)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelhl
## LM test = 0.28982, df = 1, p-value = 0.5903
karena p-value > 5% maka Tak tolak H0 sehingga tidak cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
H0: tidak ada autokorelasi
H1: ada autokorelasi
lawstat::runs.test(resid(modelhl), alternative = 'two.sided')
##
## Runs Test - Two sided
##
## data: resid(modelhl)
## Standardized Runs Statistic = -2.2162, p-value = 0.02668
karena p-value < 5% maka Tolak H0 sehingga cukup bukti untuk menyatakan bahwa terdapat autokorelasi antar amatan pada taraf 5%
cat("IPM = ", coef(modelhl)[1]/(1-0.735), "+", coef(modelhl)[2]," Tahun", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal
## IPM = 276.5185+-3.407818 Tahun
#Kesimpulan
Terdeteksi adanya autokorelasi pada data IPM Provinsi Sumatera Selatan sehingga dilakukan transpormasi cochrane-orcut dan Hildreth-Lu, setelah dilakukan transpormasi cochrane-orcut dan Hildreth-Lu hasil uji formal pada Runs Test tetap terdeteksi adanya autokorelasi, berbeda dengan dua uji lainnya (Durbin-Watson dan Breusch-Godfrey) yang tidak terdeteksi autokorelasi.