Pertemuan 2 - Regresi

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)

## Warning: package 'lmtest' was built under R version 4.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(orcutt)

## Warning: package 'orcutt' was built under R version 4.3.3

library(HoRM)

## Warning: package 'HoRM' was built under R version 4.3.3

library(readxl)

## Warning: package 'readxl' was built under R version 4.3.2

Input Data

# input data excel
data_gdp <- read_excel("C:/Users/user/Downloads/DATAMPDWPERT2.xlsx")
data_gdp

## # A tibble: 44 × 2
##    Tahun          GDP
##    <dbl>        <dbl>
##  1  1960 33090214944.
##  2  1961 34977868084.
##  3  1962 37290528805.
##  4  1963 40204671362.
##  5  1964 43801439277.
##  6  1965 47294334234.
##  7  1966 51087667657.
##  8  1967 52466543012.
##  9  1968 55815747943.
## 10  1969 62314473363.
## # ℹ 34 more rows

Eksplorasi Data

Sebelum melakukan regresi, akan diperlihatkan plot time-series dari GDP Sub SaharaPeriode 1960-2023

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

## Time Series:
## Start = 1 
## End = 44 
## Frequency = 1 
##  [1]  33090214944  34977868084  37290528805  40204671362  43801439277
##  [6]  47294334234  51087667657  52466543012  55815747943  62314473363
## [11]  70436667211  71791051606  80983186014 103240920672 131389895028
## [16] 144581341759 157959271369 171680235863 190571667082 227962103836
## [21] 290218140808 402288463246 371770082939 334284647100 299763709596
## [26] 284674940036 265263726339 298330151730 312818799588 316929212105
## [31] 373310030640 401440135787 362732355954 368788820474 377141069077
## [36] 479946737826 535109490236 564078278694 565282335759 403403922298
## [41] 427688179171 409106816914 445363669975 561066946684

#Membuat plot time series
ts.plot(data.ts, xlab="Time Period ", ylab="IPM", 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)
}

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,]  33090214944           NA           NA           NA           NA
##  [2,]  34977868084           NA           NA           NA           NA
##  [3,]  37290528805  35119537278           NA           NA           NA
##  [4,]  40204671362  37491022751           NA           NA           NA
##  [5,]  43801439277  40432213148  45934790661  43183501904   2751288756
##  [6,]  47294334234  43766814958  50173744302  46970279630   3203464672
##  [7,]  51087667657  47394480390  54454435505  50924457947   3529977558
##  [8,]  52466543012  50282848301  56552449138  53417648719   3134800418
##  [9,]  55815747943  53123319538  58836193127  55979756332   2856436795
## [10,]  62314473363  56865588106  63748927022  60307257564   3441669458
## [11,]  70436667211  62855629506  73337197084  68096413295   5240783789
## [12,]  71791051606  68180730727  79274226621  73727478674   5546747947
## [13,]  80983186014  74403634944  86250908047  80327271495   5923636552
## [14,] 103240920672  85338386097 104066657114  94702521606   9364135508
## [15,] 131389895028 105204667238 138982876195 122093771716  16889104478
## [16,] 144581341759 126404052486 167914086911 147159069699  20755017212
## [17,] 157959271369 144643502719 183095693194 163869597956  19226095238
## [18,] 171680235863 158073616330 188140067967 173106842149  15033225819
## [19,] 190571667082 173403724771 202797278434 188100501603  14696776831
## [20,] 227962103836 196738002260 238070444539 217404223400  20666221140
## [21,] 290218140808 236250637242 304490335543 270370486393  34119849151
## [22,] 402288463246 306822902630 427261013135 367041957883  60219055253
## [23,] 371770082939 354758895664 465721729969 410240312816  55481417152
## [24,] 334284647100 369447731095 420990173692 395218952394  25771221299
## [25,] 299763709596 335272813211 299498812987 317385813099 -17887000112
## [26,] 284674940036 306241098911 244748867921 275494983416 -30746115495
## [27,] 265263726339 283234125323 233203684340 258218904832 -25015220492
## [28,] 298330151730 282756272702 266781153481 274768713091  -7987559610
## [29,] 312818799588 292137559219 304327372827 298232466023   6094906804
## [30,] 316929212105 309359387808 338576016938 323967702373  14608314565
## [31,] 373310030640 334352680778 379158290463 356755485620  22402804843
## [32,] 401440135787 363893126177 419942582023 391917854100  28024727923
## [33,] 362732355954 379160840793 419211423882 399186132337  20025291544
## [34,] 368788820474 377653770738 385822820409 381738295574   4084524835
## [35,] 377141069077 369554081835 357749783260 363651932548  -5902149287
## [36,] 479946737826 408625542459 455321030689 431973286574  23347744115
## [37,] 535109490236 464065765713 564033703801 514049734757  49983969044
## [38,] 564078278694 526378168919 646421522029 586399845474  60021676555
## [39,] 565282335759 554823368230 634291902781 594557635505  39734267276
## [40,] 403403922298 510921512250 471349170485 491135341367 -19786170883
## [41,] 427688179171 465458145742 375572419746 420515282744 -44942862998
## [42,] 409106816914 413399639461 313679386747 363539513104 -49860126357
## [43,] 445363669975 427386222020 411329327911 419357774965  -8028447055
## [44,] 561066946684 471845811191 540449651792 506147731491  34301920300
## [45,]           NA           NA           NA           NA           NA
## [46,]           NA           NA           NA           NA           NA
## [47,]           NA           NA           NA           NA           NA
## [48,]           NA           NA           NA           NA           NA
## [49,]           NA           NA           NA           NA           NA
##            ramalan
##  [1,]           NA
##  [2,]           NA
##  [3,]           NA
##  [4,]           NA
##  [5,]           NA
##  [6,]  45934790661
##  [7,]  50173744302
##  [8,]  54454435505
##  [9,]  56552449138
## [10,]  58836193127
## [11,]  63748927022
## [12,]  73337197084
## [13,]  79274226621
## [14,]  86250908047
## [15,] 104066657114
## [16,] 138982876195
## [17,] 167914086911
## [18,] 183095693194
## [19,] 188140067967
## [20,] 202797278434
## [21,] 238070444539
## [22,] 304490335543
## [23,] 427261013135
## [24,] 465721729969
## [25,] 420990173692
## [26,] 299498812987
## [27,] 244748867921
## [28,] 233203684340
## [29,] 266781153481
## [30,] 304327372827
## [31,] 338576016938
## [32,] 379158290463
## [33,] 419942582023
## [34,] 419211423882
## [35,] 385822820409
## [36,] 357749783260
## [37,] 455321030689
## [38,] 564033703801
## [39,] 646421522029
## [40,] 634291902781
## [41,] 471349170485
## [42,] 375572419746
## [43,] 313679386747
## [44,] 411329327911
## [45,] 540449651792
## [46,] 574751572092
## [47,] 609053492393
## [48,] 643355412693
## [49,] 677657332994

