Tugas MPDW
Farid Yafi Suwandi
2024-09-02
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.3library(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) ## Warning: package 'lmtest' was built under R version 4.3.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(orcutt) ## Warning: package 'orcutt' was built under R version 4.3.3library(HoRM) ## Warning: package 'HoRM' was built under R version 4.3.3Input Data
Data yang digunakan dalam kesempatan kali ini adalah data GDP Sri Langka periode tahun 1960-2023.
library(readxl)
DataSL <- read_xlsx("D:/STATATISTIKA/SEM5/MPDW/Data Sri Langka.xlsx")
DataSL## # A tibble: 64 × 2
## Tahun GDP
## <chr> <dbl>
## 1 1960 1409873950.
## 2 1961 1444327731.
## 3 1962 1434156379.
## 4 1963 1240672269.
## 5 1964 1309747899.
## 6 1965 1698319328.
## 7 1966 1751470588.
## 8 1967 1859465021.
## 9 1968 1801344538.
## 10 1969 1965546218.
## # ℹ 54 more rowsEksplorasi Data
Sebelum melakukan regresi, akan diperlihatkan plot time-series dari GDP Sri Langka
#Membentuk objek time series
data.sl<-ts(DataSL$GDP)
data.sl## Time Series:
## Start = 1
## End = 64
## Frequency = 1
## [1] 1409873950 1444327731 1434156379 1240672269 1309747899 1698319328
## [7] 1751470588 1859465021 1801344538 1965546218 2296470588 2369308600
## [13] 2553936348 2875625000 3574586466 3791298146 3591319857 4104509583
## [19] 2733183857 3364611432 4024621900 4415844156 4768765017 5167913302
## [25] 6043474843 5978460972 6405210564 6682167120 6978371581 6987267684
## [31] 8032551173 9000362582 9703011636 10338679636 11717604209 13029697561
## [37] 13897738375 15091913884 15794972847 15656327860 16330814180 15749753805
## [43] 16536535647 18881765437 20662525941 24405791045 28279802406 32350238760
## [49] 40713826215 42066224093 58636161082 67753284044 70447216891 77000578167
## [55] 82528535714 85140955389 88012282206 94376237797 94493871201 89014978344
## [61] 84304298771 88609332762 74144875098 84356860421#Membuat plot time series
ts.plot(data.sl, xlab="Time Period ", ylab="GDP", main= "Time Series Plot of GDP")
points(data.sl)
Selanjutnya akan dilakukan ramalan dan pemulusan dengan metode DMA dan DES karena terlihat pada plot di atas menunjukkan adanya trend.
dt.sma <- SMA(data.sl, 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.sl,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,] 1409873950 NA NA NA NA NA
## [2,] 1444327731 NA NA NA NA NA
## [3,] 1434156379 1429452686 NA NA NA NA
## [4,] 1240672269 1373052126 NA NA NA NA
## [5,] 1309747899 1328192182 1230778550 1279485366 -48706816 NA
## [6,] 1698319328 1416246499 1503745624 1459996062 43749563 1230778550
## [7,] 1751470588 1586512605 1872236958 1729374781 142862176 1503745624
## [8,] 1859465021 1769751646 2127581104 1948666375 178914729 1872236958
## [9,] 1801344538 1804093382 1972041725 1888067553 83974171 2127581104
## [10,] 1965546218 1875451926 1993491141 1934471533 59019608 1972041725
## [11,] 2296470588 2021120448 2262917507 2142018978 120898530 1993491141
## [12,] 2369308600 2210441802 2559982623 2385212213 174770410 2262917507
## [13,] 2553936348 2406571846 2794292806 2600432326 193860480 2559982623
## [14,] 2875625000 2599623316 2987778639 2793700978 194077661 2794292806
## [15,] 3574586466 3001382605 3665762637 3333572621 332190016 2987778639
## [16,] 3791298146 3413836537 4231614639 3822725588 408889051 3665762637
## [17,] 3591319857 3652401490 4245457381 3948929435 296527946 4231614639
## [18,] 4104509583 3829042529 4223607215 4026324872 197282343 4245457381
## [19,] 2733183857 3476337766 3123825441 3300081603 -176256162 4223607215
## [20,] 3364611432 3400768291 3064872482 3232820386 -167947904 3123825441
## [21,] 4024621900 3374139063 3288253776 3331196419 -42942644 3064872482
## [22,] 4415844156 3935025829 4665122033 4300073931 365048102 3288253776
## [23,] 4768765017 4403077024 5401069795 4902073409 498996385 4665122033
## [24,] 5167913302 4784174158 5604337800 5194255979 410081821 5401069795
## [25,] 6043474843 5326717721 6304173893 5815445807 488728086 5604337800
## [26,] 5978460972 5729949706 6629288061 6179618883 449669177 6304173893
## [27,] 6405210564 6142382126 6961113344 6551747735 409365609 6629288061
## [28,] 6682167120 6355279552 6914097733 6634688642 279409091 6961113344
## [29,] 6978371581 6688583088 7274919421 6981751254 293168166 6914097733
## [30,] 6987267684 6882602128 7363496539 7123049334 240447205 7274919421
## [31,] 8032551173 7332730146 8062246863 7697488505 364758359 7363496539
## [32,] 9000362582 8006727146 9205475158 8606101152 599374006 8062246863
## [33,] 9703011636 8911975130 10568303776 9740139453 828164323 9205475158
## [34,] 10338679636 9680684618 11309129257 10494906937 814222320 10568303776
## [35,] 11717604209 10586431827 12306567764 11446499795 860067969 11309129257
## [36,] 13029697561 11695327135 13777685686 12736506410 1041179275 12306567764
## [37,] 13897738375 12881680048 15202747471 14042213760 1160533712 13777685686
## [38,] 15091913884 14006449940 16297045071 15151747505 1145297565 15202747471
## [39,] 15794972847 14928208369 16907066201 15917637285 989428916 16297045071
## [40,] 15656327860 15514404863 16910505809 16212455336 698050473 16907066201
## [41,] 16330814180 15927371629 16868791646 16398081637 470710009 16910505809
## [42,] 15749753805 15912298615 16167512440 16039905527 127606912 16868791646
## [43,] 16536535647 16205701211 16586855996 16396278603 190577393 16167512440
## [44,] 18881765437 17056018296 18385376141 17720697219 664678922 16586855996
## [45,] 20662525941 18693609009 21443941349 20068775179 1375166170 18385376141
## [46,] 24405791045 21316694141 25905868126 23611281134 2294586992 21443941349
## [47,] 28279802406 24449373131 30375001872 27412187501 2962814371 25905868126
## [48,] 32350238760 28345277404 35628269094 31986773249 3641495845 30375001872
## [49,] 40713826215 33781289127 43626574274 38703931700 4922642573 35628269094
## [50,] 42066224093 38376763023 48128069366 43252416194 4875653172 43626574274
## [51,] 58636161082 47138737130 61885018537 54511877833 7373140703 48128069366
## [52,] 67753284044 56151889740 74010742624 65081316182 8929426442 61885018537
## [53,] 70447216891 65612220672 84234763656 74923492164 9311271492 74010742624
## [54,] 77000578167 71733693034 86202543472 78968118253 7234425219 84234763656
## [55,] 82528535714 76658776924 87306537019 81982656971 5323880047 86202543472
## [56,] 85140955389 81556689757 91370629460 86463659608 4906969852 87306537019
## [57,] 88012282206 85227257770 93386623675 89306940722 4079682953 91370629460
## [58,] 94376237797 89176491797 96889182510 93032837153 3856345356 93386623675
## [59,] 94493871201 92294130401 99083804558 95688967480 3394837079 96889182510
## [60,] 89014978344 92628362447 95152430912 93890396679 1262034232 99083804558
## [61,] 84304298771 89271049438 85017453457 87144251448 -2126797991 95152430912
## [62,] 88609332762 87309536626 82455977536 84882757081 -2426779545 85017453457
## [63,] 74144875098 82352835544 74436225559 78394530551 -3958304992 82455977536
## [64,] 84356860421 82370356094 79089249439 80729802767 -1640553327 74436225559
## [65,] NA NA NA NA NA 79089249439
## [66,] NA NA NA NA NA 77448696112
## [67,] NA NA NA NA NA 75808142785
## [68,] NA NA NA NA NA 74167589458
## [69,] NA NA NA NA NA 72527036131#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="yellow",lwd=2)
lines(dt.gab[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"),
lty=8, col=c("black","yellow","red"), cex=0.8)
Selanjutnya akan dilihat keakuratan dari metode DMA
#Menghitung nilai keakuratan
error.dma = data.sl-dt.ramal[1:length(data.sl)]
SSE.dma = sum(error.dma[6:length(data.sl)]^2)
MSE.dma = mean(error.dma[6:length(data.sl)]^2)
MAPE.dma = mean(abs((error.dma[6:length(data.sl)]/data.sl[6:length(data.sl)])*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 7.148121e+20
## MSE 1.211546e+19
## MAPE 8.460205e+00Selanjutnya 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<-DataSL[1:50,2]
testing<-DataSL[51:70,2]
#data time series
training.ts<-ts(training)
testing.ts<-ts(testing,start=51)
#eksplorasi data
plot(data.sl, col="red",main="Plot semua data")
points(data.sl)
plot(training.ts, col="blue",main="Plot data training")
points(training.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: 0.5819323
## beta : 1
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 43627991175
## b 4641090788plot(des.opt)
legend("topleft", c("Data Aktual", "Peramalan"), col = c("black", "red"),
lty = c(1,1))
#ramalan
ramalandesopt<- forecast(des.opt, h=5)
ramalandesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 51 48269081963 46900703861 49637460065 46176328250 50361835676
## 52 52910172752 50810445388 55009900115 49698916864 56121428639
## 53 57551263540 54370737351 60731789728 52687068490 62415458590
## 54 62192354328 57691095358 66693613299 55308272794 69076435862
## 55 66833445116 60823968464 72842921769 57642743687 76024146546Selanjutnya akan dicari akurasi dari metode DES.
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,] -44625134
## [4,] -220625413
## [5,] 96743342
## [6,] 492622698mapedes.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 5.488986e+19
## MSE 1.097797e+18
## MAPE 7.834821e+00#Akurasi data testing
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt## Time Series:
## Start = 51
## End = 55
## Frequency = 1
## testing.ts
## [1,] -10367079118
## [2,] -14843111292
## [3,] -12895953351
## [4,] -14808223839
## [5,] -15695090597SSEtestingdesopt<-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 4.798596e+19
## MSE 4.798596e+19
## MAPE 4.807143e+00Setelah 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 7.148121e+20 5.488986e+19
## MSE 1.211546e+19 1.097797e+18
## MAPE 8.460205e+00 7.834821e+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.
Setelah melakukan peramalan, data yang telah dimasukkan kemudian dieksplorasi. Eksplorasi pertama yang dilakukan adalah dengan menggunakan scatter plot.
#Eksplorasi Data
#Pembuatan Scatter Plot
DataSL$Tahun <- as.numeric(DataSL$Tahun)
DataSL$GDP <- as.numeric(DataSL$GDP)
tahun <- DataSL$Tahun
GDP <- DataSL$GDP
plot(tahun,GDP, pch = 20, col = "red",
main = "Scatter Plot Tahun vs Nilai GDP",
xlab = "Tahun",
ylab = "Nilai GDP")
#Menampilkan Nilai Korelasi
cor(tahun,GDP)## [1] 0.8537009Berdasarkan scatter plot di atas, terlihat adanya hubungan / korelasi positif antara peubah tahun dengan nilai GDP, terlihat titik-titik pada plot yang naik ke arah kanan atas.
Setalah mengetahui adanya hubungan antar dua peubah, maka model regresi dapat ditentukan.
Regresi
#Pembuatan Model Regresi
#model regresi
model<- lm(GDP~tahun, data = DataSL)
summary(model)##
## Call:
## lm(formula = GDP ~ tahun, data = DataSL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.396e+10 -1.543e+10 -1.710e+09 1.435e+10 3.223e+10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.849e+12 2.227e+11 -12.79 <2e-16 ***
## tahun 1.444e+09 1.118e+08 12.91 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.653e+10 on 62 degrees of freedom
## Multiple R-squared: 0.7288, Adjusted R-squared: 0.7244
## F-statistic: 166.6 on 1 and 62 DF, p-value: < 2.2e-16Model yang dihasilkan adalah \[y_i=-2.849e+12+1.444e+09x_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.7288\). Artinya, sebesar 72.88% keragaman nilai GDP dapat dijelaskan oleh peubah tahun. Hasil ini menunjukkan hasil yang bagus, seolah mendapatkan hasil terbaik. Namun, kita perlu melakukan uji terhadap sisaannya seperti berikut ini.
#sisaan dan fitted value
sisaan<- residuals(model)
fitValue<- predict(model)
#Diagnostik dengan eksploratif
par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "steelblue", lwd = 2)
plot(fitValue, sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(sisaan, col = "steelblue")
plot(seq(1,64,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,64,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 di atas menunjukkan bahwa sisaan cenderung menyebar normal, tetapi histogram dari sisaan tidak menunjukkan demikian. Selanjutnya, dua plot di samping kanan digunakan untuk melihat autokorelasi. 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.
#Melihat Sisaan Menyebar Normal/Tidak
#H0: sisaan mengikuti sebaran normal
#H1: sisaan tidak mengikuti sebaran normal
shapiro.test(sisaan)##
## Shapiro-Wilk normality test
##
## data: sisaan
## W = 0.93854, p-value = 0.003246ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.1126, p-value = 0.3643
## alternative hypothesis: two-sidedBerdasarkan uji formal Saphiro-Wilk dan Kolmogorov-Smirnov didapatkan nilai p-value > \(\alpha\) (5%). Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal.
#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)
Berdasarkan plot ACF dan PACF, terlihat semua dalam rentang batas dan tidak ada lag yang signifikan. Namun, untuk lebih memastikan akan dilakukan uji formal dengan uji Durbin Watson.
#Deteksi autokorelasi dengan uji-Durbin Watson
#H0: tidak ada autokorelasi
#H1: ada autokorelasi
dwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 0.053876, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0Berdasarkan hasil DW Test, didapatkan nilai \(DW = 0.053876\) dan p-value = \(2.2e-16\). Berdasarkan tabel Durbin-Watson diperoleh nilai \(DL = 1.42\) dan \(DU = 1.58\). Nilai DW yang sangat jauh di bawah nilai DL dan DU. Artinya, menunjukkan bahwa model ini berada di daerah autokorelasi positif yang kuat. Dengan p-value yang sangat kecil (< 0.05), kita dapat menolak hipotesis nol, yang berarti ada cukup bukti untuk menyimpulkan adanya autokorelasi positif yang signifikan dalam residual model. Oleh karena itu, perlu dilakukan penanganan terhadap autokorelasi ini, misalnya dengan menggunakan metode Cochrane-Orcutt atau Hildreth-Lu untuk memperbaiki model dan menghilangkan pengaruh autokorelasi.
Penanganan Autokorelasi
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, max.iter = 1000)
modelCO## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = GDP ~ tahun, data = DataSL)
##
## number of interaction: 119
## rho 0.96566
##
## Durbin-Watson statistic
## (original): 0.05388 , p-value: 8.276e-46
## (transformed): 1.82806 , p-value: 2.063e-01
##
## coefficients:
## (Intercept) tahun
## -5.982145e+12 2.992340e+09Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=-5.982145e+12+2.992340e+09x_t\] asil juga menunjukkan bahwa nilai Durbin-Watson (DW) dan p-value meningkat secara signifikan, dengan DW berubah dari 0.05388 menjadi 1.82806 dan p-value meningkat menjadi 0.2063. Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.331 < DW < 2.669\). 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.96566\). Nilai tersebut dapat diketahui dengan syntax berikut.
#Rho optimum
rho<- modelCO$rho
rho## [1] 0.96566Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.
GDP## [1] 1409873950 1444327731 1434156379 1240672269 1309747899 1698319328
## [7] 1751470588 1859465021 1801344538 1965546218 2296470588 2369308600
## [13] 2553936348 2875625000 3574586466 3791298146 3591319857 4104509583
## [19] 2733183857 3364611432 4024621900 4415844156 4768765017 5167913302
## [25] 6043474843 5978460972 6405210564 6682167120 6978371581 6987267684
## [31] 8032551173 9000362582 9703011636 10338679636 11717604209 13029697561
## [37] 13897738375 15091913884 15794972847 15656327860 16330814180 15749753805
## [43] 16536535647 18881765437 20662525941 24405791045 28279802406 32350238760
## [49] 40713826215 42066224093 58636161082 67753284044 70447216891 77000578167
## [55] 82528535714 85140955389 88012282206 94376237797 94493871201 89014978344
## [61] 84304298771 88609332762 74144875098 84356860421GDP[-1]## [1] 1444327731 1434156379 1240672269 1309747899 1698319328 1751470588
## [7] 1859465021 1801344538 1965546218 2296470588 2369308600 2553936348
## [13] 2875625000 3574586466 3791298146 3591319857 4104509583 2733183857
## [19] 3364611432 4024621900 4415844156 4768765017 5167913302 6043474843
## [25] 5978460972 6405210564 6682167120 6978371581 6987267684 8032551173
## [31] 9000362582 9703011636 10338679636 11717604209 13029697561 13897738375
## [37] 15091913884 15794972847 15656327860 16330814180 15749753805 16536535647
## [43] 18881765437 20662525941 24405791045 28279802406 32350238760 40713826215
## [49] 42066224093 58636161082 67753284044 70447216891 77000578167 82528535714
## [55] 85140955389 88012282206 94376237797 94493871201 89014978344 84304298771
## [61] 88609332762 74144875098 84356860421#Transformasi Manual
GDP.trans<- GDP[-1]-GDP[-12]*rho
tahun.trans<- tahun[-1]-tahun[-12]*rho
modelCOmanual<- lm(GDP.trans~tahun.trans)
summary(modelCOmanual)##
## Call:
## lm(formula = GDP.trans ~ tahun.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -894667687 -464608685 -126640471 474242394 993575460
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.332e+11 9.921e+09 -13.43 <2e-16 ***
## tahun.trans 1.956e+09 1.447e+08 13.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 542500000 on 61 degrees of freedom
## Multiple R-squared: 0.7498, Adjusted R-squared: 0.7457
## F-statistic: 182.8 on 1 and 61 DF, p-value: < 2.2e-16Hasil 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)
## -3.879787e+12b1## tahun.trans
## 1955988878Hasil 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
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,1.5, by= 0.05))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4)## rho SSE
## 1 0.10 1.337501e+22
## 2 0.15 1.196781e+22
## 3 0.20 1.064429e+22
## 4 0.25 9.404476e+21
## 5 0.30 8.248356e+21
## 6 0.35 7.175934e+21
## 7 0.40 6.187208e+21
## 8 0.45 5.282178e+21
## 9 0.50 4.460846e+21
## 10 0.55 3.723210e+21
## 11 0.60 3.069270e+21
## 12 0.65 2.499028e+21
## 13 0.70 2.012482e+21
## 14 0.75 1.609633e+21
## 15 0.80 1.290481e+21
## 16 0.85 1.055025e+21
## 17 0.90 9.032659e+20
## 18 0.95 8.352036e+20
## 19 1.00 9.113604e+20
## 20 1.05 9.501692e+20
## 21 1.10 1.133197e+21
## 22 1.15 1.399922e+21
## 23 1.20 1.750343e+21
## 24 1.25 2.184461e+21
## 25 1.30 2.702276e+21
## 26 1.35 3.303787e+21
## 27 1.40 3.988995e+21
## 28 1.45 4.757900e+21
## 29 1.50 5.610502e+21Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.1. 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.2 sampai dengan 0.5.
#Rho optimal di sekitar 0.4
rOpt <- seq(0.5,1.0, 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
## 467 0.966 8.311004e+20
## 466 0.965 8.311058e+20
## 468 0.967 8.311286e+20
## 465 0.964 8.311446e+20
## 469 0.968 8.311902e+20
## 464 0.963 8.312169e+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.965, y=2.0e+21, labels = "rho=0.965", cex = 0.8)
Perhitungan yang dilakukan aplikasi R menunjukkan bahwa
nilai \(ρ\) optimum. 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.341, model)
summary(modelHL)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.615e+10 -9.980e+09 -2.183e+09 9.835e+09 2.218e+10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.958e+12 1.517e+11 -12.91 <2e-16 ***
## x 1.505e+09 1.155e+08 13.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.099e+10 on 61 degrees of freedom
## Multiple R-squared: 0.7355, Adjusted R-squared: 0.7312
## F-statistic: 169.7 on 1 and 61 DF, p-value: < 2.2e-16#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.341), "+", coef(modelHL)[2],"x", sep = "")## y = -2.971386e+12+1504521724xSetelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=-1.958e+12+1.505e+09x_t\]
#Deteksi autokorelasi
dwtest(modelHL)##
## Durbin-Watson test
##
## data: modelHL
## DW = 0.12874, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(0.12874\) berada jauh di bawah batas bawah DL dan jauh dari rentang yang menunjukkan tidak adanya autokorelasi, yaitu pada rentang DU < DW < 4-DU atau \(1.331 < DW < 2.669\). Hal tersebut juga didukung oleh p-value yang kecil dengan nilai \(2.2e-16\), di mana p-value < \(\alpha\)=5%. Oleh karena itu, hipotesis nol (H0) ditolak, yang berarti ada cukup bukti untuk menyimpulkan bahwa terdapat autokorelasi dalam data GDP dengan metode Hildreth-Lu pada taraf nyata 5%. Dengan adanya autokorelasi ini, diperlukan penanganan lebih lanjut untuk memperbaiki model dan mengatasi masalah autokorelasi.
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(GDP)
mseModelCO <- sseModelCO/length(GDP)
mseModelHL <- sseModelHL/length(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.693479e+22 1.795485e+19 7.362793e+21
## MSE 2.646060e+20 2.805446e+17 1.150436e+20Berdasarkan hasil tersebut, dapat diketahui bahwa penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memberikan hasil yang berbeda dalam hal Sum of Squared Errors (SSE) dan Mean Squared Error (MSE). Model Cochrane-Orcutt menghasilkan SSE sebesar 1.795485e+19 dan MSE sebesar 2.805446e+17, yang menunjukkan peningkatan signifikan dibandingkan model awal yang memiliki SSE sebesar 1.693479e+22 dan MSE sebesar 2.646060e+20. Model Hildreth-Lu, meskipun lebih baik daripada model awal, masih memiliki SSE sebesar 7.362793e+21dan MSE sebesar 1.150436e+20, yang lebih tinggi dibandingkan dengan model Cochrane-Orcutt. Ini menunjukkan bahwa metode Cochrane-Orcutt lebih efektif dalam mengurangi kesalahan dibandingkan dengan Hildreth-Lu, sehingga memberikan model yang lebih akurat setelah penanganan autokorelasi.