Tugas MPDW ke-2

2024-09-01

Muhammad Jodi At-Takbir

Inisialisasi Library

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(TTR)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(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.numeric
library(orcutt)   #untuk membuat model regresi Cochrane-Orcutt
## Warning: package 'orcutt' was built under R version 4.4.1
library(HoRM)     #untuk membuat model regresi Hildreth-Lu
## Warning: package 'HoRM' was built under R version 4.4.1
library(readxl)

Input Data

Data yang digunakan adalah data GDP dari East Asia & Pacific

dataGDP <- read_excel("C:/Users/jodia/Downloads/Data GDP.xlsx", col_names = TRUE)
dataGDP
## # A tibble: 64 × 2
##    Periode           GDP
##      <dbl>         <dbl>
##  1    1960  80767959484.
##  2    1961  71153218973.
##  3    1962  64897470662.
##  4    1963  70193946313.
##  5    1964  81084335638.
##  6    1965  94481207865.
##  7    1966 103367734664.
##  8    1967 100363074006.
##  9    1968 101465070644.
## 10    1969 114096748857.
## # ℹ 54 more rows

Data yang ingin dianalisis adalah data dengan periode 1960 sampai 2010

dataGDP2 <- dataGDP %>% slice(-c(52:64))
dataGDP2
## # A tibble: 51 × 2
##    Periode           GDP
##      <dbl>         <dbl>
##  1    1960  80767959484.
##  2    1961  71153218973.
##  3    1962  64897470662.
##  4    1963  70193946313.
##  5    1964  81084335638.
##  6    1965  94481207865.
##  7    1966 103367734664.
##  8    1967 100363074006.
##  9    1968 101465070644.
## 10    1969 114096748857.
## # ℹ 41 more rows

Eksplorasi Data

Sebelum melakukan regresi, akan diperlihatkan plot time-series dari GDP East Asia & Pacific

#Membentuk objek time series
data.ts<-ts(dataGDP2$GDP)
data.ts
## Time Series:
## Start = 1 
## End = 51 
## Frequency = 1 
##  [1] 8.076796e+10 7.115322e+10 6.489747e+10 7.019395e+10 8.108434e+10
##  [6] 9.448121e+10 1.033677e+11 1.003631e+11 1.014651e+11 1.140967e+11
## [11] 1.270559e+11 1.364474e+11 1.553651e+11 1.954342e+11 2.214909e+11
## [16] 2.485678e+11 2.527368e+11 2.912998e+11 2.827301e+11 3.273617e+11
## [21] 3.776751e+11 4.032023e+11 4.234608e+11 4.433967e+11 4.804188e+11
## [26] 5.260231e+11 5.246301e+11 5.185505e+11 5.755909e+11 6.222231e+11
## [31] 6.667072e+11 7.235164e+11 8.114062e+11 8.906182e+11 1.071871e+12
## [36] 1.322031e+12 1.518329e+12 1.570609e+12 1.433203e+12 1.577918e+12
## [41] 1.736143e+12 1.847120e+12 2.044360e+12 2.313863e+12 2.682937e+12
## [46] 3.107695e+12 3.742774e+12 4.731773e+12 5.987096e+12 6.491259e+12
## [51] 7.900542e+12
#Membuat plot time series
ts.plot(data.ts, xlab="Time Period ", 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.

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)
}
#Melakukan pemulusan
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
##  [1,] 8.076796e+10           NA           NA           NA           NA
##  [2,] 7.115322e+10           NA           NA           NA           NA
##  [3,] 6.489747e+10 7.227288e+10           NA           NA           NA
##  [4,] 7.019395e+10 6.874821e+10           NA           NA           NA
##  [5,] 8.108434e+10 7.205858e+10 7.412263e+10 7.309061e+10 1.032024e+09
##  [6,] 9.448121e+10 8.191983e+10 9.727507e+10 8.959745e+10 7.677621e+09
##  [7,] 1.033677e+11 9.297776e+10 1.142958e+11 1.036368e+11 1.065903e+10
##  [8,] 1.003631e+11 9.940401e+10 1.153443e+11 1.073741e+11 7.970141e+09
##  [9,] 1.014651e+11 1.017320e+11 1.091201e+11 1.054260e+11 3.694052e+09
## [10,] 1.140967e+11 1.053083e+11 1.116287e+11 1.084685e+11 3.160210e+09
## [11,] 1.270559e+11 1.142059e+11 1.284536e+11 1.213298e+11 7.123858e+09
## [12,] 1.364474e+11 1.258667e+11 1.473461e+11 1.366064e+11 1.073972e+10
## [13,] 1.553651e+11 1.396228e+11 1.657381e+11 1.526805e+11 1.305767e+10
## [14,] 1.954342e+11 1.624156e+11 2.019766e+11 1.821961e+11 1.978054e+10
## [15,] 2.214909e+11 1.907634e+11 2.437557e+11 2.172595e+11 2.649614e+10
## [16,] 2.485678e+11 2.218309e+11 2.821529e+11 2.519919e+11 3.016098e+10
## [17,] 2.527368e+11 2.409318e+11 2.871114e+11 2.640216e+11 2.308977e+10
## [18,] 2.912998e+11 2.642015e+11 3.079616e+11 2.860815e+11 2.188005e+10
## [19,] 2.827301e+11 2.755889e+11 3.062853e+11 2.909371e+11 1.534818e+10
## [20,] 3.273617e+11 3.004639e+11 3.412221e+11 3.208430e+11 2.037913e+10
## [21,] 3.776751e+11 3.292556e+11 3.842280e+11 3.567418e+11 2.748617e+10
## [22,] 4.032023e+11 3.694130e+11 4.421507e+11 4.057819e+11 3.636885e+10
## [23,] 4.234608e+11 4.014461e+11 4.709283e+11 4.361872e+11 3.474115e+10
## [24,] 4.433967e+11 4.233532e+11 4.739182e+11 4.486357e+11 2.528246e+10
## [25,] 4.804188e+11 4.490921e+11 4.980153e+11 4.735537e+11 2.446162e+10
## [26,] 5.260231e+11 4.832795e+11 5.460220e+11 5.146508e+11 3.137125e+10
## [27,] 5.246301e+11 5.103573e+11 5.692527e+11 5.398050e+11 2.944769e+10
## [28,] 5.185505e+11 5.230679e+11 5.580672e+11 5.405676e+11 1.749966e+10
## [29,] 5.755909e+11 5.395905e+11 5.700944e+11 5.548425e+11 1.525193e+10
## [30,] 6.222231e+11 5.721215e+11 6.265113e+11 5.993164e+11 2.719487e+10
## [31,] 6.667072e+11 6.215071e+11 7.090418e+11 6.652744e+11 4.376737e+10
## [32,] 7.235164e+11 6.708156e+11 7.694839e+11 7.201497e+11 4.933417e+10
## [33,] 8.114062e+11 7.338766e+11 8.508303e+11 7.923535e+11 5.847686e+10
## [34,] 8.906182e+11 8.085136e+11 9.500703e+11 8.792920e+11 7.077835e+10
## [35,] 1.071871e+12 9.246318e+11 1.129214e+12 1.026923e+12 1.022911e+11
## [36,] 1.322031e+12 1.094840e+12 1.399197e+12 1.247018e+12 1.521783e+11
## [37,] 1.518329e+12 1.304077e+12 1.696532e+12 1.500305e+12 1.962275e+11
## [38,] 1.570609e+12 1.470323e+12 1.831476e+12 1.650899e+12 1.805763e+11
## [39,] 1.433203e+12 1.507380e+12 1.667621e+12 1.587501e+12 8.012018e+10
## [40,] 1.577918e+12 1.527243e+12 1.578432e+12 1.552838e+12 2.559439e+10
## [41,] 1.736143e+12 1.582421e+12 1.669234e+12 1.625828e+12 4.340641e+10
## [42,] 1.847120e+12 1.720394e+12 1.941142e+12 1.830768e+12 1.103741e+11
## [43,] 2.044360e+12 1.875874e+12 2.175163e+12 2.025519e+12 1.496445e+11
## [44,] 2.313863e+12 2.068448e+12 2.428866e+12 2.248657e+12 1.802091e+11
## [45,] 2.682937e+12 2.347053e+12 2.846910e+12 2.596982e+12 2.499283e+11
## [46,] 3.107695e+12 2.701498e+12 3.359829e+12 3.030663e+12 3.291652e+11
## [47,] 3.742774e+12 3.177802e+12 4.049170e+12 3.613486e+12 4.356842e+11
## [48,] 4.731773e+12 3.860747e+12 5.088877e+12 4.474812e+12 6.140649e+11
## [49,] 5.987096e+12 4.820548e+12 6.555579e+12 5.688063e+12 8.675155e+11
## [50,] 6.491259e+12 5.736709e+12 7.598125e+12 6.667417e+12 9.307078e+11
## [51,] 7.900542e+12 6.792965e+12 8.812081e+12 7.802523e+12 1.009558e+12
## [52,]           NA           NA           NA           NA           NA
## [53,]           NA           NA           NA           NA           NA
## [54,]           NA           NA           NA           NA           NA
## [55,]           NA           NA           NA           NA           NA
## [56,]           NA           NA           NA           NA           NA
##            ramalan
##  [1,]           NA
##  [2,]           NA
##  [3,]           NA
##  [4,]           NA
##  [5,]           NA
##  [6,] 7.412263e+10
##  [7,] 9.727507e+10
##  [8,] 1.142958e+11
##  [9,] 1.153443e+11
## [10,] 1.091201e+11
## [11,] 1.116287e+11
## [12,] 1.284536e+11
## [13,] 1.473461e+11
## [14,] 1.657381e+11
## [15,] 2.019766e+11
## [16,] 2.437557e+11
## [17,] 2.821529e+11
## [18,] 2.871114e+11
## [19,] 3.079616e+11
## [20,] 3.062853e+11
## [21,] 3.412221e+11
## [22,] 3.842280e+11
## [23,] 4.421507e+11
## [24,] 4.709283e+11
## [25,] 4.739182e+11
## [26,] 4.980153e+11
## [27,] 5.460220e+11
## [28,] 5.692527e+11
## [29,] 5.580672e+11
## [30,] 5.700944e+11
## [31,] 6.265113e+11
## [32,] 7.090418e+11
## [33,] 7.694839e+11
## [34,] 8.508303e+11
## [35,] 9.500703e+11
## [36,] 1.129214e+12
## [37,] 1.399197e+12
## [38,] 1.696532e+12
## [39,] 1.831476e+12
## [40,] 1.667621e+12
## [41,] 1.578432e+12
## [42,] 1.669234e+12
## [43,] 1.941142e+12
## [44,] 2.175163e+12
## [45,] 2.428866e+12
## [46,] 2.846910e+12
## [47,] 3.359829e+12
## [48,] 4.049170e+12
## [49,] 5.088877e+12
## [50,] 6.555579e+12
## [51,] 7.598125e+12
## [52,] 8.812081e+12
## [53,] 9.821639e+12
## [54,] 1.083120e+13
## [55,] 1.184075e+13
## [56,] 1.285031e+13
#Plot time series
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="GDP", 
        main= "DMA N=3 Data GDP", ylim=c(64897470662
,7.90054E+12))
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.002004e+24
## MSE   4.352183e+22
## MAPE  8.236755e+00

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<-dataGDP2[1:40,2]
testing<-dataGDP2[41:51,2]

