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, union
library(TTR)
## Warning: package 'TTR' was built under R version 4.3.3
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(lmtest) #digunakan untuk uji formal pendeteksian autokorelasi
## Warning: package 'lmtest' was built under R version 4.3.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.2
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(orcutt) #untuk membuat model regresi Cochrane-Orcutt
## Warning: package 'orcutt' was built under R version 4.3.3
library(HoRM) #untuk membuat model regresi Hildreth-Lu
## Warning: package 'HoRM' was built under R version 4.3.3

Input Data

Data yang digunakan adalah data GDP Malta dalam Milyar USD pada periode tahun 2000-2023.

library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
dtMalta <-read_xlsx("C:/Users/acer/OneDrive - apps.ipb.ac.id/Semester 5/MPDW/Pertemuan2/GDP.xlsx", sheet="Malta")
dtMalta$GDP <- (dtMalta$GDP)/1000000000 #standarisasi dalam bentuk Milyar USD
dtMalta
## # A tibble: 24 × 3
##    Periode tahun   GDP
##      <dbl> <dbl> <dbl>
##  1       1  2000  4.07
##  2       2  2001  4.09
##  3       3  2002  4.47
##  4       4  2003  5.45
##  5       5  2004  6.10
##  6       6  2005  6.42
##  7       7  2006  6.78
##  8       8  2007  7.93
##  9       9  2008  9.09
## 10      10  2009  8.70
## # ℹ 14 more rows

Eksplorasi Data

Sebelum melakukan regresi, akan diperlihatkan plot time-series dari GDP Malta dalam Milyar USD pada periode tahun 2000-2023.

#Membentuk objek time series
data.ts<-ts(dtMalta$GDP)
data.ts
## Time Series:
## Start = 1 
## End = 24 
## Frequency = 1 
##  [1]  4.069516  4.088369  4.470446  5.448416  6.098093  6.416184  6.778324
##  [8]  7.925371  9.090407  8.696367  9.035824  9.638729  9.461776 10.551031
## [15] 11.625849 11.091484 11.668016 13.484542 15.404393 16.004564 15.249999
## [22] 18.123719 18.357069 20.956999
#Membuat plot time series
ts.plot(data.ts, xlab="Time Period ", ylab="GDP", main= "Time Series Plot of GDP Malta (Milyar USD)")
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)
}
dt.gab <- cbind(aktual = c(data.ts,rep(NA,5)), 
                pemulusan1 = c(dt.sma,rep(NA,5)),
                pemulusan2 = c(dt.dma, rep(NA,5)),
                At = c(At, rep(NA,5)), 
                Bt = c(Bt,rep(NA,5)),
                ramalan = c(dt.ramal, f[-1]))
dt.gab
##          aktual pemulusan1 pemulusan2        At        Bt   ramalan
##  [1,]  4.069516         NA         NA        NA        NA        NA
##  [2,]  4.088369         NA         NA        NA        NA        NA
##  [3,]  4.470446   4.209444         NA        NA        NA        NA
##  [4,]  5.448416   4.669077         NA        NA        NA        NA
##  [5,]  6.098093   5.338985   6.538617  5.938801 0.5998163        NA
##  [6,]  6.416184   5.987564   7.298942  6.643253 0.6556888  6.538617
##  [7,]  6.778324   6.430867   7.454324  6.942595 0.5117283  7.298942
##  [8,]  7.925371   7.039960   8.147619  7.593790 0.5538296  7.454324
##  [9,]  9.090407   7.931367   9.525973  8.728670 0.7973027  8.147619
## [10,]  8.696367   8.570715  10.017450  9.294083 0.7233676  9.525973
## [11,]  9.035824   8.940866   9.860632  9.400749 0.4598832 10.017450
## [12,]  9.638729   9.123640   9.614106  9.368873 0.2452330  9.860632
## [13,]  9.461776   9.378776   9.840808  9.609792 0.2310156  9.614106
## [14,] 10.551031   9.883845  10.727361 10.305603 0.4217580  9.840808
## [15,] 11.625849  10.546219  11.766096 11.156157 0.6099385 10.727361
## [16,] 11.091484  11.089454  12.255351 11.672403 0.5829483 11.766096
## [17,] 11.668016  11.461783  12.320378 11.891080 0.4292975 12.255351
## [18,] 13.484542  12.081347  13.155652 12.618499 0.5371523 12.320378
## [19,] 15.404393  13.518984  15.848875 14.683929 1.1649458 13.155652
## [20,] 16.004564  14.964500  17.850279 16.407389 1.4428896 15.848875
## [21,] 15.249999  15.552985  17.301310 16.427148 0.8741625 17.850279
## [22,] 18.123719  16.459427  18.060340 17.259884 0.8004565 17.301310
## [23,] 18.357069  17.243596  18.893448 18.068522 0.8249263 18.060340
## [24,] 20.956999  19.145929  22.205153 20.675541 1.5296118 18.893448
## [25,]        NA         NA         NA        NA        NA 22.205153
## [26,]        NA         NA         NA        NA        NA 23.734765
## [27,]        NA         NA         NA        NA        NA 25.264376
## [28,]        NA         NA         NA        NA        NA 26.793988
## [29,]        NA         NA         NA        NA        NA 28.323600
#Plot time series
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="GDP (Milyar USD)", 
        main= "DMA N=3 Data GDP Malta (Milyar USD)", ylim=c(0,30))
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      23.453268
## MSE       1.234383
## MAPE      7.002671

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<-dtMalta[1:20,3]
testing<-dtMalta[21:24,3]

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

