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/(3-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   3535112914   3241111700    294001214
##  [8,]   3370870376   3114082576   3711836528   3412959552    298876976
##  [9,]   3852147027   3349790224   4000125616   3674957920    325167696
## [10,]   4257253264   3630032464   4369589515   3999810989    369778526
## [11,]   4863526897   4085949391   5167920845   4626935118    540985727
## [12,]   5077183094   4512527570   5748432886   5130480228    617952658
## [13,]   6766743957   5241176803   6988687295   6114932049    873755246
## [14,]   8707858912   6353828215   8964743655   7659285935   1305457720
## [15,]  13209871626   8440414397  13047269699  10743842048   2303427651
## [16,]  15557902754  11060594312  17633776073  14347185193   3286590880
## [17,]  17728240932  13800968556  21575002928  17687985742   3887017186
## [18,]  20972113685  16867032249  25516591990  21191812120   4324779871
## [19,]  26364491313  20155687171  29524920369  24840303770   4684616599
## [20,]  33243706860  24577138198  36031001506  30304069852   5726931654
## [21,]  42345829079  30731535234  46028909277  38380222256   7648687021
## [22,]  44348590461  36575654428  53706955769  45141305099   8565650670
## [23,]  45207167470  41286323468  57273644739  49279984103   7993660636
## [24,]  48801369800  45175739203  58642591442  51909165322   6733426119
## [25,]  53698548293  48013919006  58515938966  53264928986   5251009980
## [26,]  57937868670  51411238558  61290105558  56350672058   4939433500
## [27,]  63692007897  56032448665  67780673279  61906560972   5874112307
## [28,]  66745818375  60518560809  73567598908  67043079858   6524519049
## [29,]  59089396860  61866272951  70684558361  66275415656   4409142705
## [30,]  55634721573  61290486177  64017574229  62654030203   1363544026
## [31,]  62048507531  60879611085  60361367744  60620489415   -259121670
## [32,]  45715676428  55622075598  47037003889  51329539744  -4292535854
## [33,]  48003133347  52850509720  43230187870  48040348795  -4810160925
## [34,]  49945584453  51428225440  43894465398  47661345419  -3766880021
## [35,]  42543176829  46551892764  36429326532  41490609648  -5061283116
## [36,]  41764291672  45564046575  38494802476  42029424525  -3534622050
## [37,]  46941554225  45298651795  41474547097  43386599446  -1912052349
## [38,]  48177612042  44856658692  43434351163  44145504928   -711153764
## [39,]  48187781984  46267809981  47809846421  47038828201    771018220
## [40,]  48640671735  47986904997  51755702258  49871303627   1884398631
## [41,]  54790398570  49949116083  55317103373  52633109728   2683993645
## [42,]  59413400924  52758063303  59793242728  56275653016   3517589712
## [43,]  61516103406  56090143659  64878316955  60484230307   4394086648
## [44,]  73482264191  62300541773  76352692909  69326617341   7026075568
## [45,]  91913680985  71581362376  93379031574  82480196975  10898834599
## [46,] 107046618670  83489666813 113738143128  98613904971  15124238158
## [47,] 123084258693  98881705635 138518478606 118700092120  19818386485
## [48,] 142482739810 116131824539 163353193936 139742509238  23610684699
## [49,] 180383848331 138249366376 196371817446 167310591911  29061225535
## [50,] 150317292079 149067034728 196036138546 172551586637  23484551909
## [51,] 177785053940 162742233540 205131471028 183936852284  21194618744
## [52,] 218331946925 181704535319 229232020975 205468278147  23763742828
## [53,] 227143746076 193394509755 236729372594 215061941175  21667431419
## [54,] 229701430292 213240544308 264180721464 238710632886  25470088578
## [55,] 238942664193 228529946871 277155072488 252842509680  24312562808
## [56,] 187493855609 220820424043 234468559639 227644491841   6824067798
## [57,] 180763839522 209225447404 191768160899 200496804152  -8728643253
## [58,] 189880896903 199270314057 168887875983 184079095020 -15191219037
## [59,] 194554483656 188173268923 155775079555 171974174239 -16199094684
## [60,] 193459662091 189664720543 175827286165 182746003354  -6918717189
## [61,] 164873415325 185692114494 175676134473 180684124483  -5007990010
## [62,] 186265418571 184788244911 180205560297 182496902604  -2291342307
## [63,] 225560256622 192539688152 201276680407 196908184279   4368496127
## [64,] 239899491128 204149645411 228864089750 216506867581  12357222169
## [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,]   3535112914
##  [9,]   3711836528
## [10,]   4000125616
## [11,]   4369589515
## [12,]   5167920845
## [13,]   5748432886
## [14,]   6988687295
## [15,]   8964743655
## [16,]  13047269699
## [17,]  17633776073
## [18,]  21575002928
## [19,]  25516591990
## [20,]  29524920369
## [21,]  36031001506
## [22,]  46028909277
## [23,]  53706955769
## [24,]  57273644739
## [25,]  58642591442
## [26,]  58515938966
## [27,]  61290105558
## [28,]  67780673279
## [29,]  73567598908
## [30,]  70684558361
## [31,]  64017574229
## [32,]  60361367744
## [33,]  47037003889
## [34,]  43230187870
## [35,]  43894465398
## [36,]  36429326532
## [37,]  38494802476
## [38,]  41474547097
## [39,]  43434351163
## [40,]  47809846421
## [41,]  51755702258
## [42,]  55317103373
## [43,]  59793242728
## [44,]  64878316955
## [45,]  76352692909
## [46,]  93379031574
## [47,] 113738143128
## [48,] 138518478606
## [49,] 163353193936
## [50,] 196371817446
## [51,] 196036138546
## [52,] 205131471028
## [53,] 229232020975
## [54,] 236729372594
## [55,] 264180721464
## [56,] 277155072488
## [57,] 234468559639
## [58,] 191768160899
## [59,] 168887875983
## [60,] 155775079555
## [61,] 175827286165
## [62,] 175676134473
## [63,] 180205560297
## [64,] 201276680407
## [65,] 228864089750
## [66,] 241221311920
## [67,] 253578534089
## [68,] 265935756259
## [69,] 278292978428
## [70,] 290650200598
## [71,] 303007422767
## [72,] 315364644936
## [73,] 327721867106
## [74,] 340079089275
## [75,] 352436311445
## [76,] 364793533614
## [77,] 377150755784
## [78,] 389507977953
ts.plot(datagabung[,1], xlab="Time Period", ylab="IPM", 
        main="DMA N=3 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.511480e+21
## MSE   1.025336e+20
## MAPE  1.083582e+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.765165e+22
## MSE   1.357819e+21
## MAPE  1.387076e+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 10.83 % dan Testing mendapat 13.61 %.

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.511480e+21                    3.830147e+21
## MSE   1.025336e+20                    7.510092e+19
## MAPE  1.083582e+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 10.83 %, 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.765165e+22                    3.329055e+21
## MSE   1.357819e+21                    3.329055e+21
## MAPE  1.387076e+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.87 % 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.

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 :

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

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

(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

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 koreksi autokorelasi dengan menggunakan metode-metode yang ada.Penanganan Autokrelasi

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 :

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

Didapat bahwa P-valuenya meningkat drasits 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

tambahin manualnya

rho<- modelCO$rho
rho
## [1] 0.9142712
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 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

Aman, Tidak terjadi autokorelasi.

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 nomral.

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.