Tugas Individu Pertemuan 2

Nama : Ghonniyu Hiban Saputra NIM : G1401221012

Input Data

Data yang digunakan dalam kesempatan kali ini adalah data GDP Africa Eastern and Southern periode tahun 1960-1999 (40 Tahun).

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
library(HoRM) #untuk membuat model regresi Hildreth-Lu

dtafrica <- read.csv("C:/Users/Ghonniyu/Downloads/Data Individu Pert 2.csv")
dtafrica$Tanggal <- as.double(dtafrica$Tanggal)
dtafrica

##    Tanggal          GDP
## 1     1960  21216962290
## 2     1961  22307471356
## 3     1962  23702472100
## 4     1963  25779376633
## 5     1964  28049537325
## 6     1965  30374910060
## 7     1966  33049155473
## 8     1967  35933757402
## 9     1968  38749864977
## 10    1969  42964345101
## 11    1970  43702318276
## 12    1971  47632993913
## 13    1972  51627511757
## 14    1973  66610917871
## 15    1974  81658052383
## 16    1975  86883552689
## 17    1976  88290755581
## 18    1977  98868135190
## 19    1978 110802000000
## 20    1979 129781000000
## 21    1980 143250000000
## 22    1981 154793000000
## 23    1982 166208000000
## 24    1983 168165000000
## 25    1984 169303000000
## 26    1985 177750000000
## 27    1986 179898000000
## 28    1987 185812000000
## 29    1988 202093000000
## 30    1989 214473000000
## 31    1990 239433000000
## 32    1991 240271000000
## 33    1992 243565000000
## 34    1993 251212000000
## 35    1994 260992000000
## 36    1995 265429000000
## 37    1996 267815000000
## 38    1997 268834000000
## 39    1998 272519000000
## 40    1999 273433000000

Eksplorasi Data

Sebelum melakukan regresi, akan diperlihatkan plot time-series dari GDP Africa Eastern and Southern periode 1960-1999 ## Membuat plot time series

#Membentuk objek time series
data.ts<-ts(dtafrica$GDP)
data.ts

## Time Series:
## Start = 1 
## End = 40 
## Frequency = 1 
##  [1]  21216962290  22307471356  23702472100  25779376633  28049537325
##  [6]  30374910060  33049155473  35933757402  38749864977  42964345101
## [11]  43702318276  47632993913  51627511757  66610917871  81658052383
## [16]  86883552689  88290755581  98868135190 110802000000 129781000000
## [21] 143250000000 154793000000 166208000000 168165000000 169303000000
## [26] 177750000000 179898000000 185812000000 202093000000 214473000000
## [31] 239433000000 240271000000 243565000000 251212000000 260992000000
## [36] 265429000000 267815000000 268834000000 272519000000 273433000000

#Membuat plot time series
ts.plot(data.ts, xlab="Time Period ", ylab="GDP", main= "Time Series Plot of GDP")
points(data.ts)

Penentuan Metode

Selanjutnya akan dilakukan ramalan dan pemulusan dengan metode DMA dan DES karena terlihat pada plot di atas menunjukkan adanya trend.

Metode DMA

dt.sma <- SMA(data.ts, n=4)
dma <- SMA(dt.sma, n = 4)
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)
}

Plot DMA

dt.gab <- cbind(aktual = c(data.ts,rep(NA,7)), #7 didapatkan dati n=4 diatas 2n-1
                pemulusan1 = c(dt.sma,rep(NA,7)),
                pemulusan2 = c(dt.dma, rep(NA,7)),
                At = c(At, rep(NA,7)), 
                Bt = c(Bt,rep(NA,7)),
                ramalan = c(dt.ramal, f[-1]))

## Warning in cbind(aktual = c(data.ts, rep(NA, 7)), pemulusan1 = c(dt.sma, :
## number of rows of result is not a multiple of vector length (arg 6)

dt.gab

