MPDW Minggu 2

Tugas Minggu 2 MPDW

Aljazair ( Algeria)

Aljazair adalah Negara Terbesar di Afrika dan Wilayah Mediterania. Negara ini terletak dibagian utara Benua Afrika. Negara kelahiran Riyad Mahrez ini memiliki sejarah yang kaya, termasuk sebagian dari kekaisaran Romawi, berpindah ke kekuasaan Arab dan sempat memasuki kolonisasi Prancis. Ekonomi Aljazair bergantung pada ekspor minyak dan gas alam yang menjadikan Alajzair sebagai salah satu negara yang memiliki cadangan energi terbesar di Dunia.

Aljazair dikategorikan sebagai negara berpenghasilan menengah ke atas oleh Bank Dunia. Meskipun demikian, Aljazair telah menunjukkan kemajuan signifikan dalam Indikator Pembangunan Manusia (IPM) dan menempati peringkat ke-83 dari 188 negara, Maka dari itu pada kesempatan kali ini saya akan melakukan analisis Regresi Deret Waktu mengenai perkembangan GDP Aljazair dari tahun 1960 sampai 2023.

Input Library

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.3
## 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
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.2

Input Data

dataalg<- read_excel("C:/Users/Admin/Downloads/Adriano Excel Putra/data algeria.xlsx")
dataalg
# Rename Data menjadi GDP dan Time menjadi Tahun 
dataalg<- dataalg %>% rename(GDP=Data, Tahun=Time)
dataalg
# Buat Tahun menjadi numeric
dataalg$GDP<- as.numeric(dataalg$GDP)
#Membentuk objek time series
dataalgts<-ts(dataalg$GDP)
dataalgts
## Time Series:
## Start = 1 
## End = 64 
## Frequency = 1 
##  [1]   2723615451   2434747056   2001444544   2702982018   2909316435
##  [6]   3136284307   3039859187   3370870376   3852147027   4257253264
## [11]   4863526897   5077183094   6766743957   8707858912  13209871626
## [16]  15557902754  17728240932  20972113685  26364491313  33243706860
## [21]  42345829079  44348590461  45207167470  48801369800  53698548293
## [26]  57937868670  63692007897  66745818375  59089396860  55634721573
## [31]  62048507531  45715676428  48003133347  49945584453  42543176829
## [36]  41764291672  46941554225  48177612042  48187781984  48640671735
## [41]  54790398570  59413400924  61516103406  73482264191  91913680985
## [46] 107046618670 123084258693 142482739810 180383848331 150317292079
## [51] 177785053940 218331946925 227143746076 229701430292 238942664193
## [56] 187493855609 180763839522 189880896903 194554483656 193459662091
## [61] 164873415325 186265418571 225560256622 239899491128
#Membuat plot time series
ts.plot(dataalgts, xlab="Time Period ", ylab="IPM", main= "Time Series Plot of GDP")
points(dataalgts)

Visualiasi yang didapat dengan menggunakan plot TIme Series menunjukkan ada kecenderungan peningkatan GDP Aljazair dari waktu ke waktu ( Dimulai dari 1960 ke 2023). Namun Tren Positif ini terkadang juga diikuti dengan fluktuasi yang signfikan.

  • Stabil Awal: Pada awal time series (periode 0 hingga sekitar 20), nilai IPM cenderung stabil dengan sedikit peningkatan.

  • Fluktuasi Sedang: Di sekitar periode 30 hingga 40, terdapat fluktuasi yang lebih sering dan intensitas perubahan nilai meningkat.

  • Peningkatan Tajam: Setelah periode ke-45, ada peningkatan tajam yang menunjukkan akselerasi pertumbuhan IPM.

  • Puncak dan Koreksi: Di periode mendekati 60, meskipun masih ada tren meningkat, terdapat beberapa koreksi atau penurunan sementara sebelum kembali naik.

Smoothing Data

Plotting data 80:20

trainingalg<- dataalg[1:51,2]
testingalg<-dataalg[52:64,2]

traints<- ts(trainingalg$GDP)
testts<- ts(testingalg$GDP)

Data yang didapat menunjukkan pola Trend, Maka akan dilakukan DMA dan DES untuk pemulusannya.

