Pemanggilan Packages

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)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

library(orcutt)
library(HoRM)
library(readxl)

Input Data

Data yang digunakan dalam kesempatan kali ini adalah data GDP Negara Sweden periode tahun 1960-2023.

data <- read_excel("C:/Users/Delita Nur Hasanah/Documents/BISMILLAH SEMESTER 5/MPDW/TUGAS/TUGAS2.xlsx")
data

## # A tibble: 64 × 2
##    Tahun          GDP
##    <dbl>        <dbl>
##  1  1960 15953203853.
##  2  1961 17354781171.
##  3  1962 18821353529.
##  4  1963 20371666154.
##  5  1964 22718426804.
##  6  1965 25000191834.
##  7  1966 27194141951.
##  8  1967 29517675437.
##  9  1968 31323283317.
## 10  1969 34016617942.
## # ℹ 54 more rows

Eksplorasi Data

Sebelum melakukan regresi, akan diperlihatkan plot time-series dari GDP Negara Sweden periode tahun 1960-2023

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

## Time Series:
## Start = 1 
## End = 64 
## Frequency = 1 
##  [1]  15953203853  17354781171  18821353529  20371666154  22718426804
##  [6]  25000191834  27194141951  29517675437  31323283317  34016617942
## [11]  38092452061  41566412923  48954145809  59404966684  66013338965
## [16]  82885397621  89362071673  94468740851 104442351223 123386410161
## [21] 142092068281 129686938223 114380557731 105014356667 109201362581
## [26] 114123537582 150498057724 183009638351 206986674501 217948315625
## [31] 261846194499 274229034312 284321115595 212953336588 229033566615
## [36] 267305875261 291743811512 268146144678 270809066781 274072182417
## [41] 262835454367 242395852494 266849061836 334337212322 385118044877
## [46] 392218088879 423093437424 491252589217 517706149201 436537014294
## [51] 495812558843 574094112973 552483727283 586841821797 581964017237
## [56] 505103781350 515654671470 541018749769 555455371487 533879529188
## [61] 547054174236 639714956069 590409594949 593267701033

