Analisis Deret Waktu Data Harga Emas 2010-2022

Kelompok 1

Paralel 2 - MPDW 2022

Anggota:

# Berliana Apriyanti
# Kenia Mauilidia
# Mohammad Abror Gustiansyah
# Oksi Al Hadi
# Raffael Julio Roger

A. Analisis Data

Library

library(googlesheets4)
library(forecast)
library(TTR)
library(imputeTS)
library(tseries)
library(ggplot2)
library(dplyr)
library(graphics)
library(TSA)
library(tidyverse)
library(lubridate)
library(gridExtra)
library(ggfortify)
library(cowplot)
library(graphics)
library(lmtest)
library(stats)
library(MASS)
library(googlesheets4)
library(tseries)
library(forecast)
library(TTR)
library(TSA)
library(imputeTS)
library(ggplot2)
library(lmtest)
library(FinTS)

Import Data

dataemas <- read.csv("D:/Kuliah/Semester 5/MPDW/Data Emas Fixx.csv")
ts.dataemas <- ts(dataemas)
kableExtra::kable(head(dataemas) ,caption = 'Data Rata-Rata Harga Emas Global 2010-2022')

Data Rata-Rata Harga Emas Global 2010-2022
Waktu	USD
1/29/2010	1078.50
2/26/2010	1108.25
3/31/2010	1115.50
4/30/2010	1179.25
5/31/2010	1207.50
6/30/2010	1244.00

str(dataemas)

## 'data.frame':    146 obs. of  2 variables:
##  $ Waktu: chr  "1/29/2010" "2/26/2010" "3/31/2010" "4/30/2010" ...
##  $ USD  : num  1078 1108 1116 1179 1208 ...

Plot Time Series

Splitting Data

y <- dataemas$USD
training<-y[1:117]
testing<-y[118:146]

training.ts<-ts(training)
testing.ts<-ts(testing,start=118)

plot(training.ts, col="blue",main="Plot Data Training")
points(training.ts)

plot(testing.ts, col="blue",main="Plot data Testing")
points(testing.ts)

B. Pemulusan Data

DMA

Metode pemulusan Double Moving Average (DMA) merupakan metode yang tepat untuk data penjualan emas karena data tersebut memiliki pola data tren. Fungsi di bawah ini akan me-looping nilai n dari 2 sampai dengan 11 sehingga didapatkan n optimum.

akurasi = data.frame()
for(n in 2:11){
  data.sma<-SMA(training.ts, n=n)
  dma <- SMA(data.sma, n=n)
  At <- 2*data.sma - dma
  Bt <- 2/(n-1)*(data.sma - dma)
  data.ramal<- c(NA, At+Bt)
  
  error.dma = training.ts-data.ramal[1:length(training.ts)]
  SSE.dma = sum(error.dma[(2*n):length(training.ts)]^2, na.rm = T)
  RMSE.dma = sqrt(mean(error.dma[(2*n):length(training.ts)]^2, na.rm = T))
  MAPE.dma = mean(abs((error.dma[(2*n):length(training.ts)]/training.ts[(2*n):length(training.ts)])*100), na.rm = T)
  
  akurasi.dma <- matrix(c(SSE.dma, RMSE.dma, MAPE.dma))
  row.names(akurasi.dma)<- c("SSE", "RMSE", "MAPE")
  colnames(akurasi.dma) <- c(paste("Akurasi m =", n))
  akurasi = rbind(akurasi, c(n=n, RMSE=akurasi.dma[2], MAPE=akurasi.dma[3]))
}
colnames(akurasi) = c("n", "RMSE", "MAPE")
akurasi[order(akurasi[,2]),]

##     n      RMSE     MAPE
## 2   3  83.51942 4.938363
## 4   5  84.28369 4.787540
## 3   4  85.47429 4.841119
## 5   6  85.50707 4.737835
## 1   2  89.76499 4.732633
## 6   7  91.45361 5.097534
## 7   8  97.75974 5.412382
## 8   9 103.57142 5.662378
## 9  10 107.28141 5.891334
## 10 11 107.80972 6.071765

Berdasarkan tabel di atas, nilai RMSE dan MAPE terkecil dimiliki oleh metode pemulusan DMA pada n = 3 yang dapat dijadikan n optimum dalam rentang nilai n dari 2 sampai 11. Dengan nilai MSE terkecil sebesar 83.29561 dan MAPE terkecil sebesar 4.918893 . Syntax di bawah ini akan memodelkan pemulusan DMA dengan n optimum.

Pemulusan Double Moving Average dengan N Optimum (3)

n=3
data.sma<-SMA(training.ts, n=n)
dma <- SMA(data.sma, n=n)
At <- 2*data.sma - dma
Bt <- 2/(n-1)*(data.sma - dma)
data.ramal<- c(NA, At+Bt)

fcast = c()
for (i in 1:20) {
  fcast[i] = At[length(At)] + Bt[length(Bt)]*(i)
} # Ramalan DMA