DMA ( Double Moving Average)

datasma <- SMA(dataalgts, n=4)
DMA <- SMA(datasma, n = 4)
At <- 2*datasma - DMA
Bt <- 2/(4-1)*(datasma - DMA)
dt.dma<- At+Bt
dt.ramal<- c(NA, dt.dma)

t = 1 :14
f = c()

for (i in t) {
  f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}
datagabung <- cbind(aktual = c(dataalgts,rep(NA,14)), 
                pemulusan1 = c(datasma,rep(NA,14)),
                pemulusan2 = c(dt.dma, rep(NA,14)),
                At = c(At, rep(NA,14)), 
                Bt = c(Bt,rep(NA,14)),
                ramalan = c(dt.ramal, f[-1]))
datagabung
##             aktual   pemulusan1   pemulusan2           At           Bt
##  [1,]   2723615451           NA           NA           NA           NA
##  [2,]   2434747056           NA           NA           NA           NA
##  [3,]   2001444544           NA           NA           NA           NA
##  [4,]   2702982018   2465697267           NA           NA           NA
##  [5,]   2909316435   2512122513           NA           NA           NA
##  [6,]   3136284307   2687506826           NA           NA           NA
##  [7,]   3039859187   2947110487   3437112509   3241111700    196000809
##  [8,]   3370870376   3114082576   3612210870   3412959552    199251317
##  [9,]   3852147027   3349790224   3891736384   3674957920    216778464
## [10,]   4257253264   3630032464   4246330007   3999810989    246519017
## [11,]   4863526897   4085949391   4987592269   4626935118    360657151
## [12,]   5077183094   4512527570   5542448667   5130480228    411968439
## [13,]   6766743957   5241176803   6697435546   6114932049    582503497
## [14,]   8707858912   6353828215   8529591082   7659285935    870305147
## [15,]  13209871626   8440414397  12279460482  10743842048   1535618434
## [16,]  15557902754  11060594312  16538245780  14347185193   2191060587
## [17,]  17728240932  13800968556  20279330533  17687985742   2591344791
## [18,]  20972113685  16867032249  24074998700  21191812120   2883186580
## [19,]  26364491313  20155687171  27963381503  24840303770   3123077733
## [20,]  33243706860  24577138198  34122024288  30304069852   3817954436
## [21,]  42345829079  30731535234  43479346936  38380222256   5099124681
## [22,]  44348590461  36575654428  50851738879  45141305099   5710433780
## [23,]  45207167470  41286323468  54609091194  49279984103   5329107090
## [24,]  48801369800  45175739203  56398116069  51909165322   4488950746
## [25,]  53698548293  48013919006  56765602306  53264928986   3500673320
## [26,]  57937868670  51411238558  59643627725  56350672058   3292955667
## [27,]  63692007897  56032448665  65822635844  61906560972   3916074871
## [28,]  66745818375  60518560809  71392759225  67043079858   4349679366
## [29,]  59089396860  61866272951  69214844126  66275415656   2939428470
## [30,]  55634721573  61290486177  63563059553  62654030203    909029351
## [31,]  62048507531  60879611085  60447741634  60620489415   -172747780
## [32,]  45715676428  55622075598  48467849174  51329539744  -2861690570
## [33,]  48003133347  52850509720  44833574845  48040348795  -3206773950
## [34,]  49945584453  51428225440  45150092072  47661345419  -2511253347
## [35,]  42543176829  46551892764  38116420904  41490609648  -3374188744
## [36,]  41764291672  45564046575  39673009826  42029424525  -2356414700
## [37,]  46941554225  45298651795  42111897880  43386599446  -1274701566
## [38,]  48177612042  44856658692  43671402418  44145504928   -474102510
## [39,]  48187781984  46267809981  47552840348  47038828201    514012147
## [40,]  48640671735  47986904997  51127569381  49871303627   1256265754
## [41,]  54790398570  49949116083  54422438824  52633109728   1789329097
## [42,]  59413400924  52758063303  58620712824  56275653016   2345059808
## [43,]  61516103406  56090143659  63413621406  60484230307   2929391099
## [44,]  73482264191  62300541773  74010667720  69326617341   4684050379
## [45,]  91913680985  71581362376  89746086708  82480196975   7265889732
## [46,] 107046618670  83489666813 108696730409  98613904971  10082825438
## [47,] 123084258693  98881705635 131912349777 118700092120  13212257657
## [48,] 142482739810 116131824539 155482965704 139742509238  15740456466
## [49,] 180383848331 138249366376 186684742268 167310591911  19374150357
## [50,] 150317292079 149067034728 188207954576 172551586637  15656367939
## [51,] 177785053940 162742233540 198066598113 183936852284  14129745829
## [52,] 218331946925 181704535319 221310773366 205468278147  15842495219
## [53,] 227143746076 193394509755 229506895454 215061941175  14444954280
## [54,] 229701430292 213240544308 255690691938 238710632886  16980059052
## [55,] 238942664193 228529946871 269050884885 252842509680  16208375205
## [56,] 187493855609 220820424043 232193870373 227644491841   4549378532
## [57,] 180763839522 209225447404 194677708650 200496804152  -5819095502
## [58,] 189880896903 199270314057 173951615662 184079095020 -10127479358
## [59,] 194554483656 188173268923 161174777783 171974174239 -10799396456
## [60,] 193459662091 189664720543 178133525228 182746003354  -4612478126
## [61,] 164873415325 185692114494 177345464476 180684124483  -3338660007
## [62,] 186265418571 184788244911 180969341066 182496902604  -1527561538
## [63,] 225560256622 192539688152 199820515031 196908184279   2912330752
## [64,] 239899491128 204149645411 224745015694 216506867581   8238148113
## [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
## [70,]           NA           NA           NA           NA           NA
## [71,]           NA           NA           NA           NA           NA
## [72,]           NA           NA           NA           NA           NA
## [73,]           NA           NA           NA           NA           NA
## [74,]           NA           NA           NA           NA           NA
## [75,]           NA           NA           NA           NA           NA
## [76,]           NA           NA           NA           NA           NA
## [77,]           NA           NA           NA           NA           NA
## [78,]           NA           NA           NA           NA           NA
##            ramalan
##  [1,]           NA
##  [2,]           NA
##  [3,]           NA
##  [4,]           NA
##  [5,]           NA
##  [6,]           NA
##  [7,]           NA
##  [8,]   3437112509
##  [9,]   3612210870
## [10,]   3891736384
## [11,]   4246330007
## [12,]   4987592269
## [13,]   5542448667
## [14,]   6697435546
## [15,]   8529591082
## [16,]  12279460482
## [17,]  16538245780
## [18,]  20279330533
## [19,]  24074998700
## [20,]  27963381503
## [21,]  34122024288
## [22,]  43479346936
## [23,]  50851738879
## [24,]  54609091194
## [25,]  56398116069
## [26,]  56765602306
## [27,]  59643627725
## [28,]  65822635844
## [29,]  71392759225
## [30,]  69214844126
## [31,]  63563059553
## [32,]  60447741634
## [33,]  48467849174
## [34,]  44833574845
## [35,]  45150092072
## [36,]  38116420904
## [37,]  39673009826
## [38,]  42111897880
## [39,]  43671402418
## [40,]  47552840348
## [41,]  51127569381
## [42,]  54422438824
## [43,]  58620712824
## [44,]  63413621406
## [45,]  74010667720
## [46,]  89746086708
## [47,] 108696730409
## [48,] 131912349777
## [49,] 155482965704
## [50,] 186684742268
## [51,] 188207954576
## [52,] 198066598113
## [53,] 221310773366
## [54,] 229506895454
## [55,] 255690691938
## [56,] 269050884885
## [57,] 232193870373
## [58,] 194677708650
## [59,] 173951615662
## [60,] 161174777783
## [61,] 178133525228
## [62,] 177345464476
## [63,] 180969341066
## [64,] 199820515031
## [65,] 224745015694
## [66,] 232983163807
## [67,] 241221311920
## [68,] 249459460033
## [69,] 257697608146
## [70,] 265935756259
## [71,] 274173904372
## [72,] 282412052485
## [73,] 290650200598
## [74,] 298888348711
## [75,] 307126496823
## [76,] 315364644936
## [77,] 323602793049
## [78,] 331840941162
ts.plot(datagabung[,1], xlab="Time Period", ylab="IPM", 
        main="DMA N=4 Data IPM")