##             aktual   pemulusan1   pemulusan2           At          Bt
##  [1,]  21216962290           NA           NA           NA          NA
##  [2,]  22307471356           NA           NA           NA          NA
##  [3,]  23702472100           NA           NA           NA          NA
##  [4,]  25779376633  23251570595           NA           NA          NA
##  [5,]  28049537325  24959714354           NA           NA          NA
##  [6,]  30374910060  26976574030           NA           NA          NA
##  [7,]  33049155473  29313244873  35689182693  32501213783  3187968910
##  [8,]  35933757402  31851840065  39004833535  35428336800  3576496735
##  [9,]  38749864977  34526921978  42246475461  38386698720  3859776742
## [10,]  42964345101  37674280738  46339698388  42006989563  4332708825
## [11,]  43702318276  40337571439  48817407207  44577489323  4239917884
## [12,]  47632993913  43262380567  51886564339  47574472453  4312091886
## [13,]  51627511757  46481792262  55567364282  51024578272  4542786010
## [14,]  66610917871  52393435454  65942716502  59168075978  6774640524
## [15,]  81658052383  61882368981  83637118311  72759743646 10877374665
## [16,]  86883552689  71695008675  98858723339  85276866007 13581857332
## [17,]  88290755581  80860819631 109166642522  95013731077 14152911446
## [18,]  98868135190  88925123961 115093711258 102009417610 13084293649
## [19,] 110802000000  96211110865 119787301029 107999205947 11788095082
## [20,] 129781000000 106935472693 134340154504 120637813598 13702340905
## [21,] 143250000000 120675283798 155652355735 138163819766 17488535969
## [22,] 154793000000 134656500000 174730316322 154693408161 20036908161
## [23,] 166208000000 148508000000 190136371755 169322185877 20814185877
## [24,] 168165000000 158104000000 193340108101 175722054051 17618054051
## [25,] 169303000000 164617250000 190908875000 177763062500 13145812500
## [26,] 177750000000 170356500000 190276625000 180316562500  9960062500
## [27,] 179898000000 173779000000 187908625000 180843812500  7064812500
## [28,] 185812000000 178190750000 191100500000 184645625000  6454875000
## [29,] 202093000000 186388250000 204807500000 195597875000  9209625000
## [30,] 214473000000 195569000000 219743500000 207656250000 12087250000
## [31,] 239433000000 210452750000 246057875000 228255312500 17802562500
## [32,] 240271000000 224067500000 263963750000 244015625000 19948125000
## [33,] 243565000000 234435500000 271044125000 252739812500 18304312500
## [34,] 251212000000 243620250000 274572750000 259096500000 15476250000
## [35,] 260992000000 249010000000 271463375000 260236687500 11226687500
## [36,] 265429000000 255299500000 274715875000 265007687500  9708187500
## [37,] 267815000000 261362000000 279440125000 270401062500  9039062500
## [38,] 268834000000 265767500000 281583000000 273675250000  7907750000
## [39,] 272519000000 268649250000 280408625000 274528937500  5879687500
## [40,] 273433000000 270650250000 278736250000 274693250000  4043000000
## [41,]           NA           NA           NA           NA          NA
## [42,]           NA           NA           NA           NA          NA
## [43,]           NA           NA           NA           NA          NA
## [44,]           NA           NA           NA           NA          NA
## [45,]           NA           NA           NA           NA          NA
## [46,]           NA           NA           NA           NA          NA
## [47,]           NA           NA           NA           NA          NA
##            ramalan
##  [1,]           NA
##  [2,]           NA
##  [3,]           NA
##  [4,]           NA
##  [5,]           NA
##  [6,]           NA
##  [7,]           NA
##  [8,]  35689182693
##  [9,]  39004833535
## [10,]  42246475461
## [11,]  46339698388
## [12,]  48817407207
## [13,]  51886564339
## [14,]  55567364282
## [15,]  65942716502
## [16,]  83637118311
## [17,]  98858723339
## [18,] 109166642522
## [19,] 115093711258
## [20,] 119787301029
## [21,] 134340154504
## [22,] 155652355735
## [23,] 174730316322
## [24,] 190136371755
## [25,] 193340108101
## [26,] 190908875000
## [27,] 190276625000
## [28,] 187908625000
## [29,] 191100500000
## [30,] 204807500000
## [31,] 219743500000
## [32,] 246057875000
## [33,] 263963750000
## [34,] 271044125000
## [35,] 274572750000
## [36,] 271463375000
## [37,] 274715875000
## [38,] 279440125000
## [39,] 281583000000
## [40,] 280408625000
## [41,] 278736250000
## [42,] 282779250000
## [43,] 286822250000
## [44,] 290865250000
## [45,] 294908250000
## [46,]           NA
## [47,]           NA