#Membuat plot time series
ts.plot(data.ts, xlab="Periode Waktu ", 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)
}

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,]  15953203853           NA           NA           NA           NA
##  [2,]  17354781171           NA           NA           NA           NA
##  [3,]  18821353529  17376446184           NA           NA           NA
##  [4,]  20371666154  18849266951           NA           NA           NA
##  [5,]  22718426804  20637148829  24002871845  22320010337   1682861508
##  [6,]  25000191834  22696761598  26634833208  24665797403   1969035805
##  [7,]  27194141951  24970920197  29376206841  27173563519   2202643322
##  [8,]  29517675437  27237336408  31775330422  29506333415   2268997007
##  [9,]  31323283317  29345033569  33666240591  31505637080   2160603511
## [10,]  34016617942  31619192232  36056535224  33837863728   2218671496
## [11,]  38092452061  34477451106  39804568715  37141009911   2663558804
## [12,]  41566412923  37891827642  44349835605  41120831623   3229003982
## [13,]  48954145809  42871003597  51786155895  47328579746   4457576149
## [14,]  59404966684  49975175139  62766854497  56371014818   6395839679
## [15,]  66013338965  58124150486  73725565310  65924857898   7800707412
## [16,]  82885397621  69434567756  89947774349  79691171053  10256603296
## [17,]  89362071673  79420269419 100274816484  89847542952  10427273532
## [18,]  94468740851  88905403382 108209383107  98557393244   9651989862
## [19,] 104442351223  96091054582 111995345492 104043200037   7952145455
## [20,] 123386410161 107432500745 127344863096 117388681920   9956181175
## [21,] 142092068281 123306943221 152033830632 137670386927  14363443705
## [22,] 129686938223 131721805555 153524583651 142623194603  10901389048
## [23,] 114380557731 128719854745 130327161887 129523508316    803653571
## [24,] 105014356667 116360617540  97880334061 107120475801  -9240141740
## [25,] 109201362581 109532092326  92187900571 100859996449  -8672095878
## [26,] 114123537582 109446418943 104779837623 107113128283  -2333290660
## [27,] 150498057724 124607652629 144765515288 134686583958  10078931329
## [28,] 183009638351 149210411219 192121578462 170665994841  21455583622
## [29,] 206986674501 180164790192 237839134549 209001962370  28837172179
## [30,] 217948315625 202648209492 253262354541 227955282017  25307072525
## [31,] 261846194499 228927061541 278954477141 253940769341  25013707800
## [32,] 274229034312 251341181478 298745909427 275043545453  23702363974
## [33,] 284321115595 273465448135 317907216969 295686332552  22220884417
## [34,] 212953336588 257167828831 250187180864 253677504848  -3490323983
## [35,] 229033566615 242102672933 211150718865 226626695899 -15475977034
## [36,] 267305875261 236430926155 218825159852 227628043003  -8802883152
## [37,] 291743811512 262694417796 293931242133 278312829964  15618412168
## [38,] 268146144678 275731943817 310624306272 293178125045  17446181228
## [39,] 270809066781 276899674324 287148332346 282024003335   5124329011
## [40,] 274072182417 271009131292 263933560920 267471346106  -3537785186
## [41,] 262835454367 269238901188 262951565695 266095233442  -3143667746
## [42,] 242395852494 259767829759 245959581119 252863705439  -6904124320
## [43,] 266849061836 257360122899 247835799466 252597961182  -4762161716
## [44,] 334337212322 281194042217 311367463402 296280752810  15086710592
## [45,] 385118044877 328768106345 408089471394 368428788870  39660682525
## [46,] 392218088879 370557782026 457993392353 414275587189  43717805163
## [47,] 423093437424 400143190393 467450185337 433796687865  33653497472
## [48,] 491252589217 435521371840 502415886013 468968628927  33447257087
## [49,] 517706149201 477350725281 556708650833 517029688057  39678962776
## [50,] 436537014294 481831917571 515693076251 498762496911  16930579340
## [51,] 495812558843 483351907446 488366022140 485858964793   2507057347
## [52,] 574094112973 502147895370 528222539185 515185217278  13037321908
## [53,] 552483727283 540796799700 604859330755 572828065227  32031265528
## [54,] 586841821797 571139887351 637363273772 604251580562  33111693211
## [55,] 581964017237 573763188772 597489649102 585626418937  11863230165
## [56,] 505103781350 557969873461 538660987328 548315430394  -9654443067
## [57,] 515654671470 534240823352 492073212999 513157018176 -21083805176
## [58,] 541018749769 520592400863 486575137471 503583769167 -17008631696
## [59,] 555455371487 537376264242 550655800421 544016032332   6639768090
## [60,] 533879529188 543451216815 562740395832 553095806323   9644589508
## [61,] 547054174236 545463024970 552195404227 548829214599   3366189628
## [62,] 639714956069 573549553165 612339462860 592944508013  19394954848
## [63,] 590409594949 592392908418 636241734219 614317321319  21924412900
## [64,] 593267701033 607797417351 640898999430 624348208390  16550791040
## [65,]           NA           NA           NA           NA           NA
## [66,]           NA           NA           NA           NA           NA
## [67,]           NA           NA           NA           NA           NA
## [68,]           NA           NA           NA           NA           NA
## [69,]           NA           NA           NA           NA           NA
##            ramalan
##  [1,]           NA
##  [2,]           NA
##  [3,]           NA
##  [4,]           NA
##  [5,]           NA
##  [6,]  24002871845
##  [7,]  26634833208
##  [8,]  29376206841
##  [9,]  31775330422
## [10,]  33666240591
## [11,]  36056535224
## [12,]  39804568715
## [13,]  44349835605
## [14,]  51786155895
## [15,]  62766854497
## [16,]  73725565310
## [17,]  89947774349
## [18,] 100274816484
## [19,] 108209383107
## [20,] 111995345492
## [21,] 127344863096
## [22,] 152033830632
## [23,] 153524583651
## [24,] 130327161887
## [25,]  97880334061
## [26,]  92187900571
## [27,] 104779837623
## [28,] 144765515288
## [29,] 192121578462
## [30,] 237839134549
## [31,] 253262354541
## [32,] 278954477141
## [33,] 298745909427
## [34,] 317907216969
## [35,] 250187180864
## [36,] 211150718865
## [37,] 218825159852
## [38,] 293931242133
## [39,] 310624306272
## [40,] 287148332346
## [41,] 263933560920
## [42,] 262951565695
## [43,] 245959581119
## [44,] 247835799466
## [45,] 311367463402
## [46,] 408089471394
## [47,] 457993392353
## [48,] 467450185337
## [49,] 502415886013
## [50,] 556708650833
## [51,] 515693076251
## [52,] 488366022140
## [53,] 528222539185
## [54,] 604859330755
## [55,] 637363273772
## [56,] 597489649102
## [57,] 538660987328
## [58,] 492073212999
## [59,] 486575137471
## [60,] 550655800421
## [61,] 562740395832
## [62,] 552195404227
## [63,] 612339462860
## [64,] 636241734219
## [65,] 640898999430
## [66,] 657449790469
## [67,] 674000581509
## [68,] 690551372548
## [69,] 707102163588

