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, unionlibrary(readxl)
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(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.numericlibrary(orcutt) #Cochrane-Orcutt
library(HoRM)#Hildreth Lu#File Data
dataregresi <- read_excel("~/Desktop/SEMESTER 5!!!/MPDW/Data Tugas Indivisu 1.xlsx", sheet = 3)
dataregresi## # A tibble: 12 × 2
## Tahun IPM
## <dbl> <dbl>
## 1 2010 67.1
## 2 2011 67.3
## 3 2012 68.4
## 4 2013 68.9
## 5 2014 69.4
## 6 2015 70.0
## 7 2016 70.7
## 8 2017 71.2
## 9 2018 71.7
## 10 2019 72.4
## 11 2020 72.4
## 12 2021 72.6#Eksplorasi Data
x <- dataregresi$Tahun
y <- dataregresi$IPM
#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")
#Korelasi x dan y
cor(x,y)## [1] 0.9908222#Model Regresi
model <- lm(y~x, data = dataregresi)
summary(model)##
## Call:
## lm(formula = y ~ x, data = dataregresi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51500 -0.13932 0.07545 0.20636 0.31045
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.024e+03 4.719e+01 -21.69 9.69e-10 ***
## x 5.427e-01 2.341e-02 23.18 5.05e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.28 on 10 degrees of freedom
## Multiple R-squared: 0.9817, Adjusted R-squared: 0.9799
## F-statistic: 537.3 on 1 and 10 DF, p-value: 5.05e-10#Deteksi Autokorelasi (Grafik)
#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 = "steelblue")
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)
#Deteksi Autokorelasi (Grafik)
#ACF dan PACF Plot
par(mfrow = c(2,2))
acf(resi1)
pacf(resi1)
Kesimpulan : Berdasarkan 2 pengujian autokorelasi melalui grafik, diperoleh bahwa sisaan pada grafik QQ Plot sudah menyebar normal, histogram yang menunjukan penyebaran sisaan pada data juga terlihat menyebar normal, serta grafik Sisaan vs Fitted Values dan Sisaan vs Order tidak membentuk pola tertentu, sehingga dapat diasumsikan bahwa kebebasan sisaan terpenuhi. Pada pengujian grafik ACF dan PACF terlihat tideak ada yang signifikan pada keduanya maka disimpulkan bahwa sisaan juga saling bebas.
#Deteksi Autokorelasi (Statistik Uji)
#Durbin-Watson
lmtest::dwtest(model, alternative = 'two.sided')##
## Durbin-Watson test
##
## data: model
## DW = 0.98038, p-value = 0.01527
## alternative hypothesis: true autocorrelation is not 0#Deteksi Autokorelasi (Statistik Uji)
#Breusch-Godfrey Test
lmtest::bgtest(y ~ x, data=dataregresi, order=1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: y ~ x
## LM test = 2.5034, df = 1, p-value = 0.1136#Deteksi Autokorelasi (Statistik Uji)
#Run's Test
lawstat::runs.test(resid(model), alternative = 'two.sided')##
## Runs Test - Two sided
##
## data: resid(model)
## Standardized Runs Statistic = -1.2111, p-value = 0.2259Kesimpulan : Berdasarkan 3 Uji Statistik yang sudah dilakukan, diperoleh 2 diantaranya memperoleh nilai P-Value diatas 0.05 yang artinya Tak Tolak H0 sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%. Namun pada pengujian Durbin-Watson diperoleh nilai P-Value sebesar 0.01527 atau lebih kecil dari 0.05. Dimana dapat diartikan sebagai Tolak H0 sehingga cukup bukti untuk menyatakan terdapat autokorelasi pada taraf 5%. Oleh karena itu, dilakukan penanganan autokorelasi untuk mengecek data lebih lanjut.
#Penanganan Autokorelasi
#Cochrane-Orcutt
modelco <- orcutt::cochrane.orcutt(model)
modelco ## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = y ~ x, data = dataregresi)
##
## number of interaction: 55
## rho 0.672673
##
## Durbin-Watson statistic
## (original): 0.98038 , p-value: 7.637e-03
## (transformed): 1.83910 , p-value: 2.496e-01
##
## coefficients:
## (Intercept) x
## -871.945087 0.467502#RHO Optimum
rho <- modelco$rho
y## [1] 67.09 67.34 68.36 68.91 69.36 69.98 70.73 71.24 71.73 72.39 72.38 72.65y[-1]## [1] 67.34 68.36 68.91 69.36 69.98 70.73 71.24 71.73 72.39 72.38 72.65y[-12]## [1] 67.09 67.34 68.36 68.91 69.36 69.98 70.73 71.24 71.73 72.39 72.38#Transformasi terhadap y dan x
(y.trans <- y[-1]-y[-12]*rho)## [1] 22.21037 23.06220 22.92607 23.00610 23.32340 23.65634 23.66184 23.80877
## [9] 24.13917 23.68520 23.96193(x.trans <- x[-1]-x[-12]*rho)## [1] 658.9273 659.2546 659.5819 659.9092 660.2366 660.5639 660.8912 661.2185
## [9] 661.5459 661.8732 662.2005#Model Baru
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho)##
## Call:
## lm(formula = y.trans ~ x.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42827 -0.14929 0.07266 0.17881 0.27632
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -285.41116 49.00486 -5.824 0.000252 ***
## x.trans 0.46750 0.07419 6.302 0.000141 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2547 on 9 degrees of freedom
## Multiple R-squared: 0.8152, Adjusted R-squared: 0.7947
## F-statistic: 39.71 on 1 and 9 DF, p-value: 0.0001407#Cek Autokorelasi
lmtest::dwtest(modelcorho,alternative = 'two.sided') ##
## Durbin-Watson test
##
## data: modelcorho
## DW = 1.8391, p-value = 0.4992
## alternative hypothesis: true autocorrelation is not 0Kesimpulan : Dengan menggunakan penanganan Cochrane-Orcutt terlihat bahwa terjadi kenaikan yang signifikan pada nilai P-Value, dimana sebelumnya 0.01527 menjadi 0.4992 atau lebih besar dari 0.05,Tak Tolak H0. Oleh karena itu, disimpulkan metode Cochrane-Orcutt berhasil menangani autokorelasi yang terjadi pada data.
#Penanganan Autokorelasi
#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,7 plg kecil## rho SSE
## 1 0.10 0.7175
## 2 0.20 0.6749
## 3 0.30 0.6404
## 4 0.40 0.6141
## 5 0.50 0.5959
## 6 0.60 0.5859
## 7 0.70 0.5841
## 8 0.80 0.5904
## 9 0.90 0.6049
## 10 0.91 0.6067
## 11 0.92 0.6087
## 12 0.93 0.6108
## 13 0.94 0.6129
## 14 0.95 0.6151
## 15 0.96 0.6174
## 16 0.97 0.6198
## 17 0.98 0.6223
## 18 0.99 0.6248#Grafik rho dan SSE
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)
#Rho optimal di sekitar 0.7
r <- seq(0.6, 0.8, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0.66 kecil## rho SSE
## 1 0.60 0.5859
## 2 0.61 0.5854
## 3 0.62 0.5849
## 4 0.63 0.5845
## 5 0.64 0.5842
## 6 0.65 0.5840
## 7 0.66 0.5838
## 8 0.67 0.5838
## 9 0.68 0.5838
## 10 0.69 0.5839
## 11 0.70 0.5841
## 12 0.71 0.5843
## 13 0.72 0.5847
## 14 0.73 0.5851
## 15 0.74 0.5856
## 16 0.75 0.5862
## 17 0.76 0.5869
## 18 0.77 0.5876
## 19 0.78 0.5885
## 20 0.79 0.5894
## 21 0.80 0.5904#Grafik SSE Optimum
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.66, model)
summary(modelhl)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42841 -0.15140 0.07235 0.18026 0.27711
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -299.10494 49.00485 -6.104 0.000178 ***
## x 0.47135 0.07143 6.599 9.94e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2547 on 9 degrees of freedom
## Multiple R-squared: 0.8287, Adjusted R-squared: 0.8097
## F-statistic: 43.55 on 1 and 9 DF, p-value: 9.936e-05# Deteksi Autokorelasi
lmtest::dwtest(modelhl,alternative = 'two.sided')##
## Durbin-Watson test
##
## data: modelhl
## DW = 1.8161, p-value = 0.4734
## alternative hypothesis: true autocorrelation is not 0# Transformasi Balik
cat("y = ", coef(modelhl)[1]/(1-0.66), "+", coef(modelhl)[2],"x", sep = "")## y = -879.7204+0.4713529xKesimpulan : Dengan menggunakan penanganan Hildreth-lu terlihat bahwa terjadi kenaikan yang signifikan pada nilai P-Value, dimana sebelumnya 0.01527 menjadi 0.4734 atau lebih besar dari 0.05, Tak Tolak H0. Oleh karena itu, disimpulkan metode Hildreth-lu berhasil menangani autokorelasi yang terjadi pada data.