#Plot time series
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="GDP", 
        main= "DMA N=3 Data GDP", xlim=c(1, 45))
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)

Menghitung nilai keakuratan

#Menghitung nilai keakuratan
error.dma = data.ts-dt.ramal[1:length(data.ts)]
SSE.dma = sum(error.dma[8:length(data.ts)]^2) #8 diambil dari n=4, nx2
MSE.dma = mean(error.dma[8:length(data.ts)]^2)
MAPE.dma = mean(abs((error.dma[8:length(data.ts)]/data.ts[8: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   4.179738e+21
## MSE   1.266587e+20
## MAPE  5.931604e+00

Metode DMA

#membagi training dan testing
training<-dtafrica[1:32,2]
testing<-dtafrica[33:40,2]

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

Plot DES

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

Menghitung nilai keakuratan

#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.2366748
##  gamma: FALSE
## 
## Coefficients:
##           [,1]
## a 240271000000
## b  10230378892

plot(des.opt)
legend("topleft", c("Data Aktual", "Peramalan"), col = c("black", "red"), 
       lty = c(1,1))

#ramalan
ramalandesopt<- forecast(des.opt, h=8)
ramalandesopt

##    Point Forecast        Lo 80        Hi 80        Lo 95        Hi 95
## 33   250501378892 242608988136 258393769647 238431008843 262571748941
## 34   260731757783 248179718568 273283796999 241535070380 279928445187
## 35   270962136675 253851616721 288072656629 244793854596 297130418754
## 36   281192515567 259399930385 302985100748 247863632645 314521398488
## 37   291422894458 264758824269 318086964648 250643718139 332202070778
## 38   301653273350 269905583247 333400963453 253099371293 350207175407
## 39   311883652242 274833353249 348933951235 255220109709 368547194775
## 40   322114031134 279542004543 364686057724 257005734980 387222327288

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 
## [1]        NA        NA 304491678 914329969 891187294 735477789

mapedes.train <- sum(abs(sisaandes[4:length(training.ts)]/training.ts[4: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.149585e+21
## MSE                     3.592454e+19
## MAPE                    3.820348e+00

#Akurasi data testing
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt

## Time Series:
## Start = 33 
## End = 40 
## Frequency = 1 
## [1]  6936378892  9519757783  9970136675 15763515567 23607894458 32819273350
## [7] 39364652242 48681031134

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                     7.550609e+20
## MSE                     7.550609e+20
## MAPE                    8.708468e+00

Perbandingan Metode DMA vs 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   4.179738e+21                    1.149585e+21
## MSE   1.266587e+20                    3.592454e+19
## MAPE  5.931604e+00                    3.820348e+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.

Regresi Time Series

Setelah melakukan peramalan, data yang telah dimasukkan kemudian dieksplorasi. Eksplorasi pertama yang dilakukan adalah dengan menggunakan scatter plot.

Eksplorasi Data Regresi

#Pembuatan Scatter Plot
plot(dtafrica$Tanggal,dtafrica$GDP, pch = 20, col = "blue",
     main = "Scatter Plot Tahun vs Nilai GDP",
     xlab = "Tahun",
     ylab = "Nilai GDP")

#Menampilkan Nilai Korelasi
cor(dtafrica$Tanggal,dtafrica$GDP)

## [1] 0.9868792

Berdasarkan 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. Hal tersebut juga diperkuat dengan hasil perhitungan aplikasi R di mana didapatkan nilai korelasi sebesar \(0.9868792\).

Pembuatan Model Regresi

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

#model regresi
modelaf<- lm(GDP~Tanggal, data = dtafrica)
summary(modelaf)

## 
## Call:
## lm(formula = GDP ~ Tanggal, data = dtafrica)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.719e+10 -1.245e+10 -2.574e+08  1.061e+10  3.442e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.504e+13  4.029e+11  -37.34   <2e-16 ***
## Tanggal      7.668e+09  2.035e+08   37.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.486e+10 on 38 degrees of freedom
## Multiple R-squared:  0.9739, Adjusted R-squared:  0.9732 
## F-statistic:  1420 on 1 and 38 DF,  p-value: < 2.2e-16

Model yang dihasilkan adalah \[y_i=-1.504e13+7.668e09x_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.9732\). Artinya, sebesar 97.32% 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.

Uji Sisaan

sisaan<- residuals(modelaf)
fitValue<- predict(modelaf)

# Q-Q Plot for residuals
qqnorm(sisaan)
qqline(sisaan, col = "steelblue", lwd = 2)

# Sisaan vs Fitted Values Plot
plot(fitValue, sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Fitted values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)

# Histogram of Sisaan
hist(sisaan, col = "steelblue", main = "Histogram of Sisaan", xlab = "Sisaan")

# Sisaan vs Order Plot
plot(seq_along(sisaan), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq_along(sisaan), 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, namun histogram dari sisaan kurang menunjukkan sisaan menyebar normal 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.

Uji Formal

Melihat Sisaan Menyebar Normal/Tidak

Hipotesis H0: sisaan mengikuti sebaran normal H1: sisaan tidak mengikuti sebaran normal

shapiro.test(sisaan)

## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan
## W = 0.98547, p-value = 0.8785

ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.076687, p-value = 0.9584
## alternative hypothesis: two-sided

Berdasarkan uji formal Saphiro-Wilk dan Kolmogorov-Smirnov didapatkan nilai p-value > \(\alpha\) (5%). Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal.

Pendeteksian Autokorelasi

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

Berdasarkan plot ACF dan PACF, terlihat semua diluar rentang batas dan terdapat lag yang signifikan. Namun, untuk lebih memastikan akan dilakukan uji formal dengan uji Durbin Watson.

Hipotesis H0: tidak ada autokorelasi H1: ada autokorelasi

#Deteksi autokorelasi dengan uji-Durbin Watson
dwtest(modelaf)

## 
##  Durbin-Watson test
## 
## data:  modelaf
## DW = 0.16924, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

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(modelaf)
modelCO

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = GDP ~ Tanggal, data = dtafrica)
## 
##  number of interaction: 17
##  rho 0.850065
## 
## Durbin-Watson statistic 
## (original):    0.16924 , p-value: 2.284e-19
## (transformed): 1.19653 , p-value: 2.384e-03
##  
##  coefficients: 
##   (Intercept)       Tanggal 
## -1.632997e+13  8.312521e+09

Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=-1.632997\text{x}10^{13}+8.312521\text{x}10^{9}x_t\] Untuk nilai \(ρ ̂\) optimum yang digunakan adalah \(0.8500649\). Nilai tersebut dapat diketahui dengan syntax berikut.

#Rho optimum
rho<- modelCO$rho
rho

## [1] 0.8500649

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

dtafrica$GDP

##  [1]  21216962290  22307471356  23702472100  25779376633  28049537325
##  [6]  30374910060  33049155473  35933757402  38749864977  42964345101
## [11]  43702318276  47632993913  51627511757  66610917871  81658052383
## [16]  86883552689  88290755581  98868135190 110802000000 129781000000
## [21] 143250000000 154793000000 166208000000 168165000000 169303000000
## [26] 177750000000 179898000000 185812000000 202093000000 214473000000
## [31] 239433000000 240271000000 243565000000 251212000000 260992000000
## [36] 265429000000 267815000000 268834000000 272519000000 273433000000

dtafrica$GDP[-1]

##  [1]  22307471356  23702472100  25779376633  28049537325  30374910060
##  [6]  33049155473  35933757402  38749864977  42964345101  43702318276
## [11]  47632993913  51627511757  66610917871  81658052383  86883552689
## [16]  88290755581  98868135190 110802000000 129781000000 143250000000
## [21] 154793000000 166208000000 168165000000 169303000000 177750000000
## [26] 179898000000 185812000000 202093000000 214473000000 239433000000
## [31] 240271000000 243565000000 251212000000 260992000000 265429000000
## [36] 267815000000 268834000000 272519000000 273433000000

#Transformasi Manual
gdp.trans<- dtafrica$GDP[-1]-dtafrica$GDP[-40]*rho
tahun.trans<- dtafrica$Tanggal[-1]-dtafrica$Tanggal[-40]*rho
modelCOmanual<- lm(gdp.trans~tahun.trans)
summary(modelCOmanual)

## 
## Call:
## lm(formula = gdp.trans ~ tahun.trans)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -8.287e+09 -3.746e+09 -1.507e+09  3.243e+09  1.827e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.448e+12  1.568e+11  -15.62   <2e-16 ***
## tahun.trans  8.313e+09  5.267e+08   15.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.55e+09 on 37 degrees of freedom
## Multiple R-squared:  0.8707, Adjusted R-squared:  0.8672 
## F-statistic: 249.1 on 1 and 37 DF,  p-value: < 2.2e-16

#Mencari Penduga Koefisien Regresi setelah Transformasi ke Persamaan Awal
b0bintang <- modelCOmanual$coefficients[-2]
b0 <- b0bintang/(1-rho)
b1 <- modelCOmanual$coefficients[-1]
b0

##   (Intercept) 
## -1.632997e+13

b1

## tahun.trans 
##  8312520660

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.8,1, by= 0.01))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, modelaf))}))
round(tab, 4)