#Plot time series
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="IPM", 
        main= "DMA N=3 Data IPM", ylim=c(10000000000,700000000000))
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   9.962723e+22
## MSE   1.688597e+21
## MAPE  1.056542e+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[1:50,2]
testing<-data[51:64,2]

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

#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.033119
##  gamma: FALSE
## 
## Coefficients:
##           [,1]
## a 436537014294
## b   8429529415

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
## 51   444966543708 410848805523 479084281894 392787965901 497145121516
## 52   453396073123 404340805671 502451340575 378372515708 528419630538
## 53   461825602538 400753913749 522897291327 368424514576 555226690500
## 54   470255131953 398585268463 541924995443 360645534498 579864729408
## 55   478684661368 397263366229 560105956506 354161536168 603207786568

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,]  64995040
## [4,] 146582736
## [5,] 938176087
## [6,] 842109009

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.424901e+22
## MSE                     6.849802e+20
## MAPE                    7.967894e+00

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

## Time Series:
## Start = 51 
## End = 55 
## Frequency = 1 
##         testing.ts
## [1,]  -50846015135
## [2,] -120698039850
## [3,]  -90658124745
## [4,] -116586689844
## [5,] -103279355870

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                     3.545094e+21
## MSE                     3.545094e+21
## MAPE                    6.092990e+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   9.962723e+22                    3.424901e+22
## MSE   1.688597e+21                    6.849802e+20
## MAPE  1.056542e+01                    7.967894e+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$Tahun,data$GDP, pch = 20, col = "blue",
     main = "Scatter Plot Tahun vs Nilai IPM",
     xlab = "Tahun",
     ylab = "Nilai IPM")

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

## [1] 0.9679317

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

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

Regresi

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

## 
## Call:
## lm(formula = GDP ~ Tahun, data = data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.170e+11 -2.328e+10 -5.507e+09  2.047e+10  1.104e+11 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.051e+13  6.845e+11  -29.96   <2e-16 ***
## Tahun        1.043e+10  3.437e+08   30.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.08e+10 on 62 degrees of freedom
## Multiple R-squared:  0.9369, Adjusted R-squared:  0.9359 
## F-statistic: 920.4 on 1 and 62 DF,  p-value: < 2.2e-16

Model yang dihasilkan adalah \[y_i=-2.051e+13 + 1.043e+10x_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.9369\). Artinya, sebesar 93.69% keragaman nilai GDP dapat dijelaskan oleh peubah tahun. Hasil ini menunjukkan hasil yang bagus, seolah mendapatkan hasil terbaik. Namun, kita perlu melakukan uji terhadap sisaannya seperti berikut ini.

#sisaan dan fitted value
sisaan<- residuals(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,64,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,64,1), sisaan, col = "red")
abline(a = 0, b = 0, lwd = 2)

Dua plot di samping kiri digunakan untuk melihat apakah sisaan menyebar normal. Normal Q-Q Plot dan histogram di atas menunjukkan bahwa sisaan cenderung 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.

#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.97723, p-value = 0.2831

Berdasarkan uji formal Saphiro-Wilk 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 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.38175, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil DW Test, didapatkan nilai p-value = \(2.2e-16\), 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 Autokorelasi Cochrane-Orcutt
modelCO<-cochrane.orcutt(model)
modelCO

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = GDP ~ Tahun, data = data)
## 
##  number of interaction: 5
##  rho 0.787154
## 
## Durbin-Watson statistic 
## (original):    0.38175 , p-value: 1.44e-17
## (transformed): 1.88290 , p-value: 2.744e-01
##  
##  coefficients: 
##   (Intercept)         Tahun 
## -2.177382e+13  1.105977e+10

Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=--2.177382e+13 + 1.105977e+10 x_t\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(1.88290\) dan \(2.744e-01\). Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.331 < DW < 2.669\). 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.3409214\). Nilai tersebut dapat diketahui dengan syntax berikut.

#Rho optimum
rho<- modelCO$rho
rho

## [1] 0.7871539

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

data$GDP

##  [1]  15953203853  17354781171  18821353529  20371666154  22718426804
##  [6]  25000191834  27194141951  29517675437  31323283317  34016617942
## [11]  38092452061  41566412923  48954145809  59404966684  66013338965
## [16]  82885397621  89362071673  94468740851 104442351223 123386410161
## [21] 142092068281 129686938223 114380557731 105014356667 109201362581
## [26] 114123537582 150498057724 183009638351 206986674501 217948315625
## [31] 261846194499 274229034312 284321115595 212953336588 229033566615
## [36] 267305875261 291743811512 268146144678 270809066781 274072182417
## [41] 262835454367 242395852494 266849061836 334337212322 385118044877
## [46] 392218088879 423093437424 491252589217 517706149201 436537014294
## [51] 495812558843 574094112973 552483727283 586841821797 581964017237
## [56] 505103781350 515654671470 541018749769 555455371487 533879529188
## [61] 547054174236 639714956069 590409594949 593267701033