data.gab <- cbind(aktual = c(training.ts,rep(NA,21)), pemulusan = c(dma, rep(NA,21)), At = c(At, rep(NA,21)), Bt = c(Bt,rep(NA,21)), ramalan = c(data.ramal, fcast[-1]))

## Warning in cbind(aktual = c(training.ts, rep(NA, 21)), pemulusan = c(dma, :
## number of rows of result is not a multiple of vector length (arg 5)

tail(data.gab)

##        aktual pemulusan At Bt  ramalan
## [133,]     NA        NA NA NA 2208.489
## [134,]     NA        NA NA NA 2251.317
## [135,]     NA        NA NA NA 2294.144
## [136,]     NA        NA NA NA 2336.972
## [137,]     NA        NA NA NA 2379.800
## [138,]     NA        NA NA NA       NA

Plot DMA

# Plot
plot.ts(dataemas$USD, xlab="Waktu", ylab="USD", main="Plot Pemulusan dan Peramalan dengan Metode DMA N Optimum (3)", ylim=c(1000,2500), type="l")
points(data.gab[,1])
lines(data.gab[,2],col="green",lwd=2)
lines(data.gab[,5],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan DMA","data peramalan"), lty=1, col=c("black","green","red"), cex=0.8)

Plot di atas menunjukkan hasil pemulusan dan peramalan metode DMA pada n = 3. Dapat dilihat bahwa hasil pemulusan, dapat dilihat pada garis pemulusan yang masih cukup menggambarkan pola data penjualan emas aktual. Kemudian, hasil peramalan menunjukkan pola naik yang nilainya berbeda jauh dengan data aktual (data testing).

Nilai Akurasi DMA

RMSE_DMA_train <- sqrt(mean((training.ts-data.ramal[-88])^2, na.rm = T))
RMSE_DMA_test <- sqrt(mean((testing.ts-fcast)^2))

## Warning in `-.default`(testing.ts, fcast): longer object length is not a
## multiple of shorter object length

MAPE_DMA_train <- mean(abs((training.ts-data.ramal[-88])/training.ts)*100, na.rm = T)
MAPE_DMA_test <- mean(abs((testing.ts-fcast)/testing.ts)*100)

## Warning in `-.default`(testing.ts, fcast): longer object length is not a
## multiple of shorter object length

Akurasi_DMA <- data.frame("DMA_Train"=c(RMSE_DMA_train,MAPE_DMA_train), "DMA_Test"=c(RMSE_DMA_test,MAPE_DMA_test))
rownames(Akurasi_DMA) <- c("RMSE","MAPE")
Akurasi_DMA

##      DMA_Train  DMA_Test
## RMSE 80.810628 252.05395
## MAPE  4.586484  10.73407

Dari syntax di atas, didapatkan bahwa nilai RMSE DMA testing sebesar 503.43018 dan MAPE DMA testing sebesar 23.17036. Kedua nilai error ini cukup besar jika dibandingkan data trainingnya.

Double Exponential Smoothing (DES)

Selain DMA, Metode pemulusan lain yang cocok untuk data yang memiliki pola data tren adalah DES. Metode ini memandang bahwa semakin besar lag, maka semakin tidak berpengaruh data tersebut terhadap data aktual saat ini. Fungsi di bawah ini akan me-looping parameter DES, yaitu alpha dan beta (dari 0.1 hingga 0.9) sehingga didapatkan parameter optimum.

#Pengulangan Melihat Nilai Keakuratan
a = seq(.1,.9,.1)
b = seq(.1, .9,.1)
output = data.frame()
for (i in a) {
  for (j in b) {
    des <- HoltWinters(training.ts, alpha = i, beta=j, gamma=F)
    sse <- des$SSE
    rmse <- sqrt(sse/length(training.ts))
    mape <- mean(abs((des$fitted[,1]-training.ts)/training.ts)*100)
    akurasi <- cbind("SSE"=sse, "RMSE"=rmse, "MAPE"=mape)
    output <- rbind(output, akurasi)
  }
}
output <- cbind("Alpha" = rep(a,each=9), "Beta"=b, output)
dplyr::arrange(.data=output,RMSE) #Mengurutkan data dari nilai SSE/MSE terkecil

##    Alpha Beta       SSE      RMSE     MAPE
## 1    0.8  0.1  539003.5  67.87391 3.619096
## 2    0.7  0.1  540056.5  67.94017 3.634656
## 3    0.9  0.1  548572.9  68.47377 3.635503
## 4    0.6  0.1  553627.6  68.78852 3.667440
## 5    0.7  0.2  560533.4  69.21621 3.746367
## 6    0.8  0.2  565512.7  69.52296 3.731267
## 7    0.6  0.2  568092.9  69.68138 3.809709
## 8    0.9  0.2  581832.0  70.51896 3.751041
## 9    0.5  0.1  584581.7  70.68539 3.702540
## 10   0.7  0.3  587644.1  70.87030 3.870450
## 11   0.6  0.3  590248.5  71.02717 3.932039
## 12   0.5  0.2  592576.2  71.16708 3.880367
## 13   0.8  0.3  598118.7  71.49913 3.858005
## 14   0.5  0.3  610657.3  72.24468 4.019681
## 15   0.6  0.4  614374.3  72.46421 4.037684
## 16   0.7  0.4  617473.8  72.64678 3.971655
## 17   0.9  0.3  621224.9  72.86710 3.886459
## 18   0.5  0.4  629016.2  73.32262 4.113274
## 19   0.8  0.4  633834.3  73.60290 3.990484
## 20   0.6  0.5  641036.3  74.01988 4.117440
## 21   0.4  0.2  643400.5  74.15625 3.965442
## 22   0.4  0.1  643543.3  74.16448 3.819144
## 23   0.5  0.5  648303.5  74.43826 4.194941
## 24   0.7  0.5  650106.0  74.54168 4.104403
## 25   0.4  0.3  660499.6  75.13518 4.168900
## 26   0.9  0.4  664514.7  75.36321 4.020308
## 27   0.5  0.6  670253.8  75.68794 4.274529
## 28   0.6  0.6  670917.0  75.72538 4.215837
## 29   0.8  0.5  672462.3  75.81254 4.129900
## 30   0.4  0.4  675741.3  75.99715 4.264309
## 31   0.7  0.6  685454.4  76.54139 4.259563
## 32   0.4  0.5  687901.3  76.67789 4.327123
## 33   0.5  0.7  695871.0  77.12078 4.356090
## 34   0.4  0.6  699491.7  77.32116 4.376990
## 35   0.6  0.7  704024.8  77.57130 4.358606
## 36   0.9  0.5  711691.0  77.99249 4.145860
## 37   0.4  0.7  713233.9  78.07699 4.424980
## 38   0.8  0.6  713732.4  78.10427 4.270901
## 39   0.7  0.7  722935.5  78.60621 4.407928
## 40   0.5  0.8  725458.4  78.74325 4.452544
## 41   0.4  0.8  730523.3  79.01765 4.483542
## 42   0.3  0.2  737446.8  79.39121 4.205519
## 43   0.6  0.8  739807.5  79.51818 4.519839
## 44   0.4  0.9  751777.5  80.15890 4.555915
## 45   0.3  0.1  755609.9  80.36295 4.168468
## 46   0.8  0.7  757275.2  80.45146 4.373580
## 47   0.5  0.9  758923.0  80.53894 4.562420
## 48   0.3  0.3  760150.6  80.60405 4.355843
## 49   0.7  0.8  761702.5  80.68629 4.531071
## 50   0.9  0.6  762954.4  80.75257 4.256133
## 51   0.6  0.9  777304.8  81.50847 4.668990
## 52   0.3  0.4  782560.0  81.78354 4.566110
## 53   0.3  0.5  797882.0  82.58029 4.663659
## 54   0.7  0.9  800942.1  82.73850 4.622105
## 55   0.8  0.8  802935.6  82.84140 4.449691
## 56   0.3  0.6  805087.4  82.95233 4.681934
## 57   0.3  0.8  805966.1  82.99759 4.737793
## 58   0.3  0.7  806091.7  83.00405 4.677841
## 59   0.3  0.9  808743.5  83.14047 4.759417
## 60   0.9  0.7  818800.9  83.65583 4.365983
## 61   0.8  0.9  851031.6  85.28643 4.505795
## 62   0.9  0.8  880237.1  86.73750 4.471834
## 63   0.2  0.2  910372.7  88.20977 4.734341
## 64   0.2  0.3  931617.5  89.23308 4.938248
## 65   0.9  0.9  948853.7  90.05477 4.608772
## 66   0.2  0.4  984065.0  91.71048 5.103562
## 67   0.2  0.1 1010835.5  92.94955 4.921396
## 68   0.2  0.5 1027445.9  93.71013 5.249791
## 69   0.2  0.6 1058656.0  95.12277 5.497555
## 70   0.2  0.7 1086746.4  96.37650 5.602457
## 71   0.2  0.8 1106589.5  97.25240 5.638832
## 72   0.2  0.9 1107906.6  97.31026 5.562416
## 73   0.1  0.4 1277729.4 104.50246 5.947659
## 74   0.1  0.5 1294114.1 105.17036 6.175418
## 75   0.1  0.3 1314939.0 106.01318 5.978746
## 76   0.1  0.6 1381545.1 108.66498 6.496542
## 77   0.1  0.7 1532968.3 114.46525 6.872032
## 78   0.1  0.2 1537704.9 114.64195 6.555855
## 79   0.1  0.8 1694500.5 120.34496 7.299353
## 80   0.1  0.9 1803994.1 124.17227 7.539577
## 81   0.1  0.1 2198413.0 137.07609 7.321322

Setelah diurutkan berdasarkan nilai RMSE, dapat terlihat bahwa nilai alpha optimum 0.8 dan beta optimum 0.1, dengan nilai RMSE terkecilnya yaitu sebesar 68.03389. Kedua parameter ini yang kemudian akan digunakan untuk pemulusan.

Pemulusan Double Exponential Smoothing dengan alpha dan beta optimum (0.8, 0.1)

df_des <- HoltWinters(training.ts, alpha = .8, beta=.1, gamma=F)
datades <- data.frame(training.ts, c(NA, NA, df_des$fitted[,1]))
colnames(datades) = c("y_train","y_train_hat")
datades

##     y_train y_train_hat
## 1   1078.50          NA
## 2   1108.25          NA
## 3   1115.50    1138.000
## 4   1179.25    1147.950
## 5   1207.50    1203.444
## 6   1244.00    1237.467
## 7   1169.00    1273.995
## 8   1246.00    1212.900
## 9   1307.00    1264.930
## 10  1346.75    1327.501
## 11  1383.50    1373.355
## 12  1405.50    1412.738
## 13  1327.00    1437.635
## 14  1411.00    1370.964
## 15  1439.00    1428.032
## 16  1535.50    1462.724
## 17  1536.50    1552.684
## 18  1505.50    1570.181
## 19  1628.50    1543.706
## 20  1813.50    1643.595
## 21  1620.00    1825.165
## 22  1722.00    1690.266
## 23  1746.00    1747.425
## 24  1531.00    1777.942
## 25  1744.00    1592.291
## 26  1770.00    1737.697
## 27  1662.50    1790.162
## 28  1651.25    1704.443
## 29  1558.00    1674.043
## 30  1598.50    1584.080
## 31  1622.00    1599.641
## 32  1648.50    1623.342
## 33  1776.00    1651.295
## 34  1719.00    1768.862
## 35  1726.00    1742.786
## 36  1657.50    1741.828
## 37  1664.75    1680.090
## 38  1588.50    1672.315
## 39  1598.25    1603.055
## 40  1469.00    1596.619
## 41  1394.50    1481.722
## 42  1192.00    1392.165
## 43  1314.50    1196.240
## 44  1394.75    1264.516
## 45  1326.50    1352.790
## 46  1324.00    1313.742
## 47  1253.00    1304.753
## 48  1204.50    1242.015
## 49  1251.00    1187.666
## 50  1326.50    1219.063
## 51  1291.75    1294.337
## 52  1288.50    1281.385
## 53  1250.50    1276.764
## 54  1315.00    1243.338
## 55  1285.25    1293.986
## 56  1285.75    1279.617
## 57  1216.50    1277.634
## 58  1164.25    1216.946
## 59  1182.75    1158.793
## 60  1206.00    1163.879
## 61  1260.25    1186.866
## 62  1214.00    1240.734
## 63  1187.00    1212.369
## 64  1180.25    1183.066
## 65  1191.40    1171.581
## 66  1171.00    1179.789
## 67  1098.40    1164.408
## 68  1135.00    1097.971
## 69  1114.00    1116.926
## 70  1142.35    1103.683
## 71  1061.90    1126.807
## 72  1060.00    1061.880
## 73  1111.80    1047.224
## 74  1234.90    1090.899
## 75  1237.00    1209.634
## 76  1285.65    1237.250
## 77  1212.10    1285.565
## 78  1320.75    1230.511
## 79  1342.00    1313.639
## 80  1309.25    1349.534
## 81  1322.50    1327.290
## 82  1272.00    1333.058
## 83  1178.10    1288.927
## 84  1145.90    1196.115
## 85  1212.80    1147.775
## 86  1255.60    1196.829
## 87  1244.85    1245.582
## 88  1266.45    1246.674
## 89  1266.20    1265.754
## 90  1242.25    1269.406
## 91  1267.55    1248.804
## 92  1311.75    1266.423
## 93  1283.10    1308.933
## 94  1270.15    1292.448
## 95  1280.20    1277.008
## 96  1291.00    1282.215
## 97  1345.05    1292.599
## 98  1317.85    1342.112
## 99  1323.85    1328.314
## 100 1313.20    1329.997
## 101 1305.35    1320.470
## 102 1250.45    1311.075
## 103 1220.95    1260.426
## 104 1202.45    1223.538
## 105 1187.25    1199.673
## 106 1214.95    1181.747
## 107 1217.55    1202.977
## 108 1279.00    1210.469
## 109 1323.25    1266.610
## 110 1319.15    1317.770
## 111 1295.40    1324.832
## 112 1282.30    1304.890
## 113 1295.55    1288.614
## 114 1409.00    1296.514
## 115 1427.55    1397.853
## 116 1528.40    1435.336
## 117 1485.30    1530.958

ramal_des <- forecast::forecast(df_des,h=21)
data.frame(ramal_des)

##     Point.Forecast    Lo.80    Hi.80     Lo.95    Hi.95
## 118       1511.950 1423.845 1600.054 1377.2058 1646.694
## 119       1529.468 1412.107 1646.829 1349.9801 1708.956
## 120       1546.986 1402.323 1691.649 1325.7433 1768.229
## 121       1564.505 1393.264 1735.745 1302.6154 1826.394
## 122       1582.023 1384.386 1779.660 1279.7633 1884.282
## 123       1599.541 1375.402 1823.680 1256.7505 1942.331
## 124       1617.059 1366.149 1867.970 1233.3248 2000.794
## 125       1634.577 1356.525 1912.630 1209.3324 2059.822
## 126       1652.096 1346.466 1957.725 1184.6756 2119.516
## 127       1669.614 1335.932 2003.296 1159.2914 2179.936
## 128       1687.132 1324.896 2049.369 1133.1391 2241.125
## 129       1704.650 1313.340 2095.961 1106.1926 2303.108
## 130       1722.169 1301.254 2143.083 1078.4361 2365.901
## 131       1739.687 1288.634 2190.740 1049.8605 2429.513
## 132       1757.205 1275.475 2238.935 1020.4619 2493.948
## 133       1774.723 1261.777 2287.669  990.2398 2559.207
## 134       1792.241 1247.543 2336.940  959.1963 2625.287
## 135       1809.760 1232.773 2386.746  927.3354 2692.184
## 136       1827.278 1217.473 2437.082  894.6623 2759.893
## 137       1844.796 1201.646 2487.946  861.1831 2828.409
## 138       1862.314 1185.296 2539.332  826.9049 2897.724

Plot DES

plot(ramal_des, xlab="periode waktu", ylab="Yt",
     main="Forecasting Using Double Exponential Smoothing (0.8, 0.1)")
lines(testing.ts, col = "black", lwd=1.8)
lines (datades[,2], col="green",lwd=1)
legend("topleft", legend = c("Data aktual", "Fitted DES", "Ramalan"), lty=1, col= c("black", "green", "blue"))

Dari plot di atas, dapat terlihat garis data ramalan DES relatif mendekati garis data aktualnya (data testing). Selain itu, hampir seluruh data aktual masuk ke dalam selang kepercayaan 95% data ramalan DES. Hal ini mengindikasikan bahwa tingkat akurasi dari pemulusan DES sudah cukup baik dalam meramalkan data.

Nilai Akurasi DES

RMSE_DES_train <- sqrt(df_des$SSE/length(training.ts))
RMSE_DES_test <- sqrt(mean((testing.ts-ramal_des$mean)^2))
MAPE_DES_train <- mean(abs((df_des$fitted[,1]-training.ts)/training.ts)*100)
MAPE_DES_test <- mean(abs((testing.ts-ramal_des$mean)/testing.ts)*100)
Akurasi_DES <- data.frame("DES_Train"=c(RMSE_DES_train,MAPE_DES_train), "DES_Test"=c(RMSE_DES_test,MAPE_DES_test))
rownames(Akurasi_DES) <- c("RMSE","MAPE")
Akurasi_DES

##      DES_Train   DES_Test
## RMSE 67.873907 122.299878
## MAPE  3.619096   5.202502

Dari syntax di atas, didapatkan bahwa nilai RMSE DES testing sebesar 113.463494 dan MAPE DES testing sebesar 5.261905 . Kedua nilai error ini juga cukup besar jika dibandingkan data trainingnya.

Perbandingan DMA vs DES

plot(ramal_des, xlab="periode waktu", ylab="Yt",
     main="Forecasting DMA vs DES")
lines(testing.ts, col = "black", lwd=1.8)
lines (datades[,2], col="green",lwd=1.5)
lines(data.gab[,5],col="red",lwd=1.5)
legend("topleft", legend = c("Aktual", "Fitted DES", "Ramalan DES", "Ramalan DMA"), lty=1, col= c("black", "green", "blue", "red"))

Plot di atas menunjukkan hasil peramalan dengan dua metode yang berbeda, yaitu DMA dan DES. Peramalan pada metode DMA yang ditunjukkan oleh garis merah memiliki pola ramalan yang naik, begitu juga peramalan pada metode DES yang ditunjukkan oleh garis biru memiliki pola naik dan terdapat daerah yang merupakan selang kepercayaan 80% dan 95%. Dari kedua metode tersebut, dapat dilihat bahwa secara eksploratif hasil ramalan metode DES cenderung lebih mendekati pola data test dibandingkan hasil ramalan data DMA yang sangat naik dan menjauhi pola data test.

Nilai akurasi DMA vs DES

Akurasi <- data.frame("DMA"=Akurasi_DMA[,2], "DES"=Akurasi_DES[,2])
rownames(Akurasi) <- c("RMSE","MAPE")
Akurasi

##            DMA        DES
## RMSE 252.05395 122.299878
## MAPE  10.73407   5.202502

Dari perbandingan nilai keakuratan kedua model di atas, terbukti bahwa model DES menunjukkan nilai error yang lebih kecil daripada model DMA. Hal ini sesuai dengan hasil eksplorasi yang juga menyatakan bahwa model DES lebih baik untuk memprediksi data bila dibandingkan model DMA.

C. Cek Kestasioneran data

Secara Eksploratif

Secara eksploratif, kestasioneran data dapat terlihat dari plot deret waktu atau melihat plot acf dari peubah yang akan diuji. Jika pada plot deret waktu menunjukan pola tren naik/turun pada kurun waktu tertentu (tidak stabil), atau, jika plot acf menunjukkan pola tails off slowly, maka data tersebut tidak stasioner.

plot.ts(training.ts, ylab=expression(Y[t]), main = "Plot Time Series Data Penjualan Emas Bulanan")

acf(training.ts)

Terlihat dari plot deret waktu dan plot acf peubah data training di atas, bahwa secara eksploratif dapat disimpulkan data tersebut tidak stasioner karena menunjukkan pola tails off slowly. Uji ADF akan dilakukan untuk meyakinkan kesimpulan tersebut.

Uji Formal

Secara formal, metode Augmented Dickey-Fuller (ADF) dapat memberikan hasil uji secara akurat untuk menentukan apakah sebuah data stasioner atau tidak. Dengan hipotesis sebagai berikut:

H0 : Data tidak stasioner

H1 : Data stasioner

adf.test(training.ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  training.ts
## Dickey-Fuller = -1.9116, Lag order = 4, p-value = 0.6133
## alternative hypothesis: stationary

Didapat nilai P−value > 0.05 maka tak tolak H0 pada taraf nyata 5%. Artinya, cukup bukti untuk mengatakan bahwa data tidak stasioner. Untuk mengatasi ketidakstasioneran ini, perlu dilakukan differencing atau pembedaan.

Differencing

y_diff <- diff(training.ts, differences = 1)
y_diff

## Time Series:
## Start = 2 
## End = 117 
## Frequency = 1 
##   [1]   29.75    7.25   63.75   28.25   36.50  -75.00   77.00   61.00   39.75
##  [10]   36.75   22.00  -78.50   84.00   28.00   96.50    1.00  -31.00  123.00
##  [19]  185.00 -193.50  102.00   24.00 -215.00  213.00   26.00 -107.50  -11.25
##  [28]  -93.25   40.50   23.50   26.50  127.50  -57.00    7.00  -68.50    7.25
##  [37]  -76.25    9.75 -129.25  -74.50 -202.50  122.50   80.25  -68.25   -2.50
##  [46]  -71.00  -48.50   46.50   75.50  -34.75   -3.25  -38.00   64.50  -29.75
##  [55]    0.50  -69.25  -52.25   18.50   23.25   54.25  -46.25  -27.00   -6.75
##  [64]   11.15  -20.40  -72.60   36.60  -21.00   28.35  -80.45   -1.90   51.80
##  [73]  123.10    2.10   48.65  -73.55  108.65   21.25  -32.75   13.25  -50.50
##  [82]  -93.90  -32.20   66.90   42.80  -10.75   21.60   -0.25  -23.95   25.30
##  [91]   44.20  -28.65  -12.95   10.05   10.80   54.05  -27.20    6.00  -10.65
## [100]   -7.85  -54.90  -29.50  -18.50  -15.20   27.70    2.60   61.45   44.25
## [109]   -4.10  -23.75  -13.10   13.25  113.45   18.55  100.85  -43.10

Cek ulang kestasioneran data

adf.test(y_diff)

## Warning in adf.test(y_diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  y_diff
## Dickey-Fuller = -4.8468, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Didapat nilai P−value = 0.01 < 0.05 maka tolak H0 pada taraf nyata 5% atau cukup bukti untuk mengatakan bahwa data sudah stasioner. Artinya, differencing pertama sukses untuk membuat data stasioner.

D. Identifikasi model Arima

Untuk mengidentifikasi model arima, nilai ACF dan PACF dapat digunakan untuk menentukan nilai q pada nilai MA(q) dan nilai p pada model AR(p). Namun, kedua nilai ini p dan q pada model campuran ARMA(p,q). Oleh karena itu, dikembangkan metode extended autocorrelation function (EACF) untuk mengindentifikasi model campuran ARMA(p,q). Pada tabel EACF, triangle of zeros akan terbentuk, dan nilai pada pojok kiri atas akan bersesuaian dengan ordo ARMA.

ACF

acf(y_diff)

Dari plot ACF di atas, dapat terlihat bahwa terdapat cuts off setelah lag 2. Oleh karena itu, model yang terbentuk adalah ARIMA(0,2,2).

PACF

pacf(y_diff)

Dari plot PACF di atas, dapat terlihat bahwa terdapat cuts off setelah lag 3. Oleh karena itu, model yang terbentuk adalah ARIMA(3,2,0).

EACF

eacf(y_diff)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o x o o o o x  o  o  o 
## 1 x o o o o x o o o o x  o  o  o 
## 2 x o o o o x o o o o x  o  o  o 
## 3 x x o o o o o o o o x  o  o  o 
## 4 x x o x o o o o o o o  o  o  o 
## 5 o x x x o o o o o o x  o  o  o 
## 6 o o x o x o o o o o x  o  o  o 
## 7 o o x o o o o o o o o  o  o  o

Dari tabel EACF di atas, dapat terlihat bahwa terdapat triangle of zero dengan titik kiri atasnya yaitu pada 0,0; 0,1; 0,2; 1,1. Maka model yang terbentuk menurut EACF adalah ARIMA(0,1,0), ARIMA(0,1,1), ARIMA(0,1,2), ARIMA(1,1,1)

Pemilihan Model Terbaik

arima010 <- Arima(training.ts, order=c(0,1,0), method="ML")
arima011 <- Arima(training.ts, order=c(0,1,1), method="ML")
arima012 <- Arima(training.ts, order=c(0,1,2), method="ML")
arima111 <- Arima(training.ts, order=c(1,1,1), method="ML")

Keakuratan Model

akurasi_model = data.frame("Model" = c("ARIMA(0,1,0)", "ARIMA(0,1,1)", "ARIMA(0,1,2)", "ARIMA(1,1,1)"),
                        "AIC" = c(arima010$aic, arima011$aic, arima012$aic,arima111$aic),
                       "BIC" = c(arima010$bic, arima011$bic, arima012$bic,arima111$bic))
akurasi_model[order(akurasi_model[,3]),]

##          Model      AIC      BIC
## 1 ARIMA(0,1,0) 1310.662 1313.415
## 2 ARIMA(0,1,1) 1309.833 1315.341
## 3 ARIMA(0,1,2) 1311.714 1319.975
## 4 ARIMA(1,1,1) 1311.767 1320.028

Perbandingan Ramalan Dua Model Terbaik

ramalan010 <- forecast(arima010, h = 21)
ramalan011 <- forecast(arima011, h = 21)
plot(ramalan010); lines(testing.ts)

plot(ramalan011); lines(testing.ts)

Berdasarkan perbandingan kedua plot di atas, dapat terlihat bahwa model ARIMA(0,1,0) jauh lebih baik dalam menduga harga emas 21 pekan setelah data training. Oleh karena itu, secara eksploratif dapat dikatakan bahwa ARIMA(0,1,0) merupakan model terbaik.

Nilai Akurasi

akurasi010 <- accuracy(ramalan010, testing.ts)
akurasi011 <- accuracy(ramalan011, testing.ts)

akurasi_model = data.frame("Model" = c("ARIMA(0,1,0)", "ARIMA(0,1,1)"),
                        "RMSE" = c(akurasi010[2,2], akurasi011[2,2]),
                       "MAPE" = c(akurasi010[2,5], akurasi011[2,5]))
akurasi_model[order(akurasi_model[,3]),]

##          Model     RMSE     MAPE
## 2 ARIMA(0,1,1) 290.7865 14.00509
## 1 ARIMA(0,1,0) 294.4020 14.21936

Dari ukuran RMSE dan MAPE, terlihat bahwa model ramalan yang memiliki nilai RMSE dan MAPE terkecil adalah ARIMA(0,1,0), yaitu sebesar 245.4816 dan 11.31409.

Kesimpulan Model terbaik

Berdasarkan hasil yang telah didapat di atas, dapat disimpulkan bahwa model terbaik untuk melakukan ramalan data harga emas bulanan adalah ARIMA(0,1,0)

D. Overfitting

ARIMA(1,1,0)

arima110 <- Arima(training.ts, order=c(1,1,0), method="ML")
summary(arima110)

## Series: training.ts 
## ARIMA(1,1,0) 
## 
## Coefficients:
##           ar1
##       -0.1472
## s.e.   0.0917
## 
## sigma^2 = 4583:  log likelihood = -653.06
## AIC=1310.11   AICc=1310.22   BIC=1315.62
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE       ACF1
## Training set 4.049471 67.11879 50.38363 0.1903495 3.685552 1.001611 -0.0153007

coeftest(arima110)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.147219   0.091658 -1.6062   0.1082

Dari model overfitting ARIMA(1,1,0), dapat terlihat bahwa tidak ada parameter yang signifikan. Dapat terlihat jelas bahwa model ini tidak lebih baik daripada model awal.

ARIMA(0,1,1)

arima011 <- Arima(training.ts, order=c(0,1,1))
summary(arima011)

## Series: training.ts 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.1615
## s.e.   0.0935
## 
## sigma^2 = 4572:  log likelihood = -652.92
## AIC=1309.83   AICc=1309.94   BIC=1315.34
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 4.193709 67.03728 50.26426 0.1979633 3.678215 0.999238
##                      ACF1
## Training set -0.000460097

coeftest(arima011)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ma1 -0.161479   0.093493 -1.7272  0.08413 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari model overfitting ARIMA(0,1,1), dapat terlihat bahwa parameter ma1 yang ditambahkan tidak signifikan. Artinya, model ini juga tidak lebih baik daripada model awal.

Nilai AIC dan BIC Model

akurasi_ovrfit = data.frame("Model" = c("ARIMA(0,1,0)", "ARIMA(1,1,0)", "ARIMA(0,1,1)"),
                        "AIC" = c(arima010$aic, arima110$aic, arima011$aic),
                       "BIC" = c(arima010$bic, arima110$bic, arima011$bic))
akurasi_ovrfit[order(akurasi_ovrfit[,3]),]

##          Model      AIC      BIC
## 1 ARIMA(0,1,0) 1310.662 1313.415
## 3 ARIMA(0,1,1) 1309.833 1315.341
## 2 ARIMA(1,1,0) 1310.111 1315.618

Dari ukuran AIC, terlihat bahwa model yang memiliki nilai AIC terkecil adalah ARIMA(0,1,1), yaitu sebesar 1299.422 . Sedangkan, model yang memiliki nilai BIC terkecil adalah ARIMA(0,1,0).

Kesimpulan

Dari signifikansi model terbaik dan ukuran BIC, dapat disimpulkan bahwa model terbaiknya adalah model awal yaitu ARIMA(0,1,0).

E. Diagnostik Model

Prosedur Eksploratif

best_model <- arima010
sisaan <- best_model$residuals
par(mfrow=c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(as.numeric(best_model$fitted), as.numeric(sisaan), xlab="Fitted", ylab="Sisaan", 
     ylim=c(-500,500), main="Residuals vs Fitted"); abline(h=0, col="red")
acf(sisaan)
plot(1:length(sisaan),sisaan,type='o', xlab="Order", ylim=c(-500,500),
     main="Residuals vs Order"); abline(h=0, col="red")

Normal Q-Q Plot Berdasarkan hasil eksplorasi di atas, terlihat bahwa banyak amatan sisaan cenderung menjauhi garis qq-plot distribusi normal. Oleh karena itu, secara eksploratif dapat disimpulkan bahwa sisaan belum cukup menyebar normal.

Plot Residual vs Fitted Berdasarkan hasil eksplorasi di atas, mayoritas titik amatan berada disekitar nol meskipun ada beberapa titik amatan yang terletak jauh dari titik nol. Selain itu, dapat terlihat bahwa lebar pita dari sisaan tidak terlalu berbeda jauh dari kiri hingga ke kanan, hal ini mengindikasikan bahwa ragam dari sisaan cukup homogen.

Plot ACF Berdasarkan hasil eksplorasi di atas, pada plot ACF terdapat garis vertikal di lag tertentu yang melebihi tinggi garis biru horizontal.Oleh karena itu, secara eksploratif plot ACF dapat disimpulkan terdapat autokorelasi pada model.

Residual vs Order Berdasarkan hasil eksplorasi di atas, titik amatan pada plot kebebasan sisaan mayoritas bergerak di sekitar titik nol. Namun, terdapat beberapa titik amatan yang terletak cukup jauh dari titik nol. Oleh karena itu, secara eksploratif plot residual vs order belum dapat disimpulkan apakah sisaan tersebut terdapat autokorelasi atau tidak.

Uji Formal

Kolmogrov-Smirnov Test H0: Sisaan menyebar secara normal

H1: Sisaan tidak menyebar secara normal

ks.test(sisaan,"pnorm")

## Warning in ks.test(sisaan, "pnorm"): ties should not be present for the
## Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.50816, p-value < 2.2e-16
## alternative hypothesis: two-sided

Berdasarkan uji normalitas dengan Kolmogorov-Smirnov Test diperoleh nilai P−value=2.2e−16 < 0.05. Sehingga dapat disimpulkan bahwa data tidak menyebar secara normal.

L-Jung Box Test H0: Sisaan saling bebas (tidak terdapat korelasi)

H1: Sisaan tidak saling bebas (terdapat korelasi)

Box.test(sisaan,type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 2.737, df = 1, p-value = 0.09805

Berdasarkan uji autokorelasi dengan L-Jung Box Test diperoleh nilai P−value=0.09805 > 0.05. Sehingga, dapat disimpulkan bahwa tidak terdapat korelasi sisaan antar lagnya.

T-Test H0: Rataan sisaan sama dengan 0

H1: Rataan sisaan tidak sama dengan 0

t.test(sisaan, mu = 0, conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.55397, df = 116, p-value = 0.5807
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -8.977942 15.950224
## sample estimates:
## mean of x 
##  3.486141

Berdasarkan uji rataan dengan T-Test diperoleh nilai P−value=0.5807 > 0.05, sehingga dapat disimpulkan bahwa sisasaannya sama dengan 0.

Peramalan

Plot

ramalan_best <- forecast(best_model, h = 21)
plot(ramalan_best); lines(testing.ts)

Nilai Akurasi

accuracy(ramalan_best, testing.ts)

##                      ME     RMSE       MAE        MPE      MAPE      MASE
## Training set   3.486141  67.8673  49.88187  0.1598729  3.640099 0.9916362
## Test set     255.557143 294.4020 257.95238 14.0553234 14.219364 5.1280143
##                    ACF1 Theil's U
## Training set -0.1510093        NA
## Test set      0.7670305  3.324298