points(datagabung[,1])

lines(datagabung[,3], col="green", lwd=2)
lines(datagabung[,6], col="red", lwd=2)
legend("topleft", c("data aktual", "data pemulusan", "data peramalan"), 
       lty=8, col=c("black", "green", "red"), cex=0.8)

Akurasi Training DMA

#Menghitung nilai keakuratan data latih

error_train.DMA = traints - dt.ramal[1:length(traints)]
SSE_train.DMA = sum(error_train.DMA[8:length(traints)]^2)
MSE_train.DMA = mean(error_train.DMA[8:length(traints)]^2)
MAPE_train.DMA = mean(abs((error_train.DMA[8:length(traints)] / traints[8:length(traints)]) * 100))

akurasi_train.DMA <- matrix(c(SSE_train.DMA, MSE_train.DMA, MAPE_train.DMA))
row.names(akurasi_train.DMA) <- c("SSE", "MSE", "MAPE")
colnames(akurasi_train.DMA) <- c("Akurasi m = 4")
akurasi_train.DMA
##      Akurasi m = 4
## SSE   4.084386e+21
## MSE   9.282695e+19
## MAPE  1.140471e+01

Akurasi Testing DMA

error_test.dma = testts-datagabung[52:64,6]
SSE_test.dma = sum(error_test.dma^2)
MSE_test.dma = mean(error_test.dma^2)
MAPE_test.dma = mean(abs((error_test.dma/testts*100)))