#data time series
training.ts<-ts(training)
testing.ts<-ts(testing,start=41)

#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.

#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.2046693
##  gamma: FALSE
## 
## Coefficients:
##           [,1]
## a 1.577918e+12
## b 8.032507e+10
plot(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
## 41   1.658243e+12 1.580069e+12 1.736418e+12 1.538685e+12 1.777801e+12
## 42   1.738568e+12 1.616175e+12 1.860962e+12 1.551384e+12 1.925753e+12
## 43   1.818893e+12 1.654216e+12 1.983570e+12 1.567041e+12 2.070745e+12
## 44   1.899218e+12 1.691761e+12 2.106676e+12 1.581940e+12 2.216497e+12
## 45   1.979543e+12 1.728044e+12 2.231043e+12 1.594908e+12 2.364179e+12

Selanjutnya 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,]  3358992200
## [4,] 14223733559
## [5,] 16906485550
## [6,] 15952729780
mapedes.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                     1.427587e+23
## MSE                     3.568967e+21
## MAPE                    6.872863e+00
#Akurasi data testing
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt
## Time Series:
## Start = 41 
## End = 45 
## Frequency = 1 
##         testing.ts
## [1,]  -77899827252
## [2,] -108551295398
## [3,] -225467006987
## [4,] -414644259650
## [5,] -703393630476
SSEtestingdesopt<-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.685269e+22
## MSE                     6.685269e+22
## MAPE                    5.957252e+00

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.002004e+24                    1.427587e+23
## MSE   4.352183e+22                    3.568967e+21
## MAPE  8.236755e+00                    6.872863e+00