#Plot time series
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="GDP", 
        main= "DMA N=3 Data GDP")
points(dt.gab[,1])
points(dt.gab[,3])
points(dt.gab[,6])
lines(dt.gab[,3],col="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   1.880773e+23
## MSE   4.822494e+21
## MAPE  1.331188e+01

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<-data_gdp[1:36,2]
testing<-data_gdp[37:45,2]

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

#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.0506884
##  gamma: FALSE
## 
## Coefficients:
##           [,1]
## a 479946737826
## b  13388718828

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
## 37   493335456654 450451987784 536218925524 427750850637 558920062671
## 38   506724175482 444521765024 568926585940 411593797956 601854553008
## 39   520112894310 442011077504 598214711117 400666473126 639559315495
## 40   533501613139 441085686145 625917540133 392163649551 674839576726
## 41   546890331967 441055464919 652725199014 385029870556 708750793378

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,]  425007581
## [4,] 1004946463
## [5,] 1636632692
## [6,] 1449801441

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                     3.846480e+22
## MSE                     1.068467e+21
## MAPE                    7.885089e+00

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

## Time Series:
## Start = 37 
## End = 41 
## Frequency = 1 
##        testing.ts
## [1,] -41774033582
## [2,] -57354103211
## [3,] -45169441448
## [4,] 130097690841
## [5,] 119202152796

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                     4.245489e+21
## MSE                     4.245489e+21
## MAPE                    9.565138e+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   1.880773e+23                    3.846480e+22
## MSE   4.822494e+21                    1.068467e+21
## MAPE  1.331188e+01                    7.885089e+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(data_gdp$Tahun,data_gdp$GDP, pch = 20, col = "blue",
     main = "Scatter Plot Tahun vs Nilai GDP",
     xlab = "Tahun",
     ylab = "Nilai GDP")

#Menampilkan Nilai Korelasi
cor(data_gdp$Tahun,data_gdp$GDP)

## [1] 0.947123

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

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

Regresi

#Pembuatan Model Regresi
#model regresi
model<- lm(data_gdp$GDP~data_gdp$Tahun, data = data_gdp)
summary(model)