akurasi_test.dma <- matrix(c(SSE_test.dma, MSE_test.dma, MAPE_test.dma))
row.names(akurasi_test.dma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi_test.dma) <- c("Akurasi m = 4")
akurasi_test.dma
##      Akurasi m = 4
## SSE   1.536171e+22
## MSE   1.181670e+21
## MAPE  1.307730e+01

Di atas adalah hasil pemulusan data menggunakan DMA dengan nilai m=4. Dari hasil tersebut, terlihat bahwa nilai Mean Absolute Percentage Error (MAPE) untuk data Training dan Data Testing cukup bervariatif namun masih berada di angka puluhan. Memang tidak dibawah 10 persen, namun angka ini menunjukkan bahwa metode ini masih dapat memprediksi dengan cukup baik dimana Data Training mendapat akurasi 11.40 % dan Testing mendapat 13.07 %.

DES

trainingalg<- dataalg[1:51,2]
testingalg<-dataalg[52:64,2]

#data time series
training.ts<-ts(trainingalg)
testing.ts<-ts(testingalg,start=51)
plot(training.ts, col="blue",main="Plot data training")
points(training.ts)

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

Disini akan dicoba untuk melakukan peramalan pada data testing yang tidak berbentuk seperti data Trend , namun lebih ke data yang fluktuatif. Tentu saja Metode ini bisa saja mendapatkan hasil prediksi yang tidak terlalu akurat karena memang Metode DES sendiri memprediksi lewat pembobotan dari data-data sebelumnya, dimana apabila model data Training dan Testingnya cukup berbeda, dapat mengakibatkan prediksi yang kurang fit.

#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: 0.5495219
##  beta : 0.5609428
##  gamma: FALSE
## 
## Coefficients:
##           [,1]
## a 179641874124
## b   9938519817
plot(des.opt)
legend("topleft", c("Data Aktual", "Peramalan"), col = c("black", "red"), 
       lty = c(1,1))