data$GDP[-1]

##  [1]  17354781171  18821353529  20371666154  22718426804  25000191834
##  [6]  27194141951  29517675437  31323283317  34016617942  38092452061
## [11]  41566412923  48954145809  59404966684  66013338965  82885397621
## [16]  89362071673  94468740851 104442351223 123386410161 142092068281
## [21] 129686938223 114380557731 105014356667 109201362581 114123537582
## [26] 150498057724 183009638351 206986674501 217948315625 261846194499
## [31] 274229034312 284321115595 212953336588 229033566615 267305875261
## [36] 291743811512 268146144678 270809066781 274072182417 262835454367
## [41] 242395852494 266849061836 334337212322 385118044877 392218088879
## [46] 423093437424 491252589217 517706149201 436537014294 495812558843
## [51] 574094112973 552483727283 586841821797 581964017237 505103781350
## [56] 515654671470 541018749769 555455371487 533879529188 547054174236
## [61] 639714956069 590409594949 593267701033

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

## 
## Call:
## lm(formula = gdp.trans ~ tahun.trans)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.462e+10 -4.696e+09  5.392e+07  5.162e+09  2.246e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.631e+12  1.424e+11  -32.53   <2e-16 ***
## tahun.trans  1.105e+10  3.357e+08   32.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.806e+09 on 61 degrees of freedom
## Multiple R-squared:  0.9467, Adjusted R-squared:  0.9459 
## F-statistic:  1084 on 1 and 61 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.175595e+13

b1

## tahun.trans 
## 11051485885

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,0.9, 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 1.288637e+23
## 2 0.2 1.084809e+23
## 3 0.3 9.129718e+22
## 4 0.4 7.731248e+22
## 5 0.5 6.652682e+22
## 6 0.6 5.894019e+22
## 7 0.7 5.455259e+22
## 8 0.8 5.336402e+22
## 9 0.9 5.537449e+22

Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.8. 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.4
rOpt <- seq(0.2,0.5, 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
## 301 0.500 6.652682e+22
## 300 0.499 6.661884e+22
## 299 0.498 6.671118e+22
## 298 0.497 6.680384e+22
## 297 0.496 6.689682e+22
## 296 0.495 6.699012e+22

#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 \(ρ=0.341\). 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 
## -8.752e+10 -2.044e+10 -2.413e+09  1.638e+10  9.469e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.378e+13  5.159e+11  -26.71   <2e-16 ***
## x            1.063e+10  3.929e+08   27.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.737e+10 on 61 degrees of freedom
## Multiple R-squared:  0.9231, Adjusted R-squared:  0.9218 
## F-statistic: 731.8 on 1 and 61 DF,  p-value: < 2.2e-16

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

## y = -2.090629e+13+10627397355x

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

#Deteksi autokorelasi
dwtest(modelHL)

## 
##  Durbin-Watson test
## 
## data:  modelHL
## DW = 0.80777, p-value = 2.856e-08
## alternative hypothesis: true autocorrelation is greater than 0

Hasil uji Durbin-Watson juga menunjukkan bawah nilai p-value sebesar \(2.856e-08\), di mana p-value < \(\alpha\)=5%. Artinya tolak \(H_0\) atau belum bukti menyatakan bahwa ada autokorelasi dalam data nilai IPM 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)

## Warning: Unknown or uninitialised column: `gdp`.

mseModelCO <- sseModelCO/length(data$gdp)

## Warning: Unknown or uninitialised column: `gdp`.

mseModelHL <- sseModelHL/length(data$gdp)

## Warning: Unknown or uninitialised column: `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.599725e+23          5.865918e+21      8.517653e+22
## MSE          Inf                   Inf               Inf

Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt memiliki SSE yang lebih rendah, sebesar \(5.865918e+21\) dan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar \(1.599725e+23\).

Simpulan

Autokorelasi yang terdapat pada data GDP terjadi karena unsur penyusunnya yang saling erat kaitannya, juga nilai yang sangat curam sehingga menyebabkan sulitnya melakukan penanganan. Walaupun ditangani dengan dua metode, yaitu Cochrane-Orcutt dan Hildreth-Lu, tetapi tidak menjadikan data tersebut menjadi tidak ada autokorelasi.

Tugas Pertemuan 2

Delita Nur Hasanah