#eksplorasi data
plot(data.ts, col="red",main="Plot semua data")
points(data.ts)

plot(training.ts, col="blue",main="Plot data training")
points(training.ts)

plot(testing.ts, col="green",main="Plot data testing")
points(testing.ts)

Selanjutnya akan dilakukan pemulusan dengan DES, kali ini langsung dicari lambda dan gamma optimum sebagai berikut. Nilai lambda dan gamma optimum dapat dilihat pada smoothing parameters alpha untuk nilai lambda dan beta untuk nilai gamma.

#Lamda dan gamma optimum
des.opt<- HoltWinters(training.ts, gamma = FALSE)
des.opt
## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = training.ts, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 1
##  beta : 0.1984021
##  gamma: FALSE
## 
## Coefficients:
##        [,1]
## a 16.004564
## b  0.868257
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
## 21       16.87282 15.96754 17.77810 15.48832 18.25733
## 22       17.74108 16.32810 19.15406 15.58011 19.90204
## 23       18.60933 16.71316 20.50551 15.70939 21.50928
## 24       19.47759 17.09413 21.86106 15.83240 23.12279
## 25       20.34585 17.46187 23.22983 15.93518 24.75652

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,]  0.36322399
## [4,]  0.88705138
## [5,]  0.38276630
## [6,] -0.02476059
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                        9.5011138
## MSE                        0.4750557
## MAPE                       5.7550040
#Akurasi data testing
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt
## Time Series:
## Start = 21 
## End = 24 
## Frequency = 1 
##      testing.ts
## [1,]  1.6228219
## [2,] -0.3826410
## [3,]  0.2522657
## [4,] -1.4794075
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                         1.258062
## MSE                         1.258062
## MAPE                        5.296549

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      23.453268                       9.5011138
## MSE       1.234383                       0.4750557
## MAPE      7.002671                       5.7550040

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(dtMalta$tahun,dtMalta$GDP, pch = 20, col = "blue",
     main = "Scatter Plot Tahun vs Nilai GDP",
     xlab = "Tahun",
     ylab = "Nilai GDP (Milyar USD)")

#Menampilkan Nilai Korelasi
cor(dtMalta$tahun,dtMalta$GDP)
## [1] 0.9753696

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.9753696\).

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

Regresi

#Pembuatan Model Regresi
#model regresi
model<- lm(GDP~tahun, data = dtMalta)
summary(model)
## 
## Call:
## lm(formula = GDP ~ tahun, data = dtMalta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8782 -0.6330  0.1580  0.5043  2.7846 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.319e+03  6.409e+01  -20.57 7.36e-16 ***
## tahun        6.609e-01  3.186e-02   20.74 6.22e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.081 on 22 degrees of freedom
## Multiple R-squared:  0.9513, Adjusted R-squared:  0.9491 
## F-statistic: 430.2 on 1 and 22 DF,  p-value: 6.223e-16

Model yang dihasilkan adalah \[y_i=-1319+0.6609x_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.9513\). Artinya, sebesar 95.13% keragaman nilai GDP Malta 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,24,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,24,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 cenderung menjulur ke kanan. 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.96441, p-value = 0.5332
ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.097548, p-value = 0.9595
## alternative hypothesis: two-sided

Berdasarkan uji formal Shapiro-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, terlihat ada lag yang berada diluar rentang batas, sehingga dapat dikatakan terdapat autokorelasi. 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.73059, p-value = 5.94e-05
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil DW Test, didapatkan nilai \(DW = 0.73059\) dan p-value < \(5.94 \times 10^{-5}\). Berdasarkan tabel Durbin-Watson diperoleh nilai \(DL = 1.273\) dan \(DU = 1.446\). Nilai DW kurang dari nilai DL. Artinya, Ada autokorelasi positif. Dan 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(model)
modelCO
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = GDP ~ tahun, data = dtMalta)
## 
##  number of interaction: 38
##  rho 0.735275
## 
## Durbin-Watson statistic 
## (original):    0.73059 , p-value: 5.94e-05
## (transformed): 2.09325 , p-value: 4.968e-01
##  
##  coefficients: 
##  (Intercept)        tahun 
## -1602.424826     0.801736