ramalandesopt<- forecast(des.opt, h=14)
ramalandesopt
##    Point Forecast        Lo 80        Hi 80        Lo 95        Hi 95
## 52   189580393941 178166186775 200994601107 172123870491 207036917391
## 53   199518913758 184480838455 214556989060 176520162259 222517665256
## 54   209457433574 189375621743 229539245405 178744952710 240169914438
## 55   219395953391 193195776853 245596129928 179326241190 259465665591
## 56   229334473207 196162003992 262506942422 178601559536 280067386878
## 57   239272993024 198408236710 280137749338 176775741438 301770244610
## 58   249211512840 200020509726 298402515955 173980365643 324442660038
## 59   259150032657 201058287091 317241778223 170306374572 347993690742
## 60   269088552474 201565405277 336611699670 165820810319 372356294628
## 61   279027072290 201575880864 356478263717 160575696757 397478447823
## 62   288965592107 201117173831 376814010382 154613030300 423318153914
## 63   298904111923 200212132613 397596091234 147967754611 449840469235
## 64   308842631740 198880218317 418805045163 140669632887 477015630593
## 65   318781151557 197138313283 440423989830 132744484439 504817818674
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,] -144434117
## [4,]  969863337
## [5,]  677666115
## [6,]  357780282

Akurasi Training DES

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.830147e+21
## MSE                     7.510092e+19
## MAPE                    9.617202e+00
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt
## Time Series:
## Start = 52 
## End = 63 
## Frequency = 1 
##         testing.ts
##  [1,] -37563352135
##  [2,] -30182516535
##  [3,] -29485230618
##  [4,]  31902097781
##  [5,]  48570633685
##  [6,]  49392096121
##  [7,]  54657029185
##  [8,]  65690370566
##  [9,] 104215137148
## [10,]  92761653720
## [11,]  63405335485
## [12,]  59004620796

Akurasi Testing DES

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.329055e+21
## MSE                     3.329055e+21
## MAPE                    2.612912e+01

Setelah didapatkan nilai akurasi untuk metode DMA dan DES, selanjutnya akan dibandingkan keakuratan antar metode keduanya.

Perbandingan Model Training DMA dan DES

cbind(akurasi_train.DMA, akurasides.opt)
##      Akurasi m = 4 Akurasi lamda dan gamma optimum
## SSE   4.084386e+21                    3.830147e+21
## MSE   9.282695e+19                    7.510092e+19
## MAPE  1.140471e+01                    9.617202e+00

Berdasarkan Perbandingan Model Training ( Data Latih), didapatkan MAPE pada Metode DES lebih kecil daripada MAPE metode DMA yaitu 9.617 % banding 11.40 %, Walaupun tidak terlalu signifikan namun tetap MAPE lebih baik secara substansial berada pada metode DES.

Perbandingan Model Testing DMA dan DES

cbind(akurasi_test.dma,akurasiDesTesting)
##      Akurasi m = 4 Akurasi lamda dan gamma optimum
## SSE   1.536171e+22                    3.329055e+21
## MSE   1.181670e+21                    3.329055e+21
## MAPE  1.307730e+01                    2.612912e+01

Hal yang berbeda didapat pada Model Testing ( Data Uji), didapatkan MAPe pada Metode DMA lebih kecil daripada MAPE dari Metode DES yaitu 13.07 % banding 26.12 %. Hal ini tentunya melawan hasil data latih tadi. Dimana dalam konteks forecasting, tentunya kita ingin mendapatkan MAPE yang kecil pada data testingnya.

nilai MAPE yang besar pada metode DES ini karena terdapat perbedaan pola pada data Training dan Testing yang signifikan. Hal ini tentunya dipengaruhi oleh pembagian data training dan data testing pada tahap awal, ketika kita merubah polanya menjadi dari 80:20 menjadi 70:30 ataupun 60:40. Hal ini terjadi karena Pada Metode SES, dilakukan pembobotan dari data training untuk memprediksi data di periode berikutnya berbeda dengan Moving Average yang melihat data secara keseluruhan.

Plot

plot(dataalg, pch = 20, col = "blue",
     main = "Scatter Plot Tahun vs Nilai GDP",
     xlab = "Tahun",
     ylab = "GDP")

Dapat dilihat dari Scatterplotnya bahwa hubungannya menyebar lebih ke kuadratik daripada linear. Namun akan dilakukan pemodelan linear terlebih dahulu pada tahap ini.

Korelasi

cor(dataalg$Tahun,dataalg$GDP)
## [1] 0.9024415

Korelasi yang didapat sangat tinggi dan bernilai positif artinya. Hubungan antara Pertambahan Tahun dan GDP sangat erat, dimana semakin bertambahnya tahun, maka GDP semakin bertambah juga.