Berdasarkan 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(dataGDP2$Periode, dataGDP2$GDP, pch = 20, col = "blue",
     main = "Scatter Plot Periode vs GDP",
     xlab = "Periode",
     ylab = "GDP")

#Menampilkan Nilai Korelasi
cor(dataGDP2$Periode,dataGDP2$GDP)
## [1] 0.7719426

Berdasarkan scatter plot di atas, terlihat adanya hubungan / korelasi positif antara peubah Periode dengan 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.7719426.

Setalah mengetahui adanya hubungan antar dua peubah, maka model regresi dapat ditentukan.

Regresi

#Pembuatan Model Regresi
#model regresi
model<- lm(GDP~Periode, data = dataGDP2)
summary(model)
## 
## Call:
## lm(formula = GDP ~ Periode, data = dataGDP2)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.052e+12 -7.931e+11 -3.039e+11  4.633e+11  4.422e+12 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.782e+14  2.111e+13  -8.442 4.06e-11 ***
## Periode      9.039e+10  1.063e+10   8.500 3.32e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.118e+12 on 49 degrees of freedom
## Multiple R-squared:  0.5959, Adjusted R-squared:  0.5876 
## F-statistic: 72.26 on 1 and 49 DF,  p-value: 3.315e-11

Model yang dihasilkan adalah yi = (−1.782e+14) + (9.039e+10)xt Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value < α (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 < α (5%) sehingga nyata dalam taraf 5%. Selanjutnya dapat dilihat juga nilai R2 = 0.5959. Artinya, sebesar 59.59% keragaman GDP dapat dijelaskan oleh peubah periode. Hasil ini menunjukkan hasil yang cukup baik. 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,51,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,51,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 tidak menyebar normal, diperkuat dengan 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
ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.17091, p-value = 0.09004
## alternative hypothesis: two-sided

Berdasarkan uji formal Kolmogorov-smirnov, diadapatkan p-value > α (5%). Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal walaupun pada histogram dan plot sebelumnya tidak cukup untuk mengatakan demikian.

#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)

Berdasarkan plot ACF dan PACF, terlihat dalam rentang batas terdapat 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.079624, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil DW Test, diperoleh 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(model)
modelCO
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = GDP ~ Periode, data = dataGDP2)
## 
##  number of interaction: 33
##  rho 1.212989
## 
## Durbin-Watson statistic 
## (original):    0.07962 , p-value: 8.874e-33
## (transformed): 2.13432 , p-value: 6.288e-01
##  
##  coefficients: 
##   (Intercept)       Periode 
## -3.098440e+13  1.582774e+10

Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. yi =  − 3.098440e + 13 + 1.582774e + 10xt Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi 2.13432 dan 0.62881. 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 1.212989. Nilai tersebut dapat diketahui dengan syntax berikut.