##     rho          SSE
## 1  0.80 1.160293e+21
## 2  0.81 1.152886e+21
## 3  0.82 1.147122e+21
## 4  0.83 1.143002e+21
## 5  0.84 1.140526e+21
## 6  0.85 1.139693e+21
## 7  0.86 1.140504e+21
## 8  0.87 1.142959e+21
## 9  0.88 1.147058e+21
## 10 0.89 1.152800e+21
## 11 0.90 1.160186e+21
## 12 0.91 1.169216e+21
## 13 0.92 1.179890e+21
## 14 0.93 1.192207e+21
## 15 0.94 1.206168e+21
## 16 0.95 1.221773e+21
## 17 0.96 1.239021e+21
## 18 0.97 1.257913e+21
## 19 0.98 1.278449e+21
## 20 0.99 1.300629e+21
## 21 1.00 1.363601e+21

#Rho optimal di sekitar 0.9
rOpt <- seq(0.9, 1, by= 0.01)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, modelaf))}))
head(tabOpt[order(tabOpt$SSE),])

##    rho          SSE
## 1 0.90 1.160186e+21
## 2 0.91 1.169216e+21
## 3 0.92 1.179890e+21
## 4 0.93 1.192207e+21
## 5 0.94 1.206168e+21
## 6 0.95 1.221773e+21