Regresi Deret Waktu

model<- lm(GDP~Tahun, data = dataalg)
summary(model)
## 
## Call:
## lm(formula = GDP ~ Tahun, data = dataalg)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -5.878e+10 -1.616e+10  1.400e+09  1.709e+10  7.611e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.278e+12  4.461e+11  -16.32   <2e-16 ***
## Tahun        3.694e+09  2.240e+08   16.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.31e+10 on 62 degrees of freedom
## Multiple R-squared:  0.8144, Adjusted R-squared:  0.8114 
## F-statistic: 272.1 on 1 and 62 DF,  p-value: < 2.2e-16

Model Yang Dihasilkan :

\[\hat{GDP} = -7.278 \times 10^{12} + 3.694 \times 10^{9} \cdot Tahun + \epsilon\]

Model ini menunjukkan GDP diperkirakan meningkat secara linear setiap tahun, dengan kenaikan sebesar 3.694 Milliar setiap tahun. Dengan NIlai konstanta -7.278 * 10 ^12 . Konstanta/ koefisien ini secara praktis lebih bersifat matematis daripada interpretatif. Koefisien yang negatif ini juga dapat menunjukkan bahwa GDP tahun-tahun sebelumnya cenderung lebih rendah daripada tahun-tahun berikutnya.

Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value (2.2*10^-16) < α (5%) artinya sangat jauh lebih kecil dari nilai alphanya. Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model.

Hasil uji-t parsial kedua parameter regresi, yaitu intersep dan koefisien regresi juga menunjukkan hal yang sama, yaitu memiliki p-value < α (5%) sehingga nyata dalam taraf 5%.

Selanjutnya dapat dilihat juga nilai R2=0.81 . Artinya, sebesar 81.44% 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.

Pengujian Asumsi

#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 sebelah kiri digunakan untuk mengevaluasi apakah residual berdistribusi normal.

Normal Q-Q Plot di atas menunjukkan bahwa residual tampaknya berdistribusi normal dimana sebagian besar data berada di garis teoritis normal. , Namukan histogram residual tidak terlalu menunjukkan distribusi normal.

Selanjutnya, dua plot di sebelah kanan digunakan untuk memeriksa autokorelasi.

Plot Residual vs Fitted Value dan Plot Residual vs Order menunjukkan adanya pola pada residual, dimana seharusnya tidak tercipat pola .

Untuk analisis lebih lanjut, akan dilakukan uji formal untuk memeriksa normalitas residual serta plot ACF dan PACF untuk menilai adanya autokorelasi.

Normalitas

ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.12496, p-value = 0.2488
## alternative hypothesis: two-sided

Dari hasil uji Kolmogorov-Smirnov, diperoleh nilai p-value lebih besar dari 0.05 sehingga residual berdistribusi normal. Uji Kolmogorov-Smirnov dipakai karena didapatkan sampel data yang cukup besar yaitu diatas 50.