#Rho optimum
rho<- modelCO$rho
rho
## [1] 1.212989

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

GDP <- dataGDP2$GDP
GDP
##  [1] 8.076796e+10 7.115322e+10 6.489747e+10 7.019395e+10 8.108434e+10
##  [6] 9.448121e+10 1.033677e+11 1.003631e+11 1.014651e+11 1.140967e+11
## [11] 1.270559e+11 1.364474e+11 1.553651e+11 1.954342e+11 2.214909e+11
## [16] 2.485678e+11 2.527368e+11 2.912998e+11 2.827301e+11 3.273617e+11
## [21] 3.776751e+11 4.032023e+11 4.234608e+11 4.433967e+11 4.804188e+11
## [26] 5.260231e+11 5.246301e+11 5.185505e+11 5.755909e+11 6.222231e+11
## [31] 6.667072e+11 7.235164e+11 8.114062e+11 8.906182e+11 1.071871e+12
## [36] 1.322031e+12 1.518329e+12 1.570609e+12 1.433203e+12 1.577918e+12
## [41] 1.736143e+12 1.847120e+12 2.044360e+12 2.313863e+12 2.682937e+12
## [46] 3.107695e+12 3.742774e+12 4.731773e+12 5.987096e+12 6.491259e+12
## [51] 7.900542e+12
Periode <- dataGDP2$Periode
Periode
##  [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
GDP[-1]
##  [1] 7.115322e+10 6.489747e+10 7.019395e+10 8.108434e+10 9.448121e+10
##  [6] 1.033677e+11 1.003631e+11 1.014651e+11 1.140967e+11 1.270559e+11
## [11] 1.364474e+11 1.553651e+11 1.954342e+11 2.214909e+11 2.485678e+11
## [16] 2.527368e+11 2.912998e+11 2.827301e+11 3.273617e+11 3.776751e+11
## [21] 4.032023e+11 4.234608e+11 4.433967e+11 4.804188e+11 5.260231e+11
## [26] 5.246301e+11 5.185505e+11 5.755909e+11 6.222231e+11 6.667072e+11
## [31] 7.235164e+11 8.114062e+11 8.906182e+11 1.071871e+12 1.322031e+12
## [36] 1.518329e+12 1.570609e+12 1.433203e+12 1.577918e+12 1.736143e+12
## [41] 1.847120e+12 2.044360e+12 2.313863e+12 2.682937e+12 3.107695e+12
## [46] 3.742774e+12 4.731773e+12 5.987096e+12 6.491259e+12 7.900542e+12
#Transformasi Manual
GDP.trans<- GDP[-1]-GDP[-51]*rho
Periode.trans<- Periode[-1]-Periode[-51]*rho
modelCOmanual<- lm(GDP.trans~Periode.trans)
summary(modelCOmanual)
## 
## Call:
## lm(formula = GDP.trans ~ Periode.trans)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.169e+11 -1.892e+10  5.650e+09  2.386e+10  3.982e+11 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   6.599e+12  2.688e+12   2.455   0.0177 *
## Periode.trans 1.583e+10  6.374e+09   2.483   0.0166 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.385e+11 on 48 degrees of freedom
## Multiple R-squared:  0.1139, Adjusted R-squared:  0.09539 
## F-statistic: 6.167 on 1 and 48 DF,  p-value: 0.01656