#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.9, y=1.160186e+21, labels = "rho=0.9", cex = 0.8)

Model Terbaik

#Model terbaik
modelHL <- hildreth.lu.func(0.9, modelaf)
summary(modelHL)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -7.940e+09 -3.929e+09 -1.871e+09  3.447e+09  1.805e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.693e+12  1.585e+11  -10.68 7.39e-13 ***
## x            8.609e+09  7.967e+08   10.81 5.32e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.6e+09 on 37 degrees of freedom
## Multiple R-squared:  0.7594, Adjusted R-squared:  0.7529 
## F-statistic: 116.8 on 1 and 37 DF,  p-value: 5.323e-13

#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.9), "+", coef(modelHL)[2],"x", sep = "")

## y = -1.692582e+13+8609002215x

#Deteksi autokorelasi
dwtest(modelHL)

## 
##  Durbin-Watson test
## 
## data:  modelHL
## DW = 1.2275, p-value = 0.003355
## alternative hypothesis: true autocorrelation is greater than 0

#Deteksi autokorelasi
dwtest(modelHL)

## 
##  Durbin-Watson test
## 
## data:  modelHL
## DW = 1.2275, p-value = 0.003355
## alternative hypothesis: true autocorrelation is greater than 0

#Perbandingan
sseModelawal <- anova(modelaf)$`Sum Sq`[-1]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-1]
sseModelHL <- anova(modelHL)$`Sum Sq`[-1]
mseModelawal <- sseModelawal/length(dtafrica$GDP)
mseModelCO <- sseModelCO/length(dtafrica$GDP)
mseModelHL <- sseModelHL/length(dtafrica$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 8.389669e+21          1.139693e+21      1.160186e+21
## MSE 2.097417e+20          2.849233e+19      2.900466e+19

Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu telah berhasil terbukti dapat mengurangi nilai SSE dan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar \(8.389669\text{x}10^{21}\) menjadi \(1.139693\text{x}10^{21}\) pada model Cochrane-Orcutt dan sebesar \(1.160186\text{x}10^{21}\) pada model Hildreth-Lu.

Akan tetapi penanganan autokorelasi masih terjadi, sehingga diperlukan metode yang lebih lanjut

Dengan model terbaik yang terbangun \[y = -1.692582\text{x}10^{13}+8609002215x\]