(Drezner, Z., Turel, O., & Zerom, D. (2010). A Modified Kolmogorov–Smirnov Test for Normality. Communications in Statistics - Simulation and Computation, 39(4), 693–704. https://doi.org/10.1080/03610911003615816)

Autokorelasi

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

Dari Plot ACF dan PACF di atas, terlihat bahwa terdapat autokorelasi pada residual, dimana terdapat lag yang melebih batas wajar, seidaknya ada beberapa lag yang melewati batas wajar. Oleh karena itu, kita perlu melakukan uji formal untuk memeriksa autokorelasi.

dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.17708, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil DW Test, didapatkan nilai DW 0.17708. Tentunya nilai DW ini sangat kecil sekali tanpa kita perlu melihat nilai tabel Durbin Watsonnya. ( NIlai Tabel DW untuk n=64 , DU = 1.6268 yang artinya DW ini berada jauh dibawah batas toleransi. Namun Nilai DW ini juga berada di antara DU dan batasnya DU-DL ( 1.6268-1.5636= 0.0629), sehingga Nilai DW berada di daerah inkonklusif.

Namun, didapat p-value yang sangat jauh lebih kecil dari batas toleransi alpha 5 persen. Dapat disimpulkan bahwa tolak H-0, bahwa cukup bukti terdapat korelasi pada model yang ada. Oleh karena itu, kita perlu melakukan penanganan pada autokorelasi dengan menggunakan metode-metode yang ada.

Metode Cochrane-Orcutt

modelCO<-cochrane.orcutt(model)
modelCO
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = GDP ~ Tahun, data = dataalg)
## 
##  number of interaction: 18
##  rho 0.914271
## 
## Durbin-Watson statistic 
## (original):    0.17708 , p-value: 4.785e-28
## (transformed): 1.74616 , p-value: 1.253e-01
##  
##  coefficients: 
##   (Intercept)         Tahun 
## -1.012853e+13  5.117985e+09

Model yang dihasilkan :

\[ \hat{GDP} = -1.0128 \times 10^{13} + 5.117 \times 10^{9} \cdot Tahun + \epsilon\]

Model Cohcrane-Orcutt meningkatan nilai DW dengan signfikan menjadi 1.74616. Tentunya ini berada di daerah Kurang cukup bukti untuk menyatakan adanya autokorelasi ( DU < DW < 4- DU) Dimana DU disini merujuk pada angka 1.6268.

Didapat bahwa P-valuenya meningkat dratis pada model yang dihasilkan menjadi 0.1 , hal ini mengindikasikan pada taraf 5 persen, kita dapat menerima H0., sehingga dapat disimpulkan bahwa model yang dihasilkan sudah tidak memiliki autokorelasi.

Normalitas CO

sisaan<- modelCO$residuals
ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.10989, p-value = 0.3937
## alternative hypothesis: two-sided

P-Value > 0.05 yang mengindikasikan bahwa Data berdistribusi normal

Perhitungan Manual

rho<- modelCO$rho
rho
## [1] 0.9142712

Didapatkan Rho optimum pada model Cocran Orchut adalah 0.91427

GDP <- dataalg$GDP
Tahun<- dataalg$Tahun
GDP.trans <- GDP[-1] - GDP[-64] * rho
Tahun.trans <- Tahun[-1] - Tahun[-64] * rho
modelCOmanual <- lm(GDP.trans ~ Tahun.trans)
summary(modelCOmanual)
## 
## Call:
## lm(formula = GDP.trans ~ Tahun.trans)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -5.144e+10 -4.467e+09  9.649e+08  3.191e+09  3.707e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.683e+11  1.881e+11  -4.617 2.06e-05 ***
## Tahun.trans  5.118e+09  1.095e+09   4.672 1.69e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.355e+10 on 61 degrees of freedom
## Multiple R-squared:  0.2636, Adjusted R-squared:  0.2515 
## F-statistic: 21.83 on 1 and 61 DF,  p-value: 1.687e-05
#Mencari Penduga Koefisien Regresi setelah Transformasi ke Persamaan Awal
b0bintang <- modelCOmanual$coefficients[-2]
b0 <- b0bintang/(1-rho)
b1 <- modelCOmanual$coefficients[-1]
b0
##   (Intercept) 
## -1.012853e+13
b1
## Tahun.trans 
##  5117985119

Dapat dilihat bahwa nilai B0 dan B1 kembali ke bentuk model awalnya.

Metode Hildreth-Lu

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

Pertama-tama akan dicari di mana kira-kira ρ yang menghasilkan SSE minimum. Pada hasil di atas terlihat ρ minimum ketika di rentang 0.6 - 0.9. Namun, hasil tersebut masih kurang teliti sehingga akan dicari kembali ρ yang lebih optimum dengan ketelitian yang lebih. Jika sebelumnya jarak antar ρ yang dicari adalah 0.1, kali ini jar ak antar ρ yang dicari adalah 0.01.

#Rho optimal di sekitar 0.4
rOpt <- seq(0.6,0.9, by= 0.001)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model))}))
head(tabOpt[order(tabOpt$SSE),])
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.899, y=1.122063e22, labels = "rho=0.341", cex = 0.8)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai ρ optimum, yaitu saat SSE terkecil terdapat pada nilai ρ=0.899 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.