Hasil 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.09844e+13
b1
## Periode.trans 
##   15827742883

Hasil 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, 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 5.054133e+25
## 2  0.2 4.202531e+25
## 3  0.3 3.431043e+25
## 4  0.4 2.739669e+25
## 5  0.5 2.128409e+25
## 6  0.6 1.597263e+25
## 7  0.7 1.146231e+25
## 8  0.8 7.753124e+24
## 9  0.9 4.845079e+24
## 10 1.0 4.656627e+24

Pertama-tama akan dicari di mana kira-kira ρ yang menghasilkan SSE minimum. Pada hasil di atas terlihat ρ minimum yang mungkin, yaitu 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.2 sampai dengan 0.5.

#Rho optimal di sekitar 0.9
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
## 200 0.999 2.755277e+24
## 199 0.998 2.772461e+24
## 198 0.997 2.789724e+24
## 197 0.996 2.807068e+24
## 196 0.995 2.824492e+24
## 195 0.994 2.841996e+24
#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.94, y=2e+25, labels = "rho=0.999", cex = 0.8)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai ρ optimum, yaitu saat SSE terkecil terdapat pada nilai ρ = 0.999. 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.999, model)
summary(modelHL)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -4.640e+11 -1.517e+11 -3.385e+10  7.720e+10  9.238e+11 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.059e+13  7.007e+12  -5.792 5.19e-07 ***
## x            1.365e+13  2.348e+12   5.815 4.80e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.396e+11 on 48 degrees of freedom
## Multiple R-squared:  0.4133, Adjusted R-squared:  0.4011 
## F-statistic: 33.81 on 1 and 48 DF,  p-value: 4.8e-07
#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.999), "+", coef(modelHL)[2],"x", sep = "")
## y = -4.059065e+16+1.365325e+13x

Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. yi = (−4.059065e+16) + (1.365325e+13)xt

#Deteksi autokorelasi
dwtest(modelHL)
## 
##  Durbin-Watson test
## 
## data:  modelHL
## DW = 0.65087, p-value = 5.494e-09
## alternative hypothesis: true autocorrelation is greater than 0

Setelah dilakukan penanganan dengan metode Hildreth-Lu didapatkan nilai p-value yang masih kurang dari 0.05, sehingga masih terjadi 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 6.121887e+25          9.210230e+23      2.755277e+24
## MSE 1.200370e+24          1.805927e+22      5.402504e+22

Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt lebih baik dibandingkan dengan metode Hildreth-Lu. Hal ini ditunjukkan dengan nilai SSE metode Cochrane-Orcutt lebih kecil dibandingkan dengan metode Hildreth-Lu, dan lebih baik pula dibandingkan model awal ketika autokorelasi masih terjadi.

Simpulan

Autokorelasi yang terdapat pada data GDP East Asia & Pacific terjadi akibat adanya korelasi di antara unsur penyusunnya. Indikator 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. Autokorelasi tersebut dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Namun, dalam penerapan data GDP East Asia & Pacific, metode Cochrane-Orcutt dapat menangani autokorelasi lebih baik dibandingkan dengan metode Hildreth-Lu yang ditunjukkan dengan SSE yang lebih kecil dan p-value yang lebih besar dari 0.05.