## 
## Call:
## lm(formula = data_gdp$GDP ~ data_gdp$Tahun, data = data_gdp)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -9.332e+10 -3.423e+10 -2.060e+10  2.927e+10  1.527e+11 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -2.479e+13  1.310e+12  -18.93   <2e-16 ***
## data_gdp$Tahun  1.264e+10  6.609e+08   19.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.567e+10 on 42 degrees of freedom
## Multiple R-squared:  0.897,  Adjusted R-squared:  0.8946 
## F-statistic: 365.9 on 1 and 42 DF,  p-value: < 2.2e-16

Model yang dihasilkan adalah \[y_i=5.567e+010x_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.8946\). Artinya, sebesar 89.46% keragaman nilai IPM 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,44,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,44,1), sisaan, col = "red")
abline(a = 0, b = 0, lwd = 2)

Dua plot di samping kiri digunakan untuk melihat apakah sisaan menyebar normal. Normal Q-Q Plot di atas menunjukkan bahwa sisaan cenderung menyebar normal, tetapi histogram dari sisaan tidak menunjukkan demikian. Selanjutnya, dua plot di samping kanan digunakan untuk melihat autokorelasi. Plot Sisaan vs Fitted Value dan Plot Sisaan vs Order menunjukkan adanya pola pada sisaan. Untuk lebih lanjut akan digunakan uji formal melihat normalitas sisaan dan plot ACF dan PACF untuk melihat apakah ada autokorelasi atau tidak.

#Melihat Sisaan Menyebar Normal/Tidak
#H0: sisaan mengikuti sebaran normal
#H1: sisaan tidak mengikuti sebaran normal
shapiro.test(sisaan)

## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan
## W = 0.92908, p-value = 0.009703

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

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

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

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

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

#Deteksi autokorelasi dengan uji-Durbin Watson
#H0: tidak ada autokorelasi
#H1: ada autokorelasi
dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.62333, p-value = 1.501e-08
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil DW Test, didapatkan nilai \(DW = 0.62333\) dan p-value = $1.501e-08. Berdasarkan tabel Durbin-Watson diperoleh nilai $DL = 1.4754 dan $DU = 1.5660. Nilai DW masih berada di bawah nilai DL dan DU. Artinya, berada di daerah konklusif, dapat dikatakan berada di daerah autokorelasi positif. Namun, 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 = data_gdp$GDP ~ data_gdp$Tahun, data = data_gdp)
## 
##  number of interaction: 5
##  rho 0.681854
## 
## Durbin-Watson statistic 
## (original):    0.62333 , p-value: 1.501e-08
## (transformed): 1.39649 , p-value: 1.361e-02
##  
##  coefficients: 
##    (Intercept) data_gdp$Tahun 
##  -2.628332e+13   1.339153e+10

#Rho optimum
rho<- modelCO$rho
rho

## [1] 0.6818543

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

data_gdp$GDP

##  [1]  33090214944  34977868084  37290528805  40204671362  43801439277
##  [6]  47294334234  51087667657  52466543012  55815747943  62314473363
## [11]  70436667211  71791051606  80983186014 103240920672 131389895028
## [16] 144581341759 157959271369 171680235863 190571667082 227962103836
## [21] 290218140808 402288463246 371770082939 334284647100 299763709596
## [26] 284674940036 265263726339 298330151730 312818799588 316929212105
## [31] 373310030640 401440135787 362732355954 368788820474 377141069077
## [36] 479946737826 535109490236 564078278694 565282335759 403403922298
## [41] 427688179171 409106816914 445363669975 561066946684

data_gdp$GDP[-1]

##  [1]  34977868084  37290528805  40204671362  43801439277  47294334234
##  [6]  51087667657  52466543012  55815747943  62314473363  70436667211
## [11]  71791051606  80983186014 103240920672 131389895028 144581341759
## [16] 157959271369 171680235863 190571667082 227962103836 290218140808
## [21] 402288463246 371770082939 334284647100 299763709596 284674940036
## [26] 265263726339 298330151730 312818799588 316929212105 373310030640
## [31] 401440135787 362732355954 368788820474 377141069077 479946737826
## [36] 535109490236 564078278694 565282335759 403403922298 427688179171
## [41] 409106816914 445363669975 561066946684

#Transformasi Manual
gdp.trans<- data_gdp$GDP[-1]-data_gdp$GDP[-12]*rho
tahun.trans<- data_gdp$Tahun[-1]-data_gdp$Tahun[-12]*rho
modelCOmanual<- lm(gdp.trans~tahun.trans)
summary(modelCOmanual)

