Tugas 2 Praktikum MPDW
Hafly Akeyla Pari - G1401221100
Pemanggilan Packages
# Memuat pustaka yang diperlukan
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(TTR)
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(lmtest) # digunakan untuk uji formal pendeteksian autokorelasi## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(orcutt) # untuk membuat model regresi Cochrane-Orcutt
library(HoRM) # untuk membuat model regresi Hildreth-Lu
library(readxl) # untuk membaca file ExcelInput Data
Data yang digunakan dalam kesempatan kali ini adalah data GDP Negara Madagaskar periode tahun 1960-2023.
Madagascar <- read_xlsx("C:/Users/Hafly Akeyla Pari/AppData/Downloads/Madagascar.xlsx")
# Menampilkan isi data yang diimpor untuk verifikasi
print(Madagascar)## # A tibble: 64 × 2
## Tahun GDP
## <chr> <dbl>
## 1 1960 673081725.
## 2 1961 699161945.
## 3 1962 739286908.
## 4 1963 759345864.
## 5 1964 802482184.
## 6 1965 833563473.
## 7 1966 900264585.
## 8 1967 956436932.
## 9 1968 1031669637.
## 10 1969 1056391056.
## # ℹ 54 more rowsMadagascar$Tahun <- as.numeric(Madagascar$Tahun)
Madagascar$GDP <- as.numeric(Madagascar$GDP)
Tahun <- Madagascar$Tahun
Tahun## [1] 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974
## [16] 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989
## [31] 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004
## [46] 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
## [61] 2020 2021 2022 2023GDP <- Madagascar$GDP
GDP## [1] 673081725 699161945 739286908 759345864 802482184 833563473
## [7] 900264585 956436932 1031669637 1056391056 1111859571 1199507631
## [13] 1341590690 1653062335 1917508190 2283049215 2181844179 2358930406
## [19] 2669755115 3463565854 5201818348 4759333998 4784977326 4686457031
## [25] 3905938481 3802557895 4347989788 3212900556 3189456965 3175638333
## [31] 3931334875 3254713056 3714966678 4063298919 3522227092 3838100904
## [37] 4931861039 4262965420 4401967633 4277903780 4629247204 5438332602
## [43] 5351701663 6372498890 5064732626 5859269753 6395712491 8524620890
## [49] 10725137724 9616879409 9982711338 11551819618 11578975062 12423555455
## [55] 12522957399 11323020701 11848613858 13176313594 13760033282 14104664679
## [61] 13051441204 14554754117 15302510500 16031702915# Menggabungkan kolom 'Tahun' dan 'GDP' dari data 'Madagascar'
Madagascar1 <- cbind(Tahun, GDP)
Madagascar1## Tahun GDP
## [1,] 1960 673081725
## [2,] 1961 699161945
## [3,] 1962 739286908
## [4,] 1963 759345864
## [5,] 1964 802482184
## [6,] 1965 833563473
## [7,] 1966 900264585
## [8,] 1967 956436932
## [9,] 1968 1031669637
## [10,] 1969 1056391056
## [11,] 1970 1111859571
## [12,] 1971 1199507631
## [13,] 1972 1341590690
## [14,] 1973 1653062335
## [15,] 1974 1917508190
## [16,] 1975 2283049215
## [17,] 1976 2181844179
## [18,] 1977 2358930406
## [19,] 1978 2669755115
## [20,] 1979 3463565854
## [21,] 1980 5201818348
## [22,] 1981 4759333998
## [23,] 1982 4784977326
## [24,] 1983 4686457031
## [25,] 1984 3905938481
## [26,] 1985 3802557895
## [27,] 1986 4347989788
## [28,] 1987 3212900556
## [29,] 1988 3189456965
## [30,] 1989 3175638333
## [31,] 1990 3931334875
## [32,] 1991 3254713056
## [33,] 1992 3714966678
## [34,] 1993 4063298919
## [35,] 1994 3522227092
## [36,] 1995 3838100904
## [37,] 1996 4931861039
## [38,] 1997 4262965420
## [39,] 1998 4401967633
## [40,] 1999 4277903780
## [41,] 2000 4629247204
## [42,] 2001 5438332602
## [43,] 2002 5351701663
## [44,] 2003 6372498890
## [45,] 2004 5064732626
## [46,] 2005 5859269753
## [47,] 2006 6395712491
## [48,] 2007 8524620890
## [49,] 2008 10725137724
## [50,] 2009 9616879409
## [51,] 2010 9982711338
## [52,] 2011 11551819618
## [53,] 2012 11578975062
## [54,] 2013 12423555455
## [55,] 2014 12522957399
## [56,] 2015 11323020701
## [57,] 2016 11848613858
## [58,] 2017 13176313594
## [59,] 2018 13760033282
## [60,] 2019 14104664679
## [61,] 2020 13051441204
## [62,] 2021 14554754117
## [63,] 2022 15302510500
## [64,] 2023 16031702915dtMadagascar <- as.data.frame(Madagascar1)
dtMadagascar## Tahun GDP
## 1 1960 673081725
## 2 1961 699161945
## 3 1962 739286908
## 4 1963 759345864
## 5 1964 802482184
## 6 1965 833563473
## 7 1966 900264585
## 8 1967 956436932
## 9 1968 1031669637
## 10 1969 1056391056
## 11 1970 1111859571
## 12 1971 1199507631
## 13 1972 1341590690
## 14 1973 1653062335
## 15 1974 1917508190
## 16 1975 2283049215
## 17 1976 2181844179
## 18 1977 2358930406
## 19 1978 2669755115
## 20 1979 3463565854
## 21 1980 5201818348
## 22 1981 4759333998
## 23 1982 4784977326
## 24 1983 4686457031
## 25 1984 3905938481
## 26 1985 3802557895
## 27 1986 4347989788
## 28 1987 3212900556
## 29 1988 3189456965
## 30 1989 3175638333
## 31 1990 3931334875
## 32 1991 3254713056
## 33 1992 3714966678
## 34 1993 4063298919
## 35 1994 3522227092
## 36 1995 3838100904
## 37 1996 4931861039
## 38 1997 4262965420
## 39 1998 4401967633
## 40 1999 4277903780
## 41 2000 4629247204
## 42 2001 5438332602
## 43 2002 5351701663
## 44 2003 6372498890
## 45 2004 5064732626
## 46 2005 5859269753
## 47 2006 6395712491
## 48 2007 8524620890
## 49 2008 10725137724
## 50 2009 9616879409
## 51 2010 9982711338
## 52 2011 11551819618
## 53 2012 11578975062
## 54 2013 12423555455
## 55 2014 12522957399
## 56 2015 11323020701
## 57 2016 11848613858
## 58 2017 13176313594
## 59 2018 13760033282
## 60 2019 14104664679
## 61 2020 13051441204
## 62 2021 14554754117
## 63 2022 15302510500
## 64 2023 16031702915Eksplorasi Data
#Membentuk objek time series
data.ts<-ts(dtMadagascar$GDP)
data.ts## Time Series:
## Start = 1
## End = 64
## Frequency = 1
## [1] 673081725 699161945 739286908 759345864 802482184 833563473
## [7] 900264585 956436932 1031669637 1056391056 1111859571 1199507631
## [13] 1341590690 1653062335 1917508190 2283049215 2181844179 2358930406
## [19] 2669755115 3463565854 5201818348 4759333998 4784977326 4686457031
## [25] 3905938481 3802557895 4347989788 3212900556 3189456965 3175638333
## [31] 3931334875 3254713056 3714966678 4063298919 3522227092 3838100904
## [37] 4931861039 4262965420 4401967633 4277903780 4629247204 5438332602
## [43] 5351701663 6372498890 5064732626 5859269753 6395712491 8524620890
## [49] 10725137724 9616879409 9982711338 11551819618 11578975062 12423555455
## [55] 12522957399 11323020701 11848613858 13176313594 13760033282 14104664679
## [61] 13051441204 14554754117 15302510500 16031702915#Membuat plot time series
ts.plot(data.ts, xlab="Tahun", ylab="GDP", main= "Time Series Plot of GDP")
points(data.ts)
Selanjutnya akan dilakukan ramalan dan pemulusan dengan metode DMA dan DES karena terlihat pada plot di atas menunjukkan adanya trend.
# Pemulusan dengan metode SMA dan DES
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)
# Ramalan ke depan
t = 1:5
f = numeric(length(t))
for (i in seq_along(t)) {
f[i] = At[length(At)] + Bt[length(Bt)] * i
}# Menggabungkan data untuk analisis
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)
)## Warning in cbind(aktual = c(data.ts, rep(NA, 5)), pemulusan1 = c(dt.sma, :
## number of rows of result is not a multiple of vector length (arg 1)print(dt.gab)## aktual pemulusan1 pemulusan2 At Bt ramalan
## [1,] 673081725 NA NA NA NA NA
## [2,] 699161945 NA NA NA NA NA
## [3,] 739286908 703843526 NA NA NA NA
## [4,] 759345864 732598239 NA NA NA NA
## [5,] 802482184 767038318 832128233 799583276 32544957 NA
## [6,] 833563473 798463840 863324589 830894215 32430374 832128233
## [7,] 900264585 845436747 929017638 887227192 41790445 863324589
## [8,] 956436932 896754997 996494600 946624799 49869802 929017638
## [9,] 1031669637 962790385 1085049735 1023920060 61129675 996494600
## [10,] 1056391056 1014832542 1128245676 1071539109 56706567 1085049735
## [11,] 1111859571 1066640088 1170411588 1118525838 51885750 1128245676
## [12,] 1199507631 1122586086 1231719114 1177152600 54566514 1170411588
## [13,] 1341590690 1217652631 1381705356 1299678993 82026363 1231719114
## [14,] 1653062335 1398053552 1701965811 1550009681 151956129 1381705356
## [15,] 1917508190 1637387072 2076765712 1857076392 219689320 1701965811
## [16,] 2283049215 1951206580 2529188271 2240197426 288990845 2076765712
## [17,] 2181844179 2127467195 2571694353 2349580774 222113579 2529188271
## [18,] 2358930406 2274607933 2588302661 2431455297 156847364 2571694353
## [19,] 2669755115 2403509900 2673473014 2538491457 134981557 2588302661
## [20,] 3463565854 2830750458 3486339181 3158544819 327794361 2673473014
## [21,] 5201818348 3778379772 5326712563 4552546168 774166395 3486339181
## [22,] 4759333998 4474906067 6035360668 5255133368 780227301 5326712563
## [23,] 4784977326 4915376557 5967021408 5441198983 525822425 6035360668
## [24,] 4686457031 4743589452 4808186972 4775888212 32298760 5967021408
## [25,] 3905938481 4459124279 3965312646 4212218463 -246905817 4808186972
## [26,] 3802557895 4131651136 3505376829 3818513982 -313137153 3965312646
## [27,] 4347989788 4018828721 3650083406 3834456064 -184372657 3505376829
## [28,] 3212900556 3787816080 3404584281 3596200180 -191615899 3650083406
## [29,] 3189456965 3583449103 3156951373 3370200238 -213248865 3404584281
## [30,] 3175638333 3192665285 2535375542 2864020413 -328644871 3156951373
## [31,] 3931334875 3432143391 3490924987 3461534189 29390798 2535375542
## [32,] 3254713056 3453895421 3642550199 3548222810 94327389 3490924987
## [33,] 3714966678 3633671536 3887874377 3760772957 127101420 3642550199
## [34,] 4063298919 3677659551 3856160981 3766910266 89250715 3887874377
## [35,] 3522227092 3766830897 3915051367 3840941132 74110235 3856160981
## [36,] 3838100904 3807875638 3922049524 3864962581 57086943 3915051367
## [37,] 4931861039 4097396345 4510787115 4304091730 206695385 3922049524
## [38,] 4262965420 4344309121 4866539960 4605424540 261115419 4510787115
## [39,] 4401967633 4532264697 4947480649 4739872673 207607976 4866539960
## [40,] 4277903780 4314278944 4148934991 4231606968 -82671976 4947480649
## [41,] 4629247204 4436372872 4453840941 4445106907 8734034 4148934991
## [42,] 5438332602 4781827862 5323830467 5052829165 271001302 4453840941
## [43,] 5351701663 5139760490 5847307320 5493533905 353773415 5323830467
## [44,] 6372498890 5720844385 6734244664 6227544524 506700139 5847307320
## [45,] 5064732626 5596311060 5817655890 5706983475 110672415 6734244664
## [46,] 5859269753 5765500423 5908064023 5836782223 71281800 5817655890
## [47,] 6395712491 5773238290 5896348355 5834793322 61555032 5908064023
## [48,] 8524620890 6926534378 8469421073 7697977725 771443348 5896348355
## [49,] 10725137724 8548490368 11479962414 10014226391 1465736023 8469421073
## [50,] 9616879409 9622212674 12135146409 10878679542 1256466868 11479962414
## [51,] 9982711338 10108242824 11472097894 10790170359 681927535 12135146409
## [52,] 11551819618 10383803455 11075237730 10729520593 345717137 11472097894
## [53,] 11578975062 11037835339 12093584939 11565710139 527874800 11075237730
## [54,] 12423555455 11851450045 13372290909 12611870477 760420432 12093584939
## [55,] 12522957399 12175162639 13149189234 12662175937 487013298 13372290909
## [56,] 11323020701 12089844519 12191895421 12140869970 51025451 13149189234
## [57,] 11848613858 11898197320 11585788974 11741993147 -156204173 12191895421
## [58,] 13176313594 12115982718 12278598449 12197290584 81307866 11585788974
## [59,] 13760033282 12928320245 14156627213 13542473729 614153484 12278598449
## [60,] 14104664679 13680337185 15224584789 14452460987 772123802 14156627213
## [61,] 13051441204 13638713055 14084558842 13861635948 222922893 15224584789
## [62,] 14554754117 13903620000 14229079839 14066349920 162729920 14084558842
## [63,] 15302510500 14302901940 15011882491 14657392215 354490275 14229079839
## [64,] 16031702915 15296322511 16887071231 16091696871 795374360 15011882491
## [65,] NA NA NA NA NA 16887071231
## [66,] NA NA NA NA NA 16887071231
## [67,] NA NA NA NA NA 17682445592
## [68,] NA NA NA NA NA 18477819952
## [69,] NA NA NA NA NA 19273194312
## [70,] 673081725 NA NA NA NA 20068568673# Plot time series dengan ramalan
ts.plot(dt.gab[,1], xlab="Tahun", ylab="GDP", main= "DMA N=3 Data GDP")
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 DMA
error.dma = data.ts-dt.ramal[1:length(data.ts)]
SSE.dma = sum(error.dma[32:length(data.ts)]^2)
MSE.dma = mean(error.dma[32:length(data.ts)]^2)
MAPE.dma = mean(abs((error.dma[32:length(data.ts)]/data.ts[32: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")
print(akurasi.dma)## Akurasi m = 3
## SSE 4.296413e+19
## MSE 1.301943e+18
## MAPE 1.026692e+01Selanjutnya akan digunakan metode Double Exponential Smoothing dengan cara sebagai berikut.
Pertama akan data akan dibagi menjadi data training dan data testing.
# Metode Double Exponential Smoothing (DES)
# Membagi data menjadi training dan testing
training <- dtMadagascar$GDP[1:51]
testing <- dtMadagascar$GDP[52:64]
# Data time series
training.ts <- ts(training, start = dtMadagascar$Tahun[1], frequency = 1)
testing.ts <- ts(testing, start = dtMadagascar$Tahun[52], frequency = 1)
# 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)
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.
# Pemulusan dengan DES dan peramalan
des.opt <- HoltWinters(training.ts, gamma = FALSE)
plot(des.opt)
legend("topleft", c("Data Aktual", "Peramalan"), col = c("black", "red"),
lty = c(1, 1))
# Ramalan
ramalandesopt <- forecast(des.opt, h=5)
print(ramalandesopt)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2011 10202227601 9304990315 11099464886 8830021587 11574433615
## 2012 10425931442 9174775389 11677087495 8512453382 12339409502
## 2013 10649635284 9106648971 12192621597 8289841356 13009429212
## 2014 10873339125 9069674800 12677003451 8114872581 13631805670
## 2015 11097042967 9051435342 13142650592 7968556078 14225529857Selanjutnya akan dicari akurasi dari metode DES.
# Akurasi metode DES untuk data training
ssedes.train <- des.opt$SSE
msedes.train <- ssedes.train / length(training.ts)
sisaandes <- ramalandesopt$residuals
mapedes.train <- sum(abs(sisaandes[!is.na(sisaandes)] / training.ts[!is.na(sisaandes)]) * 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")
print(akurasides.opt)## Akurasi lamda dan gamma optimum
## SSE 2.419479e+19
## MSE 4.744077e+17
## MAPE 1.013714e+01#Akurasi data testing
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt## Time Series:
## Start = 2011
## End = 2015
## Frequency = 1
## [1] -1349592017 -1153043620 -1773920172 -1649618274 -225977734SSEtestingdesopt<-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 6.976929e+17
## MSE 6.976929e+17
## MAPE 3.929861e+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 4.296413e+19 2.419479e+19
## MSE 1.301943e+18 4.744077e+17
## MAPE 1.026692e+01 1.013714e+01Berdasarkan 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
plot(Tahun, GDP, pch = 20, col = "blue",
main = "Scatter Plot Tahun vs Nilai GDP",
xlab = "Tahun",
ylab = "Nilai GDP")
#Menampilkan Nilai Korelasi
cor(Tahun,GDP)## [1] 0.9194147Berdasarkan scatter plot di atas, terlihat adanya hubungan / korelasi
positif antara peubah tahun dengan nilai IPM, 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.9194147\).
Setalah mengetahui adanya hubungan antar dua peubah, maka model regresi dapat ditentukan.
Regresi
#Pembuatan Model Regresi
#model regresi
model1<- lm(GDP~Tahun, data = dtMadagascar)
summary(model1)##
## Call:
## lm(formula = GDP ~ Tahun, data = dtMadagascar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.305e+09 -1.698e+09 2.086e+08 1.441e+09 3.467e+09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.340e+11 2.388e+10 -18.17 <2e-16 ***
## Tahun 2.208e+08 1.199e+07 18.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.772e+09 on 62 degrees of freedom
## Multiple R-squared: 0.8453, Adjusted R-squared: 0.8428
## F-statistic: 338.8 on 1 and 62 DF, p-value: < 2.2e-16Model yang dihasilkan adalah \[y_i=-434000000000+220800000x_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.8453\). Artinya, sebesar 84.53% 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(model1)
sisaan## 1 2 3 4 5 6
## 2016007098 1821334740 1640707125 1440013504 1262397246 1072725958
## 7 8 9 10 11 12
## 918674492 754094262 608574389 412543230 247259168 114154650
## 13 14 15 16 17 18
## 35485132 126204199 169897476 314685924 -7271690 -50938040
## 19 20 21 22 23 24
## 39134091 612192252 2129692169 1466455242 1271345992 952073119
## 25 26 27 28 29 30
## -49198009 -373331172 -48651857 -1404493667 -1648689835 -1883261045
## 31 32 33 34 35 36
## -1348317080 -2245691477 -2006190432 -1878610769 -2640435173 -2545313939
## 37 38 39 40 41 42
## -1672306382 -2561954579 -2643704943 -2988521373 -2857930527 -2269597707
## 43 44 45 46 47 48
## -2576981223 -1776936574 -3305455415 -2731670867 -2415980706 -507824885
## 49 50 51 52 53 54
## 1471939372 142928480 288007831 1636363533 1442766400 2066594215
## 55 56 57 58 59 60
## 1945243582 524554306 829394886 1936342043 2299309154 2423187973
## 61 62 63 64
## 1149211921 2431772256 2958776062 3467215899fitValue<- predict(model1)
#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)
Terlihat normal, tapi bisa terlihat polanya, artinya ada
autokorelasi
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.94629, p-value = 0.007545ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.14469, p-value = 0.1241
## 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)
Terdapat autokorelasi di amatan ke-0 saja, artinya dari eksplorasi ini
tidak terdapat auto korelasi
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(model1)##
## Durbin-Watson test
##
## data: model1
## DW = 0.16973, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0karena nilai p-value (\(< 2.2 \times 10^{-16}\)) < 0.05 maka terdapat autokorelasi dari uji durbin tersebut. Namun dalam tabel Durbinnya itu berada pada DL dan DU yang artinya INKONKLUSIF = Tidak tahu apakah benar-benar autokorelasi atau tidak
Berdasarkan hasil DW Test, didapatkan nilai \(DW = 0.16973\) dan p-value = \(< 2.2 \times 10^{-16}\). Berdasarkan tabel Durbin-Watson, diperoleh nilai \(DL = 1.5635\) dan \(DU = 1.6268\). Nilai DW masih berada di antara nilai DL dan DU. Artinya, berada di daerah inkonklusif, tidak dapat dikatakan berada di daerah autokorelasi positif maupun bebas dari autokorelasi. Dengan 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
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(model1)
modelCO## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = GDP ~ Tahun, data = dtMadagascar)
##
## number of interaction: 44
## rho 0.936607
##
## Durbin-Watson statistic
## (original): 0.16973 , p-value: 1.205e-28
## (transformed): 2.12800 , p-value: 6.477e-01
##
## coefficients:
## (Intercept) Tahun
## -689653286676 348292384Dari hasil tersebut dengan kohren orcutt dapat tertangani autokorelasinya
Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=-689653286676+348292384x_t\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(2.12800\) dan \(0.6477\). Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.6268 < DW < 2.3732\). 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.9366069\). Nilai tersebut dapat diketahui dengan syntax berikut.
#Rho optimum
rho<- modelCO$rho
rho## [1] 0.9366069Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.
GDP## [1] 673081725 699161945 739286908 759345864 802482184 833563473
## [7] 900264585 956436932 1031669637 1056391056 1111859571 1199507631
## [13] 1341590690 1653062335 1917508190 2283049215 2181844179 2358930406
## [19] 2669755115 3463565854 5201818348 4759333998 4784977326 4686457031
## [25] 3905938481 3802557895 4347989788 3212900556 3189456965 3175638333
## [31] 3931334875 3254713056 3714966678 4063298919 3522227092 3838100904
## [37] 4931861039 4262965420 4401967633 4277903780 4629247204 5438332602
## [43] 5351701663 6372498890 5064732626 5859269753 6395712491 8524620890
## [49] 10725137724 9616879409 9982711338 11551819618 11578975062 12423555455
## [55] 12522957399 11323020701 11848613858 13176313594 13760033282 14104664679
## [61] 13051441204 14554754117 15302510500 16031702915GDP[-1]## [1] 699161945 739286908 759345864 802482184 833563473 900264585
## [7] 956436932 1031669637 1056391056 1111859571 1199507631 1341590690
## [13] 1653062335 1917508190 2283049215 2181844179 2358930406 2669755115
## [19] 3463565854 5201818348 4759333998 4784977326 4686457031 3905938481
## [25] 3802557895 4347989788 3212900556 3189456965 3175638333 3931334875
## [31] 3254713056 3714966678 4063298919 3522227092 3838100904 4931861039
## [37] 4262965420 4401967633 4277903780 4629247204 5438332602 5351701663
## [43] 6372498890 5064732626 5859269753 6395712491 8524620890 10725137724
## [49] 9616879409 9982711338 11551819618 11578975062 12423555455 12522957399
## [55] 11323020701 11848613858 13176313594 13760033282 14104664679 13051441204
## [61] 14554754117 15302510500 16031702915ipm-1 artinya datanya dimundurin satu
#Transformasi Manual
ipm.trans<- GDP[-1]-GDP[-64]*rho
tahun.trans<- Tahun[-1]-Tahun[-64]*rho
modelCOmanual<- lm(ipm.trans~tahun.trans)
summary(modelCOmanual)##
## Call:
## lm(formula = ipm.trans ~ tahun.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.758e+09 -2.807e+08 7.072e+07 2.647e+08 1.799e+09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.372e+10 9.882e+09 -4.424 4.07e-05 ***
## tahun.trans 3.483e+08 7.768e+07 4.484 3.30e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 710700000 on 61 degrees of freedom
## Multiple R-squared: 0.2479, Adjusted R-squared: 0.2356
## F-statistic: 20.11 on 1 and 61 DF, p-value: 3.295e-05Hasil 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)
## -689653286676b1## tahun.trans
## 348292384Hasil 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.0, by= 0.1))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model1))}))
round(tab, 4)## rho SSE
## 1 0.1 1.581543e+20
## 2 0.2 1.295309e+20
## 3 0.3 1.045465e+20
## 4 0.4 8.320081e+19
## 5 0.5 6.549400e+19
## 6 0.6 5.142602e+19
## 7 0.7 4.099687e+19
## 8 0.8 3.420655e+19
## 9 0.9 3.105507e+19
## 10 1.0 3.302031e+19“rho” = r, “SSE” = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}) fungsi untuk menghitung SSE dari setiap R
Pertama-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.0.
#Rho optimal di sekitar 0.9
rOpt <- seq(0.8,1.0, by= 0.001)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model1))}))
head(tabOpt[order(tabOpt$SSE),])## rho SSE
## 138 0.937 3.081128e+19
## 137 0.936 3.081132e+19
## 139 0.938 3.081161e+19
## 136 0.935 3.081173e+19
## 140 0.939 3.081230e+19
## 135 0.934 3.081249e+19#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.9366069, y=30811260000000000000, labels = "rho=0.9366069", cex = 0.8)
Perhitungan yang dilakukan aplikasi R menunjukkan bahwa
nilai \(ρ\) optimum, yaitu saat SSE
terkecil terdapat pada nilai \(ρ=0.9366069\). 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.9366069, model1)
summary(modelHL)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.758e+09 -2.807e+08 7.072e+07 2.647e+08 1.799e+09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.372e+10 9.882e+09 -4.424 4.07e-05 ***
## x 3.483e+08 7.768e+07 4.484 3.30e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 710700000 on 61 degrees of freedom
## Multiple R-squared: 0.2479, Adjusted R-squared: 0.2356
## F-statistic: 20.11 on 1 and 61 DF, p-value: 3.295e-05#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.9366069), "+", coef(modelHL)[2],"x", sep = "")## y = -689653134862+348292308xSetelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=-689653134862+348292308x_t\]
#Deteksi autokorelasi
dwtest(modelHL)##
## Durbin-Watson test
##
## data: modelHL
## DW = 2.128, p-value = 0.6477
## alternative hypothesis: true autocorrelation is greater than 0Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(2.128\) berada pada selang daerah tidak ada autokorelasi, yaitu pada rentang DU < DW < 4-DU atau \(1.6268 < DW < 2.3732\). Hal tersebut juga didukung oleh p-value sebesar \(0.6477\), di mana p-value > \(\alpha\)=5%. Artinya tak tolak \(H_0\) atau tidak cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai GDP dengan metode Hildreth-Lu pada taraf nyata 5%.
Terakhir, bandingkan nilai SSE dari ketiga metode (metode awal, metode Cochrane-Orcutt, dan Hildreth-Lu).
#Perbandingan
sseModelawal <- anova(model1)$`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.94745e+20 3.081126e+19 3.081126e+19
## MSE 3.04289e+18 4.814259e+17 4.814259e+17Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hidreth-Lu memiliki SSE yang sama, sebesar \(3.081126 \times 10^{19}\) dan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar \(1.94745 \times 10^{20}\).
Simpulan
Autokorelasi yang terdapat pada data GDP terjadi akibat adanya korelasi di antara unsur penyusunnya. Nilai 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.