modelHL <- hildreth.lu.func(0.899, model)
summary(modelHL)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -5.024e+10 -4.958e+09  1.049e+09  3.619e+09  3.756e+10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -9.772e+11  1.880e+11  -5.197 2.48e-06 ***
## x            4.893e+09  9.304e+08   5.259 1.97e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.356e+10 on 61 degrees of freedom
## Multiple R-squared:  0.3119, Adjusted R-squared:  0.3007 
## F-statistic: 27.65 on 1 and 61 DF,  p-value: 1.968e-06
cat("y = ", coef(modelHL)[1]/(1-0.899), "+", coef(modelHL)[2],"x", sep = "")
## y = -9.675112e+12+4892573024x
dwtest(modelHL)
## 
##  Durbin-Watson test
## 
## data:  modelHL
## DW = 1.7179, p-value = 0.1033
## alternative hypothesis: true autocorrelation is greater than 0

Dari hasil pengujian Durbin-Watson pada Model Hilderth-Lu , Didapatkan Nilai DW juga ikut meningkat menjadi 1.7179 dan p-valuenya menjadi 0.103. Tentunya kedua hal ini mengindikasikan bahwa cukup bukti untuk menyatakan model ini tidak memiliki autokorelasi. Dimana p-valuenya lebih besar dari alpha 0.05 sehingga cukup bukti untuk menyatakan tidak ada autokorelasi pada model Hilderth-Lu.

Dari sini dapat ditarik kesimpulan bahwa dengan menggunakan Model Cochrane-Orcutt dan Hilderth-Lu, masalah autokorelasi dapat terselesaikan.

Normalitas HL

sisaan<- modelHL$residuals
ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.1681, p-value = 0.05027
## alternative hypothesis: two-sided

Dapat dilihat bahwa P-value yang dihasilkan masih lebih besar dari batas toleransi 0.05 yang artinya masih dapat dikatakan Model Hilderth Lu bertisribusi normal.

Perbandingan ketiga model

sseModelawal <- anova(model)$`Sum Sq`[-1]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-1]
sseModelHL <- anova(modelHL)$`Sum Sq`[-1]
mseModelawal <- sseModelawal/length(GDP)
mseModelCO <- sseModelCO/length(GDP)
mseModelHL <- sseModelHL/length(GDP)
akurasi <- matrix(c(sseModelawal,sseModelCO,sseModelHL,
                    mseModelawal,mseModelCO,mseModelHL),nrow=2,ncol=3,byrow = T)
colnames(akurasi) <- c("Model Awal", "Model Cochrane-Orcutt", "Model Hildreth-Lu")
row.names(akurasi) <- c("SSE","MSE")
akurasi
##       Model Awal Model Cochrane-Orcutt Model Hildreth-Lu
## SSE 6.793502e+22          1.120526e+22      1.122063e+22
## MSE 1.061485e+21          1.750822e+20      1.753223e+20

Berdasarkan Hasil Perbandingan ketiga model, dapat disimpulkan bahwa Model Ochrance-Orcutt dan Hilderth-Lu mampu menurunkan SSE dari model awal, hal ini juga diikuti bahwa setelah dilakukan transformasi model, autokorelasi dapat teratasi. Jika dilihat Model CO dan Hildert-Lu hampir mirip di nilai SSE nya namun nilai yang lebih kecil didiapat di Model Cochrane-Orcutt

Lihat saja bahwa Model dengan Autokrelasi mendapat nilai r-squared 81 persen dimana model dapat dipercaya, padahal model yang aman sendiri mendapat r-squared yang rendah 30 % dan 24%, dimana terjadi overfittng pada model awal.

Simpulan

Autokorelasi pada GDP sering terjadi karena memang pertumbuhan ekonomi sendiri tersusun atas berbagai hal, utamanya dalam GDP, seperti ekspor , impor , jumlah tenaga kerja dan sebagainya. Hal ini berpotensi untuk memunculkan autokorelasi, Autokorelasi ini menjadi masalah karena akan membuat model menjadi tidak fit dan meningkatkan galat, Masalah ini mampu diatasi dengan menggunakan Model Cohrane- Orcutt dan Model Hilderth-Lu.