## 
## Call:
## lm(formula = gdp.trans ~ tahun.trans)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -3.178e+10 -1.085e+10 -5.526e+09  7.954e+09  5.088e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.325e+12  4.627e+11  -17.99   <2e-16 ***
## tahun.trans  1.333e+10  7.336e+08   18.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.794e+10 on 41 degrees of freedom
## Multiple R-squared:  0.8896, Adjusted R-squared:  0.8869 
## F-statistic: 330.2 on 1 and 41 DF,  p-value: < 2.2e-16

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) 
## -2.616685e+13

b1

## tahun.trans 
## 13331259904

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(-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  -1.0 4.322857e+23
## 2  -0.9 3.902052e+23
## 3  -0.8 3.507034e+23
## 4  -0.7 3.137803e+23
## 5  -0.6 2.794359e+23
## 6  -0.5 2.476702e+23
## 7  -0.4 2.184832e+23
## 8  -0.3 1.918748e+23
## 9  -0.2 1.678452e+23
## 10 -0.1 1.463942e+23
## 11  0.0 1.275220e+23
## 12  0.1 1.112284e+23
## 13  0.2 9.751353e+22
## 14  0.3 8.637734e+22
## 15  0.4 7.781985e+22
## 16  0.5 7.184104e+22
## 17  0.6 6.844093e+22
## 18  0.7 6.761950e+22
## 19  0.8 6.937677e+22
## 20  0.9 7.371273e+22
## 21  1.0 8.112305e+22

Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.1. Namun, hasil tersebut masih kurang teliti sehingga akan dicari kembali \(ρ\) yang lebih optimum dengan ketelitian yang lebih. Jika sebelumnya jarak antar \(ρ\) yang dicari adalah 0.1, kali ini jarak antar \(ρ\) adalah 0.001 dan dilakukan pada selang 0.0 sampai dengan 0.1.

#Rho optimal di sekitar 0.4
rOpt <- seq(0.0,0.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
## 101 0.100 1.112284e+23
## 100 0.099 1.113786e+23
## 99  0.098 1.115290e+23
## 98  0.097 1.116797e+23
## 97  0.096 1.118306e+23
## 96  0.095 1.119818e+23

#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.341, y=0.2397500, labels = "rho=0.341", cex = 0.8)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai \(ρ\) optimum, yaitu saat SSE terkecil terdapat pada nilai \(ρ=1.112284e+23\). 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.341, model)
summary(modelHL)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.106e+11 -2.483e+10 -1.180e+10  2.093e+10  1.356e+11 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.674e+13  1.093e+12  -15.31   <2e-16 ***
## x            1.295e+10  8.368e+08   15.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.487e+10 on 41 degrees of freedom
## Multiple R-squared:  0.8538, Adjusted R-squared:  0.8502 
## F-statistic: 239.4 on 1 and 41 DF,  p-value: < 2.2e-16

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

## y = -2.53997e+13+12946738938x

Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=-2.53997e+13+12946738938x_t\]

#Deteksi autokorelasi
dwtest(modelHL)

## 
##  Durbin-Watson test
## 
## data:  modelHL
## DW = 0.963, p-value = 5.394e-05
## alternative hypothesis: true autocorrelation is greater than 0

Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(0.963\) berada pada selang daerah tidak ada autokorelasi. Hal tersebut juga didukung oleh p-value sebesar \(5.394e-05\), di mana p-value >\(\alpha\)=5%. Artinya tak tolak \(H_0\) atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai GDP 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(data_gdp$GDP)
mseModelCO <- sseModelCO/length(data_gdp$GDP)
mseModelHL <- sseModelHL/length(data_gdp$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.301541e+23          1.318888e+22      8.255688e+22
## MSE 2.958048e+21          2.997472e+20      1.876293e+21

Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki SSE sebesar $ 1.318888e+22 $ dan $ 8.255688e+22$ lebih buruk dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar $1.301541e+23 $.

Simpulan

Autokorelasi yang terdapat pada data GDP 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. Namun, autokorelasi tersebut tidak dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode menghasilkan nilai SSE yang sama, artinya keduanya tidak baik untuk digunakan.

Daftar Pustaka

Aprianto A, Debataraja NN, Imro’ah N. 2020. Metode cochrane-orcutt untuk mengatasi autokorelasi pada estimasi parameter ordinary least squares. Bimaster : Buletin Ilmiah Matematika, Statistika dan Terapannya. 9(1):95–102. doi:10.26418/bbimst.v9i1.38590.

BPS. 2021a. Indeks Pembangunan Manusia 2020. Jakarta (ID): Badan Pusat Statistik.

BPS. 2021b. Indeks Pembangunan Manusia (IPM) 2021. Berita Resmi Statistik., siap terbit.