Hasil keluaran model setelah dilakukan penanganan Cochrane-Orcutt adalah sebagai berikut. \[y_i=-1602.424826+0.801736x_t\] .Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(2.09325\) dan \(0,4968\). Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.446 < DW < 2.554\). 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.735275\). Nilai tersebut dapat diketahui dengan syntax berikut.

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

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

GDP <- dtMalta$GDP
tahun <- dtMalta$tahun
GDP
##  [1]  4.069516  4.088369  4.470446  5.448416  6.098093  6.416184  6.778324
##  [8]  7.925371  9.090407  8.696367  9.035824  9.638729  9.461776 10.551031
## [15] 11.625849 11.091484 11.668016 13.484542 15.404393 16.004564 15.249999
## [22] 18.123719 18.357069 20.956999
#Transformasi Manual
GDP.trans<- GDP[-1]-GDP[-24]*rho
tahun.trans<- tahun[-1]-tahun[-24]*rho
modelCOmanual<- lm(GDP.trans~tahun.trans, data=dtMalta)
summary(modelCOmanual)
## 
## Call:
## lm(formula = GDP.trans ~ tahun.trans, data = dtMalta)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.62914 -0.42605  0.00462  0.50992  1.71140 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -424.2025    54.1362  -7.836 1.15e-07 ***
## tahun.trans    0.8017     0.1015   7.899 1.01e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8548 on 21 degrees of freedom
## Multiple R-squared:  0.7482, Adjusted R-squared:  0.7362 
## F-statistic: 62.39 on 1 and 21 DF,  p-value: 1.009e-07

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) 
##   -1602.425
b1
## tahun.trans 
##   0.8017358

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 21.9992
## 2  0.2 20.0687
## 3  0.3 18.4680
## 4  0.4 17.1971
## 5  0.5 16.2562
## 6  0.6 15.6450
## 7  0.7 15.3637
## 8  0.8 15.4123
## 9  0.9 15.7907
## 10 1.0 18.6432

Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.7. 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.7
rOpt <- seq(0.6,0.8, 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
## 136 0.735 15.34321
## 137 0.736 15.34322
## 135 0.734 15.34324
## 138 0.737 15.34326
## 134 0.733 15.34330
## 139 0.738 15.34333
#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.7352746, y=15.35, labels = "rho=0.7352746", cex = 0.8)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai \(ρ\) optimum, yaitu saat SSE terkecil terdapat pada nilai \(ρ=0.7352746\). 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.7352746, model)
summary(modelHL)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.62914 -0.42605  0.00462  0.50992  1.71140 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -424.2026    54.1362  -7.836 1.15e-07 ***
## x              0.8017     0.1015   7.899 1.01e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8548 on 21 degrees of freedom
## Multiple R-squared:  0.7482, Adjusted R-squared:  0.7362 
## F-statistic: 62.39 on 1 and 21 DF,  p-value: 1.009e-07
#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.7352746), "+", coef(modelHL)[2],"x", sep = "")
## y = -1602.425+0.8017358x

Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i = -1602.425+0.8017358x_t\]

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

Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(2.0932\) berada pada selang daerah tidak terdapat autokorelasi, yaitu pada rentang DU < DW < 4-DU atau \(1.446 < DW < 2.554\). Hal tersebut juga didukung oleh p-value sebesar \(0.4968\) > \(\alpha\)=5%. Artinya tak tolak \(H_0\) atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai GDP Malta dengan metode Hildreth-Lu pada taraf nyata 5%.

Terakhir, akan dibandingkan nilai SSE dari ketiga metode (metode awal, metode Cochrane-Orcutt, dan Hildreth-Lu).

#Perbandingan
sseModelawal <- anova(model)$`Sum Sq`[-1]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-1]
sseModelHL <- anova(modelHL)$`Sum Sq`[-1]
mseModelawal <- sseModelawal/length(dtMalta$GDP)
mseModelCO <- sseModelCO/length(dtMalta$GDP)
mseModelHL <- sseModelHL/length(dtMalta$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  25.687332            15.3432103        15.3432103
## MSE   1.070305             0.6393004         0.6393004

Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt memiliki SSE dan MSE yang lebih rendah dari Hildreth-Lu dan model awal yaitu \(15.3432103\) dan \(0.6393004\).

Simpulan

Autokorelasi yang terdapat pada data GDP Malta 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, metode Cochrane-Orcutt menghasilkan nilai SSE yang lebih rendah dari metode Hildreth-Lu, sehingga dalam kasus ini, metode Cochrane-Orcutt lebih baik untuk digunakan.