Tugas Praktikum 2 - MPDW
Sandra Berliana Putri
2024-09-01
Penanganan Autokorelasi
Pemanggilan Packages
library (dplyr)## Warning: package 'dplyr' was built under R version 4.3.2##
## 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.3.2library (forecast)## Warning: package 'forecast' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary (lmtest) #digunakan untuk uji formal pendeteksian autokorelasi## Warning: package 'lmtest' was built under R version 4.3.2## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary (orcutt) #untuk membuat model regresi Cochrane-Orcutt## Warning: package 'orcutt' was built under R version 4.3.3library (HoRM) #untuk membuat model regresi Hildreth-Lu## Warning: package 'HoRM' was built under R version 4.3.3Import Data
Data yang digunakan dalam kesempatan kali ini adalah data PDB (GDP) Negara Tunisia periode tahun 1961-2023.
library(readxl)## Warning: package 'readxl' was built under R version 4.3.2dtTunisia <- read_excel("D:/Sem 5/MPDW/Tugas/Individu/2/Tunisia.xlsx")
dtTunisia## # A tibble: 63 × 2
## Tahun GDP
## <dbl> <dbl>
## 1 1961 866155429.
## 2 1962 880027733.
## 3 1963 1026737600
## 4 1964 1025866792.
## 5 1965 991047619.
## 6 1966 1040952381.
## 7 1967 1085714286.
## 8 1968 1214666667.
## 9 1969 1289904762.
## 10 1970 1439238095.
## # ℹ 53 more rowsEksplorasi Data 1
Sebelum melakukan regresi, akan diperlihatkan plot time-series dari PDB (GDP) Negara Tunisia Periode 1961-2023.
#Membentuk objek time series
data.ts<-ts(dtTunisia$GDP)
data.ts## Time Series:
## Start = 1
## End = 63
## Frequency = 1
## [1] 866155429 880027733 1026737600 1025866792 991047619 1040952381
## [7] 1085714286 1214666667 1289904762 1439238095 1685162272 2237556149
## [13] 2730813385 3545868575 4328965588 4508191942 5109324009 5968460080
## [19] 7188863904 8744134354 8428445294 8133580052 8350582748 8254541195
## [25] 8410226053 9017806654 9696715911 10096245762 10101851745 12290568182
## [31] 13074782609 15496708060 14608335608 15633174304 18030876599 19587161807
## [37] 20746210354 21802893587 22943202175 21473528161 22065832449 23141616605
## [43] 27453902261 31183885241 32272186695 34376664601 38915353867 44859439902
## [49] 43455740497 46206091938 48123325825 47311401813 48685446414 50271812921
## [55] 45779494042 44360072680 42163530591 42686580021 41905540184 42491816476
## [61] 46812290143 44579761023 48529595417#Membuat plot time series
ts.plot(data.ts, xlab="Time Period ", ylab="GDP", main= "Time Series Plot of GDP")
points(data.ts)
Pemulusan
Selanjutnya akan dilakukan ramalan dan pemulusan dengan metode DMA dan DES karena terlihat pada plot di atas menunjukkan adanya trend.
DMA
dt.sma <- SMA(data.ts, n=3)
dma <- SMA(dt.sma, n = 3)
At <- 2*dt.sma - dma
Bt <- 2/(3-1)*(dt.sma - dma)
dt.dma<- At+Bt
dt.ramal<- c(NA, dt.dma)
t = 1:5
f = c()
for (i in t) {
f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}dt.gab <- cbind(aktual = c(data.ts,rep(NA,5)),
pemulusan1 = c(dt.sma,rep(NA,5)),
pemulusan2 = c(dt.dma, rep(NA,5)),
At = c(At, rep(NA,5)),
Bt = c(Bt,rep(NA,5)),
ramalan = c(dt.ramal, f[-1]))
dt.gab## aktual pemulusan1 pemulusan2 At Bt ramalan
## [1,] 866155429 NA NA NA NA NA
## [2,] 880027733 NA NA NA NA NA
## [3,] 1026737600 924306921 NA NA NA NA
## [4,] 1025866792 977544042 NA NA NA NA
## [5,] 991047619 1014550670 1099384256 1056967463 42416793 NA
## [6,] 1040952381 1019288931 1050277697 1034783314 15494383 1099384256
## [7,] 1085714286 1039238095 1068995821 1054116958 14878863 1050277697
## [8,] 1214666667 1113777778 1226463464 1170120621 56342843 1068995821
## [9,] 1289904762 1196761905 1357100529 1276931217 80169312 1226463464
## [10,] 1439238095 1314603175 1527047619 1420825397 106222222 1357100529
## [11,] 1685162272 1471435043 1759105048 1615270046 143835002 1527047619
## [12,] 2237556149 1787318839 2313051812 2050185325 262866487 1759105048
## [13,] 2730813385 2217843935 3002466595 2610155265 392311330 2313051812
## [14,] 3545868575 2838079370 3952076680 3395078025 556998655 3002466595
## [15,] 4328965588 3535215849 4878221445 4206718647 671502798 3952076680
## [16,] 4508191942 4127675368 5382379047 4755027207 627351839 4878221445
## [17,] 5109324009 4648827180 5738669274 5193748227 544921047 5382379047
## [18,] 5968460080 5195325344 6271424104 5733374724 538049380 5738669274
## [19,] 7188863904 6088882664 7644624534 6866753599 777870935 6271424104
## [20,] 8744134354 7300486113 9511662257 8406074185 1105588072 7644624534
## [21,] 8428445294 8120481184 10021543578 9071012381 950531197 9511662257
## [22,] 8133580052 8435386567 9401923791 8918655179 483268612 10021543578
## [23,] 8350582748 8304202698 8339227795 8321715246 17512548 9401923791
## [24,] 8254541195 8246234665 8081488041 8163861353 -82373312 8339227795
## [25,] 8410226053 8338449998 8422758421 8380604210 42154211 8081488041
## [26,] 9017806654 8560857967 8918878815 8739868391 179010424 8422758421
## [27,] 9696715911 9041582873 9830821392 9436202132 394619260 8918878815
## [28,] 10096245762 9603589443 10673414806 10138502124 534912682 9830821392
## [29,] 10101851745 9964937806 10821406670 10393172238 428234432 10673414806
## [30,] 12290568182 10829555230 12223277370 11526416300 696861070 10821406670
## [31,] 13074782609 11822400845 13722606615 12772503730 950102885 12223277370
## [32,] 15496708060 13620686284 16680297279 15150491781 1529805498 13722606615
## [33,] 14608335608 14393275426 16622251241 15507763333 1114487908 16680297279
## [34,] 15633174304 15246072658 16898195061 16072133859 826061202 16622251241
## [35,] 18030876599 16090795504 17785624120 16938209812 847414308 16898195061
## [36,] 19587161807 17750404237 20526364445 19138384341 1387980104 17785624120
## [37,] 20746210354 19454749587 22833615875 21144182731 1689433144 20526364445
## [38,] 21802893587 20712088583 23524770810 22118429697 1406341114 22833615875
## [39,] 22943202175 21830768705 24160568200 22995668452 1164899747 23524770810
## [40,] 21473528161 22073207974 23142247081 22607727528 534519554 24160568200
## [41,] 22065832449 22160854262 22439342157 22300098209 139243948 23142247081
## [42,] 23141616605 22226992405 22373607454 22300299929 73307525 22439342157
## [43,] 27453902261 24220450438 26922486578 25571468508 1351018070 22373607454
## [44,] 31183885241 27259801369 32641241299 29950521334 2690719965 26922486578
## [45,] 32272186695 30303324733 36387589838 33345457285 3042132553 32641241299
## [46,] 34376664601 32610912179 37716711016 35163811598 2552899419 36387589838
## [47,] 38915353867 35188068387 40162668296 37675368342 2487299954 37716711016
## [48,] 44859439902 39383819457 46696258354 43040038905 3656219449 40162668296
## [49,] 43455740497 42410178089 49242490311 45826334200 3416156111 46696258354
## [50,] 46206091938 44840424112 50098324566 47469374339 2628950227 49242490311
## [51,] 48123325825 45928386087 48999166068 47463776078 1535389991 50098324566
## [52,] 47311401813 47213606525 49652541759 48433074142 1219467617 48999166068
## [53,] 48685446414 48040058017 49998806966 49019432491 979374474 49652541759
## [54,] 50271812921 48756220382 50262071197 49509145790 752925407 49998806966
## [55,] 45779494042 48245584459 48042178138 48143881298 -101703161 50262071197
## [56,] 44360072680 46803793214 44540980939 45672387077 -1131406138 48042178138
## [57,] 42163530591 44101032438 39536157239 41818594839 -2282437599 44540980939
## [58,] 42686580021 43070061097 39893592126 41481826612 -1588234486 39536157239
## [59,] 41905540184 42251883599 40473666040 41362774819 -889108779 39893592126
## [60,] 42491816476 42361312227 41961765399 42161538813 -199773414 40473666040
## [61,] 46812290143 43736548934 45643150297 44689849616 953300681 41961765399
## [62,] 44579761023 44627955881 46733322947 45680639414 1052683533 45643150297
## [63,] 48529595417 46640548861 49918277465 48279413163 1638864302 46733322947
## [64,] NA NA NA NA NA 49918277465
## [65,] NA NA NA NA NA 51557141768
## [66,] NA NA NA NA NA 53196006070
## [67,] NA NA NA NA NA 54834870372
## [68,] NA NA NA NA NA 56473734675#Plot time series
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="GDP",
main= "DMA N=3 Data GDP")
points(dt.gab[,1])
points(dt.gab[,3])
points(dt.gab[,6])
lines(dt.gab[,3],col="green",lwd=2)
lines(dt.gab[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"),
lty=8, col=c("black","green","red"), cex=0.8)
Selanjutnya akan dilihat keakuratan dari metode DMA
#Menghitung nilai keakuratan
error.dma = data.ts-dt.ramal[1:length(data.ts)]
SSE.dma = sum(error.dma[6:length(data.ts)]^2)
MSE.dma = mean(error.dma[6:length(data.ts)]^2)
MAPE.dma = mean(abs((error.dma[6:length(data.ts)]/data.ts[6:length(data.ts)])*100))
akurasi.dma <- matrix(c(SSE.dma, MSE.dma, MAPE.dma))
row.names(akurasi.dma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi.dma) <- c("Akurasi m = 3")
akurasi.dma## Akurasi m = 3
## SSE 2.210433e+20
## MSE 3.811091e+18
## MAPE 7.750508e+00DES
Selanjutnya akan digunakan metode Double Exponential Smoothing dengan cara sebagai berikut.
Pertama akan data akan dibagi menjadi data training dan data testing.
#membagi training dan testing
training<-dtTunisia[1:50,2]
testing<-dtTunisia[51:63,2]
#data time series
training.ts<-ts(training)
testing.ts<-ts(testing,start=51)
#eksplorasi data
plot(data.ts, col="red",main="Plot semua data")
points(data.ts)
plot(training.ts, col="blue",main="Plot data training")
points(training.ts)
plot(testing.ts, col="green",main="Plot data testing")
points(testing.ts)
Selanjutnya akan dilakukan pemulusan dengan DES, kali ini langsung dicari lambda dan gamma optimum sebagai berikut. Nilai lambda dan gamma optimum dapat dilihat pada smoothing parameters alpha untuk nilai lambda dan beta untuk nilai gamma.
#Lamda dan gamma optimum
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: 1
## beta : 0.1614014
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 46206091938
## b 2192615677#ramalan
ramalandesopt<- forecast(des.opt, h=5)
ramalandesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 51 48398707615 46714011054 50083404176 45822186566 50975228664
## 52 50591323292 48009362648 53173283937 46642555245 54540091339
## 53 52783938969 49373245871 56194632067 47567734008 58000143931
## 54 54976554646 50747489144 59205620149 48508757098 61444352195
## 55 57169170323 52112402954 62225937693 49435512007 64902828640plot(des.opt)
legend("topleft", c("Data Aktual", "Peramalan"), col = c("black", "red"),
lty = c(1,1))
Selanjutnya akan dicari akurasi dari metode DES.
#Akurasi data training
ssedes.train<-des.opt$SSE
msedes.train<-ssedes.train/length(training.ts)
sisaandes<-ramalandesopt$residuals
head(sisaandes)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## x
## [1,] NA
## [2,] NA
## [3,] 132837562
## [4,] -36183275
## [5,] -64291611
## [6,] 30809077mapedes.train <- sum(abs(sisaandes[3:length(training.ts)]/training.ts[3:length(training.ts)])*100)/length(training.ts)
akurasides.opt <- matrix(c(ssedes.train,msedes.train,mapedes.train))
row.names(akurasides.opt)<- c("SSE", "MSE", "MAPE")
colnames(akurasides.opt) <- c("Akurasi lamda dan gamma optimum")
akurasides.opt## Akurasi lamda dan gamma optimum
## SSE 8.501740e+19
## MSE 1.700348e+18
## MAPE 6.408108e+00#Akurasi data testing
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt## Time Series:
## Start = 51
## End = 55
## Frequency = 1
## testing.ts
## [1,] 275381790
## [2,] 3279921479
## [3,] 4098492555
## [4,] 4704741726
## [5,] 11389676282SSEtestingdesopt<-sum(selisihdesopt^2)
SSEtestingdesopt<-SSEtestingdesopt/length(testing.ts)
MAPEtestingdesopt<-sum(abs(selisihdesopt/testing.ts)*100)/length(testing.ts)
akurasiDesTesting <- matrix(c(SSEtestingdesopt,SSEtestingdesopt,MAPEtestingdesopt))
row.names(akurasiDesTesting)<- c("SSE", "MSE", "MAPE")
colnames(akurasiDesTesting) <- c("Akurasi lamda dan gamma optimum")
akurasiDesTesting## Akurasi lamda dan gamma optimum
## SSE 1.380698e+19
## MSE 1.380698e+19
## MAPE 3.858555e+00Bandingkan DMA dan DES
Setelah didapatkan nilai akurasi untuk metode DMA dan DES, selanjutnya akan dibandingkan keakuratan antar metode keduanya.
cbind(akurasi.dma, akurasides.opt)## Akurasi m = 3 Akurasi lamda dan gamma optimum
## SSE 2.210433e+20 8.501740e+19
## MSE 3.811091e+18 1.700348e+18
## MAPE 7.750508e+00 6.408108e+00Berdasarkan perbandingan akurasi tersebut, terlihat nilai SSE, MSE, dan MAPE metode DES lebih kecil dibandingkan dengan metode DMA. Oleh karena itu, metode peramalan dan pemulusan yang terbaik antara keduanya adalah dengan metode DES.
Eksplorasi Data 2
Setelah melakukan peramalan, data yang telah dimasukkan kemudian dieksplorasi. Eksplorasi pertama yang dilakukan adalah dengan menggunakan scatter plot.
#Eksplorasi Data
#Pembuatan Scatter Plot
plot(dtTunisia$Tahun,dtTunisia$GDP, pch = 20, col = "blue",
main = "Scatter Plot Tahun vs Nilai GDP",
xlab = "Tahun",
ylab = "Nilai GDP")
#Menampilkan Nilai Korelasi
cor(dtTunisia$Tahun,dtTunisia$GDP)## [1] 0.9544216Berdasarkan scatter plot di atas, terlihat adanya hubungan / korelasi
positif antara peubah tahun dengan nilai PDB (GDP), terlihat titik-titik
pada plot yang naik ke arah kanan atas. Hal tersebut juga diperkuat
dengan hasil perhitungan aplikasi R di mana didapatkan
nilai korelasi sebesar \(0.9544216\).
Setalah mengetahui adanya hubungan antar dua peubah, maka model regresi dapat ditentukan.
Regresi
#Pembuatan Model Regresi
#model regresi
model<- lm(dtTunisia$GDP~dtTunisia$Tahun, data = dtTunisia)
summary(model)##
## Call:
## lm(formula = dtTunisia$GDP ~ dtTunisia$Tahun, data = dtTunisia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.354e+09 -4.075e+09 -6.242e+08 3.145e+09 1.076e+10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.782e+12 7.216e+10 -24.70 <2e-16 ***
## dtTunisia$Tahun 9.047e+08 3.622e+07 24.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.228e+09 on 61 degrees of freedom
## Multiple R-squared: 0.9109, Adjusted R-squared: 0.9095
## F-statistic: 623.8 on 1 and 61 DF, p-value: < 2.2e-16Model yang dihasilkan adalah \[y_i=(-1.782e+12)+(9.047e+08)x_t\]. Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value < \(\alpha\) (5%). Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model. Hasil uji-t parsial kedua parameter regresi, yaitu intersep dan koefisien regresi juga menunjukkan hal yang sama, yaitu memiliki p-value < \(\alpha\) (5%) sehingga nyata dalam taraf 5%. Selanjutnya dapat dilihat juga nilai \(R^2=0.9109\). Artinya, sebesar 91.09% keragaman nilai PDB (GDP) dapat dijelaskan oleh peubah tahun. Hasil ini menunjukkan hasil yang bagus, seolah mendapatkan hasil terbaik.
Uji Asumsi Regresi dengan Uji Non Formal
Namun, kita perlu melakukan uji terhadap sisaannya seperti berikut ini.
# Hitung sisaan dan fitted values dari model
sisaan <- residuals(model)
fitValue <- predict(model)
# Panjang sisaan (jumlah observasi)
n <- length(sisaan)
# Diagnostik dengan eksploratif
par(mfrow = c(2,2))
# QQ plot untuk memeriksa normalitas residuals
qqnorm(sisaan)
qqline(sisaan, col = "steelblue", lwd = 2)
# Plot Sisaan vs Fitted Values
plot(fitValue, sisaan, col = "steelblue", pch = 20,
xlab = "Fitted Values", ylab = "Sisaan",
main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
# Histogram sisaan
hist(sisaan, col = "steelblue", main = "Histogram Sisaan", xlab = "Sisaan")
# Plot Sisaan vs Order (berdasarkan jumlah observasi)
plot(seq(1, n, 1), sisaan, col = "steelblue", pch = 20,
xlab = "Order", ylab = "Sisaan",
main = "Sisaan vs Order")
lines(seq(1, n, 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 dan histogram di atas menunjukkan bahwa sisaan cenderung menyebar normal. 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.
#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)
Berdasarkan plot ACF, terlihat lag 1 hingga 6 berada di luar rentang
batas dan signifikan. Sedangkan berdasarkan plot PACF, terlihat lag 1
berada di luar rentang batas dan signifikan. Namun, untuk lebih
memastikan akan dilakukan uji formal dengan uji Durbin Watson.
Uji Asumsi Regresi dengan Uji Formal
#Melihat Sisaan Menyebar Normal/Tidak
#H0: sisaan mengikuti sebaran normal
#H1: sisaan tidak mengikuti sebaran normal
ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan)) #menggunakan ks test karena sampel > 50##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.1272, p-value = 0.2388
## alternative hypothesis: two-sidedBerdasarkan uji formal Kolmogorov-Smirnov didapatkan nilai p-value > \(\alpha\) (5%). Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal.
#Deteksi autokorelasi dengan uji-Durbin Watson
#H0: tidak ada autokorelasi
#H1: ada autokorelasi
dwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 0.10873, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0Berdasarkan hasil DW Test, didapatkan nilai \(DW = 0.10873\) dan p-value = \(< 2.2e-16\). Berdasarkan tabel Durbin-Watson diperoleh nilai \(DL = 1.5599\) dan \(DU = 1.6243\). Nilai DW berada di antara nilai 0 dan DL. Artinya, dapat dikatakan berada di daerah autokorelasi positif. Pernyataan tersebut didukung juga oleh nilai p-value < 0.05 dapat disimpulkan bahwa tolak H0, cukup bukti mengatakan adanya autokorelasi. Oleh karena itu, diperlukan penangan autokorelasi. Penanganan yang akan digunakan menggunakan dua metode, yaitu Cochrane-Orcutt dan Hildret-Lu.
Penanganan Autokorelasi
Sisaan saling bebas atau tidak ada autokorelasi adalah keadaan ketika error yang satu bukan merupakan error yang lainnya. Ketika error yang satu merupakan fungsi error dari yang lain maka errornya tidak akan pernah menjadi minimum.
Metode Cochrane-Orcutt
Penanganan metode Cochrane-Orcutt dapat dilakukan dengan bantuan
packages Orcutt pada aplikasi R maupun secara manual.
Berikut ini ditampilkan cara menggunakan bantuan library
packages Orcutt.
#Penanganan Autokorelasi Cochrane-Orcutt
modelCO<-cochrane.orcutt(model)
modelCO## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = dtTunisia$GDP ~ dtTunisia$Tahun, data = dtTunisia)
##
## number of interaction: 15
## rho 0.922803
##
## Durbin-Watson statistic
## (original): 0.10873 , p-value: 1.407e-34
## (transformed): 1.71846 , p-value: 1.051e-01
##
## coefficients:
## (Intercept) dtTunisia$Tahun
## -2.155802e+12 1.090309e+09Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=-2.156766000000+1090784250x_t\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(1.71846\) dan \(0.1051\). Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.6243 < DW < 2.3757\). Hal tersebut juga didukung dengan nilai p-value > 0.05, artinya belum cukup bukti menyatakan bahwa sisaan terdapat autokorelasi pada taraf nyata 5%. Untuk nilai \(ρ ̂\) optimum yang digunakan adalah \(0.922803\). Nilai tersebut dapat diketahui dengan syntax berikut.
#Rho optimum
rho<- modelCO$rho
rho## [1] 0.9228032Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.
dtTunisia$GDP #yt## [1] 866155429 880027733 1026737600 1025866792 991047619 1040952381
## [7] 1085714286 1214666667 1289904762 1439238095 1685162272 2237556149
## [13] 2730813385 3545868575 4328965588 4508191942 5109324009 5968460080
## [19] 7188863904 8744134354 8428445294 8133580052 8350582748 8254541195
## [25] 8410226053 9017806654 9696715911 10096245762 10101851745 12290568182
## [31] 13074782609 15496708060 14608335608 15633174304 18030876599 19587161807
## [37] 20746210354 21802893587 22943202175 21473528161 22065832449 23141616605
## [43] 27453902261 31183885241 32272186695 34376664601 38915353867 44859439902
## [49] 43455740497 46206091938 48123325825 47311401813 48685446414 50271812921
## [55] 45779494042 44360072680 42163530591 42686580021 41905540184 42491816476
## [61] 46812290143 44579761023 48529595417dtTunisia$GDP[-1] #y(t-1)## [1] 880027733 1026737600 1025866792 991047619 1040952381 1085714286
## [7] 1214666667 1289904762 1439238095 1685162272 2237556149 2730813385
## [13] 3545868575 4328965588 4508191942 5109324009 5968460080 7188863904
## [19] 8744134354 8428445294 8133580052 8350582748 8254541195 8410226053
## [25] 9017806654 9696715911 10096245762 10101851745 12290568182 13074782609
## [31] 15496708060 14608335608 15633174304 18030876599 19587161807 20746210354
## [37] 21802893587 22943202175 21473528161 22065832449 23141616605 27453902261
## [43] 31183885241 32272186695 34376664601 38915353867 44859439902 43455740497
## [49] 46206091938 48123325825 47311401813 48685446414 50271812921 45779494042
## [55] 44360072680 42163530591 42686580021 41905540184 42491816476 46812290143
## [61] 44579761023 48529595417#Transformasi Manual
GDP.trans<- dtTunisia$GDP[-1]-dtTunisia$GDP[-63]*rho
Tahun.trans<- dtTunisia$Tahun[-1]-dtTunisia$Tahun[-63]*rho
modelCOmanual<- lm(GDP.trans~Tahun.trans)
summary(modelCOmanual)##
## Call:
## lm(formula = GDP.trans ~ Tahun.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.796e+09 -7.617e+08 -1.737e+07 4.969e+08 5.353e+09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.664e+11 2.367e+10 -7.031 2.22e-09 ***
## Tahun.trans 1.090e+09 1.530e+08 7.128 1.51e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.664e+09 on 60 degrees of freedom
## Multiple R-squared: 0.4585, Adjusted R-squared: 0.4495
## F-statistic: 50.81 on 1 and 60 DF, p-value: 1.514e-09Hasil model transformasi bukan merupakan model sesungguhnya. Koefisien regresi masih perlu dicari kembali mengikuti \(β_0^*=β_0+ρ ̂β_0\) dan \(β_1^*=β_1\).
#Mencari Penduga Koefisien Regresi setelah Transformasi ke Persamaan Awal
b0bintang <- modelCOmanual$coefficients[-2]
b0 <- b0bintang/(1-rho)
b1 <- modelCOmanual$coefficients[-1]
b0## (Intercept)
## -2.155802e+12b1## Tahun.trans
## 1090308580Hasil perhitungan koefisien regresi tersebut akan menghasilkan hasil yang sama dengan model yang dihasilkan menggunakan packages.
Metode Hildreth-Lu
Penanganan kedua adalah menggunakan 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 #karena jika mulai dari satu maka t-1 = 0 dan ini tidak ada
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))}))
round(tab, 4)## rho SSE
## 1 0.1 1.294893e+21
## 2 0.2 1.037193e+21
## 3 0.3 8.128388e+20
## 4 0.4 6.218309e+20
## 5 0.5 4.641692e+20
## 6 0.6 3.398536e+20
## 7 0.7 2.488842e+20
## 8 0.8 1.912610e+20
## 9 0.9 1.669839e+20Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.9. 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.8 sampai dengan 1.
#Rho optimal di sekitar 0.8 hingga 1
rOpt <- seq(0.8,1, by= 0.001)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model))}))
head(tabOpt[order(tabOpt$SSE),])## rho SSE
## 124 0.923 1.661170e+20
## 123 0.922 1.661180e+20
## 125 0.924 1.661193e+20
## 122 0.921 1.661224e+20
## 126 0.925 1.661250e+20
## 121 0.920 1.661300e+20#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.923, y=1.661170e+20 , labels = "rho=0.923", cex = 0.8)
Perhitungan yang dilakukan aplikasi R menunjukkan bahwa
nilai \(ρ\) optimum, yaitu saat SSE
terkecil terdapat pada nilai \(ρ=0.923\). Hal tersebut juga ditunjukkan
pada plot. Selanjutnya, model dapat didapatkan dengan mengevaluasi nilai
\(ρ\) ke dalam fungsi
hildreth.lu.func, serta dilanjutkan dengan pengujian
autokorelasi dengan uji Durbin-Watson. Namun, setelah pengecekan
tersebut tidak lupa koefisien regresi tersebut digunakan untuk
transformasi balik. Persamaan hasil transformasi itulah yang menjadi
persamaan sesungguhnya.
#Model terbaik
modelHL <- hildreth.lu.func(0.923, model)
summary(modelHL)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.798e+09 -7.606e+08 -1.815e+07 4.968e+08 5.352e+09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.661e+11 2.367e+10 -7.016 2.35e-09 ***
## x 1.091e+09 1.534e+08 7.113 1.61e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.664e+09 on 60 degrees of freedom
## Multiple R-squared: 0.4575, Adjusted R-squared: 0.4484
## F-statistic: 50.59 on 1 and 60 DF, p-value: 1.606e-09#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.923), "+", coef(modelHL)[2],"x", sep = "")## y = -2.156766e+12+1090784250xSetelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=-2.156766000000+1090784250x_t\]
#Deteksi autokorelasi
dwtest(modelHL)##
## Durbin-Watson test
##
## data: modelHL
## DW = 1.7188, p-value = 0.1053
## alternative hypothesis: true autocorrelation is greater than 0Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(1.7188\) berada pada selang daerah tidak ada autokorelasi, yaitu pada rentang DU < DW < 4-DU atau \(1.6243 < DW < 2.3757\). Hal tersebut juga didukung oleh p-value sebesar \(0.1053\), di mana p-value > \(\alpha\)=5%. Artinya tak tolak \(H_0\) atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai IPM dengan metode Hildreth-Lu pada taraf nyata 5%.
Terakhir, akan dibandingkan nilai SSE dari ketiga metode (metode awal, metode Cochrane-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(dtTunisia$GDP)
mseModelCO <- sseModelCO/length(dtTunisia$GDP)
mseModelHL <- sseModelHL/length(dtTunisia$GDP)
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 1.667413e+21 1.661169e+20 1.661170e+20
## MSE 2.646687e+19 2.636777e+18 2.636778e+18Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki SSE yang sama, sebesar \(1.661169e+20\) dan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar \(1.667413e+21\).
Simpulan
Autokorelasi yang terdapat pada data PDB (GDP) terjadi akibat adanya korelasi di antara unsur penyusunnya. Indikator PDB (GDP) yang erat hubungannya dengan perekonomian sangat rawan menjadi penyebab adanya autokorelasi. Adanya autokorelasi menyebabkan model regresi kurang baik karena akan meingkatkan galatnya. Autokorelasi dapat dideteksi secara eksploratif melalui plot sisaan, ACF, dan PACF, serta dengan uji formal Durbin-Watson. Namun, autokorelasi tersebut dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode menghasilkan nilai SSE yang sama, artinya keduanya baik untuk digunakan.
Daftar Pustaka
Aprianto A, Debataraja NN, Imro’ah N. 2020. Metode cochrane-orcutt untuk mengatasi autokorelasi pada estimasi parameter ordinary least squares. Bimaster : Buletin Ilmiah Matematika, Statistika dan Terapannya. 9(1):95–102. doi:10.26418/bbimst.v9i1.38590.