Tugas Praktikum 2 MPDW
Hafidz Candra Ramadhan_G1401221079
2024-09-01
Pemanggilan Packages
library(readxl)## Warning: package 'readxl' was built under R version 4.3.2library(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, unionlibrary(TTR)## Warning: package 'TTR' was built under R version 4.3.2library(forecast)## Warning: package 'forecast' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(lmtest) #digunakan untuk uji formal pendeteksian autokorelasi## Warning: package 'lmtest' was built under R version 4.3.2## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(orcutt) #untuk membuat model regresi Cochrane-Orcutt## Warning: package 'orcutt' was built under R version 4.3.3library(HoRM) #untuk membuat model regresi Hildreth-Lu## Warning: package 'HoRM' was built under R version 4.3.3Input Data
Data yang digunakan dalam kesempatan kali ini adalah data Gross Domestic Product (GDP) pada South Asia periode 1960-2023.
#membaca data
southasia <- read_xlsx("D:\\IPB\\Statistika Sem 5\\MPDW (Deret Waktu)\\API_NY.GDP.MKTP.CD_DS2_en_excel_v2_3401502.xlsx", sheet = "Data")## New names:
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`
## • `` -> `...30`
## • `` -> `...31`
## • `` -> `...32`
## • `` -> `...33`
## • `` -> `...34`
## • `` -> `...35`
## • `` -> `...36`
## • `` -> `...37`
## • `` -> `...38`
## • `` -> `...39`
## • `` -> `...40`
## • `` -> `...41`
## • `` -> `...42`
## • `` -> `...43`
## • `` -> `...44`
## • `` -> `...45`
## • `` -> `...46`
## • `` -> `...47`
## • `` -> `...48`
## • `` -> `...49`
## • `` -> `...50`
## • `` -> `...51`
## • `` -> `...52`
## • `` -> `...53`
## • `` -> `...54`
## • `` -> `...55`
## • `` -> `...56`
## • `` -> `...57`
## • `` -> `...58`
## • `` -> `...59`
## • `` -> `...60`
## • `` -> `...61`
## • `` -> `...62`
## • `` -> `...63`
## • `` -> `...64`
## • `` -> `...65`
## • `` -> `...66`
## • `` -> `...67`
## • `` -> `...68`head(southasia)## # A tibble: 6 × 68
## `Data Source` World Development In…¹ ...3 ...4 ...5 ...6 ...7 ...8 ...9
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 Last Updated… 45471 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 3 Country Name Country Code Indi… Indi… 1960 1961 1962 1963 1964
## 4 Aruba ABW GDP … NY.G… <NA> <NA> <NA> <NA> <NA>
## 5 Africa Easte… AFE GDP … NY.G… 2121… 2230… 2370… 2577… 2804…
## 6 Afghanistan AFG GDP … NY.G… <NA> <NA> <NA> <NA> <NA>
## # ℹ abbreviated name: ¹`World Development Indicators`
## # ℹ 59 more variables: ...10 <chr>, ...11 <chr>, ...12 <chr>, ...13 <chr>,
## # ...14 <chr>, ...15 <chr>, ...16 <chr>, ...17 <chr>, ...18 <chr>,
## # ...19 <chr>, ...20 <chr>, ...21 <chr>, ...22 <chr>, ...23 <chr>,
## # ...24 <chr>, ...25 <chr>, ...26 <chr>, ...27 <chr>, ...28 <chr>,
## # ...29 <chr>, ...30 <chr>, ...31 <chr>, ...32 <chr>, ...33 <chr>,
## # ...34 <chr>, ...35 <chr>, ...36 <chr>, ...37 <chr>, ...38 <chr>, …#mengambil baris Country Name dan South Asia saja dan menghapus kolom 1,2,3,4
dtsas <- southasia[c(3,208), -c(1,2,3,4)]
#membuat data frame
dtsas <- as.data.frame(t(dtsas))
#mengubah nama kolom
colnames(dtsas) <- c("tahun", "gdp")
#mengubah tipe data
dtsas$tahun <- as.numeric(dtsas$tahun)
dtsas$gdp <- as.numeric(dtsas$gdp)
head(dtsas)## tahun gdp
## ...5 1960 47274328379
## ...6 1961 50467431650
## ...7 1962 53905757721
## ...8 1963 60496393964
## ...9 1964 69320094074
## ...10 1965 74300329652Eksplorasi Data
Sebelum melakukan regresi, akan diperlihatkan plot time-series dari GDP pada South Asia periode 1960-2023
#Membentuk objek time series
data.ts<-ts(dtsas$gdp)
data.ts## Time Series:
## Start = 1
## End = 64
## Frequency = 1
## [1] 4.727433e+10 5.046743e+10 5.390576e+10 6.049639e+10 6.932009e+10
## [6] 7.430033e+10 6.163414e+10 6.798891e+10 7.164250e+10 7.886021e+10
## [11] 8.514925e+10 9.060167e+10 9.133792e+10 1.045136e+11 1.265509e+11
## [16] 1.353672e+11 1.318736e+11 1.526997e+11 1.737963e+11 1.946675e+11
## [21] 2.355434e+11 2.500767e+11 2.587272e+11 2.738818e+11 2.726071e+11
## [26] 2.964355e+11 3.139552e+11 3.486019e+11 3.745842e+11 3.780169e+11
## [31] 4.069634e+11 3.620632e+11 3.845216e+11 3.808902e+11 4.320856e+11
## [36] 4.795855e+11 5.246813e+11 5.504088e+11 5.581171e+11 5.979953e+11
## [41] 6.476817e+11 6.625136e+11 6.954560e+11 8.116688e+11 9.416009e+11
## [46] 1.075831e+12 1.220740e+12 1.536037e+12 1.559694e+12 1.702715e+12
## [51] 2.082358e+12 2.294091e+12 2.327939e+12 2.389895e+12 2.614916e+12
## [56] 2.733619e+12 3.011037e+12 3.433921e+12 3.534067e+12 3.658210e+12
## [61] 3.493121e+12 4.079883e+12 4.327289e+12 4.475568e+12ts.plot(data.ts,xlab="Time Period ", 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:13
f = c()
for (i in t) {
f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}dt.gab <- cbind(aktual = c(data.ts,rep(NA,13)),
pemulusan1 = c(dt.sma,rep(NA,13)),
pemulusan2 = c(dt.dma, rep(NA,13)),
At = c(At, rep(NA,13)),
Bt = c(Bt,rep(NA,13)),
ramalan = c(dt.ramal, f[-1]))
dt.gab## aktual pemulusan1 pemulusan2 At Bt
## [1,] 4.727433e+10 NA NA NA NA
## [2,] 5.046743e+10 NA NA NA NA
## [3,] 5.390576e+10 5.054917e+10 NA NA NA
## [4,] 6.049639e+10 5.495653e+10 NA NA NA
## [5,] 6.932009e+10 6.124075e+10 7.255795e+10 6.689935e+10 5658598937
## [6,] 7.430033e+10 6.803894e+10 8.129267e+10 7.466581e+10 6626867365
## [7,] 6.163414e+10 6.841819e+10 7.345598e+10 7.093709e+10 2518896919
## [8,] 6.798891e+10 6.797446e+10 6.763566e+10 6.780506e+10 -169402384
## [9,] 7.164250e+10 6.708852e+10 6.561144e+10 6.634998e+10 -738537601
## [10,] 7.886021e+10 7.283054e+10 7.989594e+10 7.636324e+10 3532701241
## [11,] 8.514925e+10 7.855066e+10 9.000549e+10 8.427807e+10 5727417721
## [12,] 9.060167e+10 8.487038e+10 9.711009e+10 9.099023e+10 6119853549
## [13,] 9.133792e+10 8.902961e+10 9.878841e+10 9.390901e+10 4879396963
## [14,] 1.045136e+11 9.548438e+10 1.068636e+11 1.011740e+11 5689589047
## [15,] 1.265509e+11 1.074674e+11 1.277480e+11 1.176077e+11 10140298278
## [16,] 1.353672e+11 1.221439e+11 1.497012e+11 1.359225e+11 13778644946
## [17,] 1.318736e+11 1.312639e+11 1.532082e+11 1.422361e+11 10972157051
## [18,] 1.526997e+11 1.399802e+11 1.576819e+11 1.488310e+11 8850863091
## [19,] 1.737963e+11 1.527899e+11 1.756803e+11 1.642351e+11 11445229908
## [20,] 1.946675e+11 1.737212e+11 2.101694e+11 1.919453e+11 18224092086
## [21,] 2.355434e+11 2.013357e+11 2.521093e+11 2.267225e+11 25386798477
## [22,] 2.500767e+11 2.267625e+11 2.790746e+11 2.529186e+11 26156047171
## [23,] 2.587272e+11 2.481158e+11 2.935380e+11 2.708269e+11 22711099173
## [24,] 2.738818e+11 2.608952e+11 2.921700e+11 2.765326e+11 15637396727
## [25,] 2.726071e+11 2.684054e+11 2.869385e+11 2.776719e+11 9266570810
## [26,] 2.964355e+11 2.809748e+11 3.027407e+11 2.918578e+11 10882981961
## [27,] 3.139552e+11 2.943326e+11 3.205226e+11 3.074276e+11 13094995240
## [28,] 3.486019e+11 3.196642e+11 3.623448e+11 3.410045e+11 21340327715
## [29,] 3.745842e+11 3.457137e+11 3.973343e+11 3.715240e+11 25810255256
## [30,] 3.780169e+11 3.670677e+11 4.129060e+11 3.899868e+11 22919148335
## [31,] 4.069634e+11 3.865215e+11 4.266959e+11 4.066087e+11 20087196111
## [32,] 3.620632e+11 3.823478e+11 3.897522e+11 3.860500e+11 3702165825
## [33,] 3.845216e+11 3.845160e+11 3.846246e+11 3.845703e+11 54252292
## [34,] 3.808902e+11 3.758250e+11 3.656824e+11 3.707537e+11 -5071305085
## [35,] 4.320856e+11 3.991658e+11 4.244928e+11 4.118293e+11 12663499092
## [36,] 4.795855e+11 4.308537e+11 4.886649e+11 4.597593e+11 28905575911
## [37,] 5.246813e+11 4.787841e+11 5.638166e+11 5.213004e+11 42516243602
## [38,] 5.504088e+11 5.182252e+11 6.027669e+11 5.604961e+11 42270851911
## [39,] 5.581171e+11 5.444024e+11 6.055994e+11 5.750009e+11 30598490878
## [40,] 5.979953e+11 5.688404e+11 6.188758e+11 5.938581e+11 25017725328
## [41,] 6.476817e+11 6.012647e+11 6.607890e+11 6.310269e+11 29762185100
## [42,] 6.625136e+11 6.360635e+11 7.040781e+11 6.700708e+11 34007315158
## [43,] 6.954560e+11 6.685504e+11 7.350655e+11 7.018080e+11 33257552418
## [44,] 8.116688e+11 7.232128e+11 8.177539e+11 7.704833e+11 47270549930
## [45,] 9.416009e+11 8.162419e+11 9.767223e+11 8.964821e+11 80240187390
## [46,] 1.075831e+12 9.430337e+11 1.174109e+12 1.058571e+12 115537559663
## [47,] 1.220740e+12 1.079391e+12 1.345728e+12 1.212560e+12 133168729092
## [48,] 1.536037e+12 1.277536e+12 1.632635e+12 1.455085e+12 177549285533
## [49,] 1.559694e+12 1.438824e+12 1.785970e+12 1.612397e+12 173573381667
## [50,] 1.702715e+12 1.599482e+12 1.921218e+12 1.760350e+12 160867934624
## [51,] 2.082358e+12 1.781589e+12 2.131504e+12 1.956546e+12 174957390528
## [52,] 2.294091e+12 2.026388e+12 2.474191e+12 2.250290e+12 223901735410
## [53,] 2.327939e+12 2.234796e+12 2.675873e+12 2.455334e+12 220538344463
## [54,] 2.389895e+12 2.337308e+12 2.612930e+12 2.475119e+12 137810933857
## [55,] 2.614916e+12 2.444250e+12 2.655180e+12 2.549715e+12 105465159104
## [56,] 2.733619e+12 2.579477e+12 2.831073e+12 2.705275e+12 125798313438
## [57,] 3.011037e+12 2.786524e+12 3.152738e+12 2.969631e+12 183107117373
## [58,] 3.433921e+12 3.059526e+12 3.561560e+12 3.310543e+12 251017060163
## [59,] 3.534067e+12 3.326342e+12 3.864098e+12 3.595220e+12 268877952962
## [60,] 3.658210e+12 3.542066e+12 4.007576e+12 3.774821e+12 232754816938
## [61,] 3.493121e+12 3.561799e+12 3.731926e+12 3.646863e+12 85063506410
## [62,] 4.079883e+12 3.743738e+12 3.999478e+12 3.871608e+12 127870135376
## [63,] 4.327289e+12 3.966764e+12 4.385425e+12 4.176095e+12 209330481947
## [64,] 4.475568e+12 4.294246e+12 4.879574e+12 4.586910e+12 292663670557
## [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
## ramalan
## [1,] NA
## [2,] NA
## [3,] NA
## [4,] NA
## [5,] NA
## [6,] 7.255795e+10
## [7,] 8.129267e+10
## [8,] 7.345598e+10
## [9,] 6.763566e+10
## [10,] 6.561144e+10
## [11,] 7.989594e+10
## [12,] 9.000549e+10
## [13,] 9.711009e+10
## [14,] 9.878841e+10
## [15,] 1.068636e+11
## [16,] 1.277480e+11
## [17,] 1.497012e+11
## [18,] 1.532082e+11
## [19,] 1.576819e+11
## [20,] 1.756803e+11
## [21,] 2.101694e+11
## [22,] 2.521093e+11
## [23,] 2.790746e+11
## [24,] 2.935380e+11
## [25,] 2.921700e+11
## [26,] 2.869385e+11
## [27,] 3.027407e+11
## [28,] 3.205226e+11
## [29,] 3.623448e+11
## [30,] 3.973343e+11
## [31,] 4.129060e+11
## [32,] 4.266959e+11
## [33,] 3.897522e+11
## [34,] 3.846246e+11
## [35,] 3.656824e+11
## [36,] 4.244928e+11
## [37,] 4.886649e+11
## [38,] 5.638166e+11
## [39,] 6.027669e+11
## [40,] 6.055994e+11
## [41,] 6.188758e+11
## [42,] 6.607890e+11
## [43,] 7.040781e+11
## [44,] 7.350655e+11
## [45,] 8.177539e+11
## [46,] 9.767223e+11
## [47,] 1.174109e+12
## [48,] 1.345728e+12
## [49,] 1.632635e+12
## [50,] 1.785970e+12
## [51,] 1.921218e+12
## [52,] 2.131504e+12
## [53,] 2.474191e+12
## [54,] 2.675873e+12
## [55,] 2.612930e+12
## [56,] 2.655180e+12
## [57,] 2.831073e+12
## [58,] 3.152738e+12
## [59,] 3.561560e+12
## [60,] 3.864098e+12
## [61,] 4.007576e+12
## [62,] 3.731926e+12
## [63,] 3.999478e+12
## [64,] 4.385425e+12
## [65,] 4.879574e+12
## [66,] 5.172237e+12
## [67,] 5.464901e+12
## [68,] 5.757565e+12
## [69,] 6.050228e+12
## [70,] 6.342892e+12
## [71,] 6.635556e+12
## [72,] 6.928219e+12
## [73,] 7.220883e+12
## [74,] 7.513547e+12
## [75,] 7.806211e+12
## [76,] 8.098874e+12
## [77,] 8.391538e+12#Plot time series
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="GDP", main= "DMA N=3 Data GDP")
points(dt.gab[,1])
points(dt.gab[,3])
points(dt.gab[,6])
lines(dt.gab[,3],col="green",lwd=2)
lines(dt.gab[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"),
lty=8, col=c("blue","green","red"), cex=0.8)
Karena grafik di atas terlihat terpotong sebab data aktual dan data peramalan tidak sama panjang, maka akan dilakukan perbaikan dengan cara menambahkan ylim untuk menyesuaikan rentang sumbu y, sebagai berikut
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="GDP", main= "DMA N=3 Data GDP", ylim=c(0e+00,9e+12))
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("blue","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.209042e+23
## MSE 1.560855e+22
## MAPE 7.024005e+00Dari hasil diatas diperoleh nilai akurasi dengan metode Double Moving
Average (DMA) dengan parameter m=3, yaitu nilai MAPE kurang
dari 10%, dimana nilai ini dapat dikategorikan sebagai nilai akurasi
yang sangat baik. Nilai MAPE yang baik adalah dibawah 10%.
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<-dtsas[1:51,2]
testing<-dtsas[52:64,2]
#data time series
training.ts<-ts(training)
testing.ts<-ts(testing,start=52)
#eksplorasi data
plot(data.ts, col="red",main="Plot semua data")
points(data.ts)
plot(training.ts, col="blue",main="Plot data training")
points(training.ts)
plot(testing.ts, col="green",main="Plot data testing")
points(testing.ts)
Selanjutnya akan dilakukan pemulusan dengan DES, kali ini langsung dicari lambda dan gamma optimum sebagai berikut. Nilai lambda dan gamma optimum dapat dilihat pada smoothing parameters alpha untuk nilai lambda dan beta untuk nilai gamma.
#Lamda dan gamma optimum
des.opt<- HoltWinters(training.ts, gamma = FALSE)
des.opt## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = training.ts, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.539927
## beta : 0.752379
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 1.992062e+12
## b 2.210768e+11plot(des.opt)
legend("topleft", c("Data Aktual", "Peramalan"), col = c("black", "red"),
lty = c(1,1))
#ramalan
ramalandesopt<- forecast(des.opt, h=13)
ramalandesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 52 2.213139e+12 2.145389e+12 2.280889e+12 2.109524e+12 2.316754e+12
## 53 2.434216e+12 2.340946e+12 2.527486e+12 2.291572e+12 2.576860e+12
## 54 2.655293e+12 2.524547e+12 2.786038e+12 2.455335e+12 2.855250e+12
## 55 2.876370e+12 2.699479e+12 3.053260e+12 2.605839e+12 3.146901e+12
## 56 3.097446e+12 2.867659e+12 3.327234e+12 2.746017e+12 3.448876e+12
## 57 3.318523e+12 3.030175e+12 3.606872e+12 2.877533e+12 3.759514e+12
## 58 3.539600e+12 3.187701e+12 3.891499e+12 3.001417e+12 4.077784e+12
## 59 3.760677e+12 3.340698e+12 4.180656e+12 3.118375e+12 4.402979e+12
## 60 3.981754e+12 3.489507e+12 4.474000e+12 3.228927e+12 4.734580e+12
## 61 4.202831e+12 3.634392e+12 4.771269e+12 3.333479e+12 5.072182e+12
## 62 4.423907e+12 3.775569e+12 5.072246e+12 3.432359e+12 5.415456e+12
## 63 4.644984e+12 3.913217e+12 5.376752e+12 3.525842e+12 5.764126e+12
## 64 4.866061e+12 4.047489e+12 5.684633e+12 3.614163e+12 6.117959e+12Selanjutnya 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
## [1] NA NA 245222799 3410736579 5714625406 609667340mapedes.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 1.400215e+23
## MSE 2.745520e+21
## MAPE 6.429495e+00#Akurasi data testing
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt## Time Series:
## Start = 52
## End = 64
## Frequency = 1
## [1] -80951869718 106277058566 265397670148 261454046070 363827443131
## [6] 307486202619 105678762030 226610049754 323543976019 709709849417
## [11] 344024649799 317695260166 390493442765SSEtestingdesopt<-sum(selisihdesopt^2)
SSEtestingdesopt<-SSEtestingdesopt/length(testing.ts)
MAPEtestingdesopt<-sum(abs(selisihdesopt/testing.ts)*100)/length(testing.ts)
akurasiDesTesting <- matrix(c(SSEtestingdesopt,SSEtestingdesopt,MAPEtestingdesopt))
row.names(akurasiDesTesting)<- c("SSE", "MSE", "MAPE")
colnames(akurasiDesTesting) <- c("Akurasi lamda dan gamma optimum")
akurasiDesTesting## Akurasi lamda dan gamma optimum
## SSE 1.09709e+23
## MSE 1.09709e+23
## MAPE 8.91301e+00Setelah 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.209042e+23 1.400215e+23
## MSE 1.560855e+22 2.745520e+21
## MAPE 7.024005e+00 6.429495e+00Berdasarkan perbandingan akurasi tersebut, terlihat nilai SSE dan MAPE metode DES lebih kecil dibandingkan dengan metode DMA, tetapi pada nilai MSE metode DES lebih besar 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
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.2ggplot(dtsas, aes(x=tahun, y=gdp)) +
geom_point() +
labs(title = "Scatter Plot Tahun vs Nilai GDP",
x = "Tahun",
y = "Nilai GDP")
#Menampilkan Nilai Korelasi
cor(dtsas$tahun, dtsas$gdp)## [1] 0.8549233Berdasarkan scatter plot di atas, terlihat adanya hubungan / korelasi
positif antara peubah tahun dengan nilai IPM, 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.8549233\).
Setalah mengetahui adanya hubungan antar dua peubah, maka model regresi dapat ditentukan.
Regresi
#Pembuatan Model Regresi
#model regresi
model<- lm(gdp~tahun, data = dtsas)
summary(model)##
## Call:
## lm(formula = gdp ~ tahun, data = dtsas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.267e+11 -5.461e+11 -4.045e+10 4.315e+11 1.638e+12
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.143e+14 8.886e+12 -12.86 <2e-16 ***
## tahun 5.790e+10 4.462e+09 12.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.594e+11 on 62 degrees of freedom
## Multiple R-squared: 0.7309, Adjusted R-squared: 0.7266
## F-statistic: 168.4 on 1 and 62 DF, p-value: < 2.2e-16Model yang dihasilkan adalah \[y_i=-1.143+5.790x_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.7309\). Artinya, sebesar 73.09% 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 di atas menunjukkan bahwa sisaan cenderung menyebar normal, tetapi histogram dari sisaan tidak menunjukkan demikian. Selanjutnya, dua plot di samping kanan digunakan untuk melihat autokorelasi. Plot Sisaan vs Fitted Value dan Plot Sisaan vs Order menunjukkan adanya pola pada sisaan. Untuk lebih lanjut akan digunakan uji formal melihat normalitas sisaan dan plot ACF dan PACF untuk melihat apakah ada autokorelasi atau tidak.
#Melihat Sisaan Menyebar Normal/Tidak
#H0: sisaan mengikuti sebaran normal
#H1: sisaan tidak mengikuti sebaran normal
shapiro.test(sisaan)##
## Shapiro-Wilk normality test
##
## data: sisaan
## W = 0.95376, p-value = 0.0176ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.087069, p-value = 0.6842
## alternative hypothesis: two-sidedBerdasarkan uji formal Saphiro-Wilk didapatkan p-value < \(\alpha\) (5%). Artinya, tidak cukup bukti untuk menyatakan sisaan berdistribusi normal. Akan tetapi, pada uji formal Kolmogorov-Smirnov didapatkan nilai p-value > \(\alpha\) (5%). Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal. Namun, karena uji formal Kolmogorov-Smirnov lebih cocok terhadap data yang besar, maka akan lebih dipercaya hasil dari uji Kolmogorov-Smirnov. Selanjutnya akan dilakukan uji autokorelasi.
#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)
Berdasarkan plot ACF dan PACF, terlihat hampir sebagian berada diluar rentang batas dan 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.033641, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0Berdasarkan hasil DW Test, didapatkan nilai \(DW = 0.033641\) dan p-value < \(\alpha\) (5%). Berdasarkan tabel Durbin-Watson diperoleh nilai \(DL = 1.5635\) dan \(DU = 1.6268\). Nilai DW < DL. Artinya, terdapat autokorelasi positif pada pemodelan regresi. Hal ini diperkuat dengan nilai p-value < 0.05 sehingga dapat disimpulkan bahwa tolak H0, cukup bukti mengatakan adanya autokorelasi. Oleh karena itu, diperlukan penangan autokorelasi. Penanganan yang akan digunakan menggunakan dua metode, yaitu Cochrane-Orcutt dan Hildret-Lu.
Penanganan Autokorelasi
Metode Cochrane-Orcutt
Penanganan metode Cochrane-Orcutt dapat dilakukan dengan bantuan
packages Orcutt pada aplikasi R maupun secara manual.
Berikut ini ditampilkan cara menggunakan bantuan library
packages Orcutt.
#Penanganan Autokorelasi Cochrane-Orcutt
modelCO<-cochrane.orcutt(model, max.iter = 1000)
modelCO## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = gdp ~ tahun, data = dtsas)
##
## number of interaction: 765
## rho 1.03077
##
## Durbin-Watson statistic
## (original): 0.03364 , p-value: 6.338e-54
## (transformed): 2.24555 , p-value: 8.026e-01
##
## coefficients:
## (Intercept) tahun
## 1.340921e+14 -6.914329e+10Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=(1.340921e+14)-(6.914329e+10)x_t\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(2.24555\) dan \(8.026e-01\). Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.6268 < DW < 2.3732\). 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 \(1.03077\). Nilai tersebut dapat diketahui dengan syntax berikut.
#Rho optimum
rho<- modelCO$rho
rho## [1] 1.03077Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.
gdp <- dtsas$gdp
gdp## [1] 4.727433e+10 5.046743e+10 5.390576e+10 6.049639e+10 6.932009e+10
## [6] 7.430033e+10 6.163414e+10 6.798891e+10 7.164250e+10 7.886021e+10
## [11] 8.514925e+10 9.060167e+10 9.133792e+10 1.045136e+11 1.265509e+11
## [16] 1.353672e+11 1.318736e+11 1.526997e+11 1.737963e+11 1.946675e+11
## [21] 2.355434e+11 2.500767e+11 2.587272e+11 2.738818e+11 2.726071e+11
## [26] 2.964355e+11 3.139552e+11 3.486019e+11 3.745842e+11 3.780169e+11
## [31] 4.069634e+11 3.620632e+11 3.845216e+11 3.808902e+11 4.320856e+11
## [36] 4.795855e+11 5.246813e+11 5.504088e+11 5.581171e+11 5.979953e+11
## [41] 6.476817e+11 6.625136e+11 6.954560e+11 8.116688e+11 9.416009e+11
## [46] 1.075831e+12 1.220740e+12 1.536037e+12 1.559694e+12 1.702715e+12
## [51] 2.082358e+12 2.294091e+12 2.327939e+12 2.389895e+12 2.614916e+12
## [56] 2.733619e+12 3.011037e+12 3.433921e+12 3.534067e+12 3.658210e+12
## [61] 3.493121e+12 4.079883e+12 4.327289e+12 4.475568e+12gdp[-1]## [1] 5.046743e+10 5.390576e+10 6.049639e+10 6.932009e+10 7.430033e+10
## [6] 6.163414e+10 6.798891e+10 7.164250e+10 7.886021e+10 8.514925e+10
## [11] 9.060167e+10 9.133792e+10 1.045136e+11 1.265509e+11 1.353672e+11
## [16] 1.318736e+11 1.526997e+11 1.737963e+11 1.946675e+11 2.355434e+11
## [21] 2.500767e+11 2.587272e+11 2.738818e+11 2.726071e+11 2.964355e+11
## [26] 3.139552e+11 3.486019e+11 3.745842e+11 3.780169e+11 4.069634e+11
## [31] 3.620632e+11 3.845216e+11 3.808902e+11 4.320856e+11 4.795855e+11
## [36] 5.246813e+11 5.504088e+11 5.581171e+11 5.979953e+11 6.476817e+11
## [41] 6.625136e+11 6.954560e+11 8.116688e+11 9.416009e+11 1.075831e+12
## [46] 1.220740e+12 1.536037e+12 1.559694e+12 1.702715e+12 2.082358e+12
## [51] 2.294091e+12 2.327939e+12 2.389895e+12 2.614916e+12 2.733619e+12
## [56] 3.011037e+12 3.433921e+12 3.534067e+12 3.658210e+12 3.493121e+12
## [61] 4.079883e+12 4.327289e+12 4.475568e+12#Transformasi Manual
tahun <- dtsas$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
## -3.780e+11 -3.321e+10 -4.101e+08 1.958e+10 3.768e+11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.126e+12 1.310e+12 -3.150 0.00253 **
## tahun.trans -6.914e+10 2.174e+10 -3.181 0.00231 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.653e+10 on 61 degrees of freedom
## Multiple R-squared: 0.1423, Adjusted R-squared: 0.1282
## F-statistic: 10.12 on 1 and 61 DF, p-value: 0.002309Hasil 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)
## 1.340921e+14b1## tahun.trans
## -69143293834Hasil 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,2.0, by= 0.05))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4)## rho SSE
## 1 0.10 2.144548e+25
## 2 0.15 1.926274e+25
## 3 0.20 1.720049e+25
## 4 0.25 1.525873e+25
## 5 0.30 1.343747e+25
## 6 0.35 1.173669e+25
## 7 0.40 1.015640e+25
## 8 0.45 8.696609e+24
## 9 0.50 7.357305e+24
## 10 0.55 6.138492e+24
## 11 0.60 5.040170e+24
## 12 0.65 4.062338e+24
## 13 0.70 3.204998e+24
## 14 0.75 2.468148e+24
## 15 0.80 1.851789e+24
## 16 0.85 1.355921e+24
## 17 0.90 9.805437e+23
## 18 0.95 7.256571e+23
## 19 1.00 8.970918e+23
## 20 1.05 5.773566e+23
## 21 1.10 6.839425e+23
## 22 1.15 9.110192e+23
## 23 1.20 1.258587e+24
## 24 1.25 1.726645e+24
## 25 1.30 2.315194e+24
## 26 1.35 3.024234e+24
## 27 1.40 3.853765e+24
## 28 1.45 4.803787e+24
## 29 1.50 5.874299e+24
## 30 1.55 7.065302e+24
## 31 1.60 8.376796e+24
## 32 1.65 9.808781e+24
## 33 1.70 1.136126e+25
## 34 1.75 1.303422e+25
## 35 1.80 1.482768e+25
## 36 1.85 1.674163e+25
## 37 1.90 1.877607e+25
## 38 1.95 2.093100e+25
## 39 2.00 2.320642e+25Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 1.05. 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.05, kali ini jarak antar \(ρ\) adalah 0.01 dan dilakukan pada selang 0.09 sampai dengan 1.20.
#Rho optimal di sekitar 0.4
rOpt <- seq(0.09, 1.20, by= 0.01)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model))}))
#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=1.03, y=1.6e+24, labels = "rho=1.03", cex = 0.9)
Perhitungan yang dilakukan aplikasi R menunjukkan bahwa
nilai \(ρ\) optimum, yaitu saat SSE
terkecil terdapat pada nilai \(ρ=1.03\). 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(1.03, model)
summary(modelHL)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.771e+11 -3.360e+10 -9.046e+07 1.977e+10 3.775e+11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.207e+12 1.309e+12 -3.213 0.00210 **
## x -7.234e+10 2.229e+10 -3.245 0.00191 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.654e+10 on 61 degrees of freedom
## Multiple R-squared: 0.1472, Adjusted R-squared: 0.1332
## F-statistic: 10.53 on 1 and 61 DF, p-value: 0.001911#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-1.03), coef(modelHL)[2],"x", sep = "")## y = 1.402288e+14-72337409559xSetelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=(1.402288e+14)-(72337409559)x_t\]
#Deteksi autokorelasi
dwtest(modelHL)##
## Durbin-Watson test
##
## data: modelHL
## DW = 2.2437, p-value = 0.8006
## alternative hypothesis: true autocorrelation is greater than 0Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(2.2437\) berada pada selang daerah tidak ada autokorelasi, yaitu pada rentang DU < DW < 4-DU atau \(1.6268 < DW < 2.3732\). Hal tersebut juga didukung oleh p-value sebesar \(0.8006\), di mana p-value > \(\alpha\)=5%. Artinya tak tolak \(H_0\) atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai GDP dengan metode Hildreth-Lu pada taraf nyata 5%.
Terakhir, akan dibandingkan nilai SSE dari ketiga metode (metode awal, metode Cochrane-Orcutt, dan Hildreth-Lu).
#Perbandingan
sseModelawal <- anova(model)$`Sum Sq`[-1]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-1]
sseModelHL <- anova(modelHL)$`Sum Sq`[-1]
mseModelawal <- sseModelawal/length(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 2.695421e+25 5.684453e+23 5.684596e+23
## MSE 4.211595e+23 8.881958e+21 8.882182e+21Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki SSE berturut-turut \(5.684453e+23\) dan \(5.927685e+23\). Nilai tersebut lebih baik dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar \(2.695421e+25\).
Simpulan
Autokorelasi yang terdapat pada data GDP terjadi akibat adanya korelasi di antara unsur penyusunnya. Indikator GDP yang erat hubungannya dengan perekonomian sangat rawan menjadi penyebab adanya autokorelasi. Adanya autokorelasi menyebabkan model regresi kurang baik karena akan meingkatkan galatnya. Autokorelasi dapat dideteksi secara eksploratif melalui plot sisaan, ACF, dan PACF, serta dengan uji formal Durbin-Watson. Namun, autokorelasi tersebut dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode menghasilkan nilai SSE yang lebih baik dibandingkan model awal, artinya keduanya baik untuk digunakan.
Daftar Pustaka
Aprianto A, Debataraja NN, Imro’ah N. 2020. Metode cochrane-orcutt untuk mengatasi autokorelasi pada estimasi parameter ordinary least squares. Bimaster : Buletin Ilmiah Matematika, Statistika dan Terapannya. 9(1):95–102. doi:10.26418/bbimst.v9i1.38590.
BPS. 2021a. Indeks Pembangunan Manusia 2020. Jakarta (ID): Badan Pusat Statistik.
BPS. 2021b. Indeks Pembangunan Manusia (IPM) 2021. Berita Resmi Statistik., siap terbit.