Analisis Deret Waktu pada Data Harga Minyak Mentah WTI Periode Januari 2012 - Agustus 2022
KELOMPOK 1
- Angelika Anggreni Batubara (G1401201040)
- Fidira Dwi Andari (G1401201059)
- Tahniah Salsabila Elenaputri (G1401201067)
- Nadia Aulia (G1401201079)
- Irma Wakhidatul Setyowati (G1401201089)
Package
Terdapat 9 packages yang digunakan, yaitu:
readxl, TTR, TSA,
tseries, graphics, ggplot2,
MASS, cowplot dan forecast
Latar Belakang
Minyak mentah awalnya merupakan campuran antara minyak bumi dengan gas alam, lalu dipisahkan dengan gas alam tersebut setelah proses penambangan (Wardhani dan Titah 2020). Minyak bumi menjadi salah satu faktor penggerak perekonomian dunia dan merupakan sumber daya alam yang tidak dapat diperbaharui. Permintaan masyarakat akan minyak bumi terus meningkat sehingga membuat ketersediaannya semakin menipis. Jika tidak diimbangi dengan ketersediaan minyak yang ada, maka minyak bumi akan menjadi sumber daya alam yang langka. Dampak yang terjadi antara ketidakseimbangan antara permintaan masyarakat dan ketersediaan minyak adalah harga minyak yang berfluktuatif.
AlMadi dan Zhang (2011), Pasar minyak mentah dunia memiliki 4 tolak ukur (benchmark) internasional utama yang digunakan untuk penetapan harga minyak mentah lainnya, yaitu 2 di West - West Texas Intermediate (WTI) di AS dan Brent di Eropa, dan 2 di Asia - Dubai dan Oman. Benchmark digunakan untuk menentukan harga ekspor minyak mentah setelah melakukan penyesuaian terkait perbedaan kualitas dan lokasi.
Untuk saat ini, harga minyak mentah umumnya mengikuti acuan West Texas Intermediate (WTI) atau light sweet dan Brent Crude Oil. WTI Crude Oil lebih ringan dan manis jika dibandingkan dengan Brent Crude Oil. Walaupun begitu, kedua jenis minyak tersebut termasuk ke dalam jenis minyak ringan dan manis. Sebelum terjadi krisis di Libya, minyak mentah WTI dijual lebih mahal dibandingkan dengan Brent Crude Oil karena minyak mentah WTI lebih ringan dan manis. Namun, setelah terjadinya krisis tersebut harga jual minyak mental WTI menjadi lebih murah dibandingkan Brent Crude Oil. Produksi minyak mentah WTI juga mengalami peningkatan sehingga produksi dan peran Brent Crude Oil semakin menurun. Selisih yang terjadi antara harga minyak mentah WTI dengan Brent Crude Oil membuat pembeli minyak internasional lebih tertarik untuk membeli minyak mentah WTI.
Analisis Data
Deskripsi Data dan Time Series Plot
Data yang digunakan dalam pemodelan deret waktu ini adalah sebanyak 131 amatan yaitu Harga Minyak WTI untuk Januari 2012 hingga November 2022 dengan periode pengamatan bulanan. Dengan Time Series Plot sebagai berikut:
data <- read_excel("C:\\Users\\ASUS\\Downloads\\PROJECT MPDW\\NEW WTI Crude Oil Data1.xlsx", sheet = "Sheet2")
dataa <- read_excel("C:\\Users\\ASUS\\Downloads\\PROJECT MPDW\\NEW WTI Crude Oil Data1.xlsx", sheet = "Sheet6")
ts.dt <- ts(data)
ggplot(dataa, aes(x = Tanggal, y = Terakhir)) +
geom_line(size = 1) +
ggtitle("Time Series Plot of WTI Crude Oil Prices")+
labs(x = "Time Period", y = "Crude Oil Prices")+
theme_grey()+
theme(plot.title = element_text(hjust = 0.5))
Pergerakan harga minyak mentah cenderung nonstasioner serta dipengaruhi oleh banyak faktor. Berbagai metode peramalan dapat memberikan hasil yang baik jika deretan harga yang terbentuk merepresentasikan fungsi yang cenderung linier. Namun, fluktuasi harga riil minyak mentah cenderung memiliki karakteristik nonlinier dan nonstasioner. Hal ini disebabkan oleh kompleksitas intrinsik dari pasar minyak (Herawati dan Djunaidy 2014).
Eksplorasi Data
par(mfrow = c(2,2))
gp1 <- ggplot(dataa[1:24,], aes(x = Tanggal, y = Terakhir)) +
geom_line(size = 1) +
ggtitle("Januari 2012 - June 2014")+
labs(x = "Time Period", y = "Crude Oil Prices")+
theme_grey()+
theme(plot.title = element_text(hjust = 0.5))
gp2 <- ggplot(dataa[31:96,], aes(x = Tanggal, y = Terakhir)) +
geom_line(size = 1) +
ggtitle("Juli 2014 - Desember 2019")+
labs(x = "Time Period", y = "Crude Oil Prices")+
theme_grey()+
theme(plot.title = element_text(hjust = 0.5))
gp3 <- ggplot(dataa[97:108,], aes(x = Tanggal, y = Terakhir)) +
geom_line(size = 1) +
ggtitle("Januari 2020 - Desember 2020")+
labs(x = "Time Period", y = "Crude Oil Prices")+
theme_grey()+
theme(plot.title = element_text(hjust = 0.5))
gp4 <- ggplot(dataa[109:131,], aes(x = Tanggal, y = Terakhir)) +
geom_line(size = 1) +
ggtitle("Januari 2021 - November 2022")+
labs(x = "Time Period", y = "Crude Oil Prices")+
theme_grey()+
theme(plot.title = element_text(hjust = 0.5))
plot_grid(gp1, gp2, gp3, gp4)
Berdasarkan plot time series di atas, terlihat bahwa data merupakan data gabungan stasioner dan trend naik di akhir periode.
Terlihat pula dari ketiga plot diatas bahwa harga minyak WTI berfluktuatif. Seperti yang ditunjukkan pada plot untuk periode waktu Juli 2014-Desember 2019, dimana pada awal periode terjadi penurunan yang begitu signifikan. Lonjakan pasokan minyak dan minimnya permintaan merupakan faktor utama di belakang kejatuhan harga minyak semenjak pertengahan 2014. Ditambah lagi bangkitnya minyak shale dari Amerika Serikat dan juga penolakan OPEC untuk menurunkan produksinya. Selain itu pada awal tahun 2016 juga terjadi penurunan, dimana pada Februari 2016 harga minyak WTI hanya sebesar 33.75 USD. Yang mana jika dibandingkan sejak Juli 2014, di mana harga minyak berada di puncaknya, nilai ini sudah jatuh 70 persen (Bareksa.com)
Dan pada plot untuk periode waktu Januari 2020-November 2022, dimana masa pandemic yang ikut memengaruhi harga minyak WTI. Menurut Widyastuti N L dan Nugroho Hanan (Juni 2020), yaitu Permintaan minyak mentah yang sudah mulai terjadi pada kuartal pertama tahun 2020, akan terus menurun hingga mencapai titik terbawahnya pada akhir kuartal ke-2, kemudian diperkirakan naik lagi dan mencapai titik keseimbangannya dengan pasokan (supply) minyak bumi pada akhir kuartal ke-3 tahun 2020. Dalam masa penurunan permintaan tersebut, minyak yang dikonsumsi akan berkurang hingga lebih 10 juta barel per hari, atau bahkan 15-20 juta barel per hari. Atau sekitar 20%-30% menurut beberapa analisis yang lain. Produksi minyak mentah tidak mungkin bisa diturunkan seketika, dalam masa penurunan permintaan yang tajam tersebut akan terbentuk stock minyak yang cukup besar. Hal ini membuat harga minyak turun.
Namun turunnya harga minyak bukan saja karena dampak Covid-19, tapi juga didorong oleh “pertikaian” (conflict) dalam industri minyak sendiri. Pertikaian, khususnya dalam kelompok produsen, diawali oleh Saudi Arabia dan Rusia, menyangkut kesepakatan mengenai berapa jumlah minyak mentah yang seharusnya diproduksi.
Splitting Data
Dilakukan splitting data atau membagi data menjadi data training dan data testing. Data training digunakan dalam proses konstruksi model atau digunakan untuk melatih dan membangun model, sedangkan data testing digunakan untuk menguji performa model yang didapatkan pada saat tahapan training.
Fluktuasi yang terjadi pada minyak mentah WTI terutama pada tahun 2012 hingga 2020, sedangkan setelah jatuhnya harga pada pada 2020, harga minyak mentah WTI terus mengalami kenaikan. Untuk itu harga minyak mentah WTI untuk periode waktu Januari 2012 hingga Juni 2021 akan digunakan untuk memodelkan, di mana pada data pemodelan ini sudah meliputi pola kenaikan data yang terjadi (diwakili oleh data pada bulan Mei 2020 hingga Juni 2021). Dan data pada periode waktu Juli 2021 hingga November 2022 akan digunakan untuk menguji model peramalan nantinya.
train <- 1:114
y <- data$Terakhir
train.ts <- ts(y[train])
test.ts <- ts(y[-train], start = 115)plot.train <- ggplot(dataa[1:114,], aes(x = Tanggal, y = Terakhir)) +
geom_line(size = 1) +
ggtitle("Time Series Plot of Data Train")+
labs(x = "Time Period", y = "Crude Oil Prices")+
theme_grey()+
theme(plot.title = element_text(hjust = 0.5))
plot.test <- ggplot(dataa[115:131,], aes(x = Tanggal, y = Terakhir)) +
geom_line(size = 1) +
ggtitle("Time Series Plot of Data Test")+
labs(x = "Time Period", y = "Crude Oil Prices")+
theme_grey()+
theme(plot.title = element_text(hjust = 0.5))
plot_grid(plot.train, plot.test)
Pemulusan
Double Moving Average
Metode peramalan moving average dilakukan dengan mengambil rata-rata sekelompok nilai pengamatan, lalu menggunakan rata-rata tersebut sebagai ramalan untuk periode berikutnya. Istilah rata-rata bergerak digunakan karena setiap kali data observasi baru tersedia, angka rata-rata yang baru dihitung dan dipergunakan sebagi ramalan. Salah satu cara untuk meramalkan data time series yang tidak stasioner atau memiliki kecenderungan trend adalah menggunakan double moving average atau rata-rata bergerak ganda. Dasar metode ini adalah menghitung rata-rata bergerak (moving average) sebanyak dua kali.
Pola data yang cenderung menunjukkan pola trend, sehingga untuk pemulusan digunakan Double Moving Average(DMA).
Penentuan n optimum
Dalam pemulusan dengan menggunakan Double Moving Average perlu ditentukan nilai n optimum agar model yang didapatkan menjadi model yang baik. Parameter yang optimum ditunjukkan dengan ukuran keakuratan SSE, RMSE, atau MAPE yang minimum.
colnames(akurasi) = c("n", "SSE", "RMSE", "MAPE")
akurasi[order(akurasi[,4]),]## n SSE RMSE MAPE
## 1 2 6655.627 7.743423 10.89529
## 2 3 8473.291 8.816837 12.33606
## 3 4 9746.306 9.543950 14.32823
## 4 5 10050.303 9.783515 14.55041
## 5 6 10499.016 10.096148 14.69279
## 6 7 11033.105 10.451730 15.77016
## 7 8 11927.440 10.976301 16.92944
## 8 9 12969.381 11.563086 18.20430
## 9 10 13795.998 12.050769 19.10420Fungsi tersebut menampilkan nilai SSE, RMSE, dan MAPE untuk setiap nilai n dalam rentang 2 sampai 10. Didapatkan kesimpulan bahwa nilai SSE, RMSE, dan MAPE minimum terjadi pada saat n = 2.
Pemulusan dengan n Optimum
n = 2
data.sma <- SMA(train.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:17) {
fcast[i] = At[length(At)] + Bt[length(Bt)]*(i)
} # Ramalan DMA
join.data <- cbind(aktual = c(train.ts, rep(NA, 17)),
pemulusan = c(dma, rep(NA, 17)),
At = c(At, rep(NA, 17)), Bt = c(Bt, rep(NA, 17)),
ramalan = c(data.ramal, fcast[-1]))Plot DMA dengan n optimum
plot.ts(data, xlab = "Periode Waktu", ylab="Price",
main = "Double Moving Average n = 2")
lines(join.data[,2], col = "maroon", lwd = 2)
lines(join.data[,5], col = "green", lwd = 1.5)
legend("bottomleft",
c("Aktual", "Pemulusan", "Peramalan"), lty=1,
col = c("black", "maroon", "green"), cex=0.8)
Plot di atas menunjukkan hasil pemulusan dan peramalan metode DMA pada n = 2. Dapat dilihat bahwa garis pemulusan yang berwarna merah sangat menggambarkan pola data aktual. Sementara itu, hasil peramalan menunjukkan pola yang semakin menaik.
Perbandingan Akurasi Data Training dan Testing DMA
baris.dma <- c("Training", "Testing")
rmse.dma.train <- sqrt(mean((train.ts-data.ramal[-115])^2, na.rm = T))
rmse.dma.test <- sqrt(mean((test.ts-fcast)^2))
mape.dma.train <- mean(abs((train.ts-data.ramal[-115])/train.ts)*100, na.rm = T)
mape.dma.test <- mean(abs((test.ts-fcast)/test.ts)*100)
rmse_dma <- c(rmse.dma.train, rmse.dma.test)
mape_dma <- c(mape.dma.train, mape.dma.test)
akurasi.dma <- data.frame(baris.dma, rmse_dma, mape_dma)
kableExtra::kable(akurasi.dma, align = c(rep('c', times = 3)),
col.names = c("DMA", "RMSE", "MAPE"), digits = 3)| DMA | RMSE | MAPE |
|---|---|---|
| Training | 7.743 | 10.895 |
| Testing | 35.316 | 33.355 |
Berdasarkan tabel di atas, didapatkan bahwa nilai RMSE dan MAPE data testing dengan metode DMA sangat besar dibandingkan dengan nilai RMSE dan MAPE data training.
Double Exponential Smoothing
Metode peramalan exponential smooting merupakan metode yang digunakan untuk meramalkan data yang mengalami trend kenaikan dan apabila data yang digunakan semakin banyak dalam perhitungan peramalannya maka percentase error peramalannya akan semakin kecil, begitu juga sebaliknya (Purwanto dan Hanief 2017).
Penentuan α dan β optimum
Dalam pemulusan dengan menggunakan Double Exponential Smoothing perlu ditentukan nilai parameter (α dan β) optimum agar model yang didapatkan menjadi model yang baik. Parameter yang optimum ditunjukkan dengan ukuran keakuratan SSE, RMSE, atau MAPE yang minimum.
hasil <- cbind("Alpha" = rep(a, each=9), "Beta" = b, hasil)
dplyr::arrange(.data = hasil, MAPE)## Alpha Beta SSE RMSE MAPE
## 1 0.9 0.1 5194.014 6.749928 9.241108
## 2 0.9 0.2 5106.813 6.693026 9.305813
## 3 0.9 0.3 5197.202 6.751999 9.420838
## 4 0.9 0.4 5320.670 6.831731 9.487223
## 5 0.8 0.2 5486.979 6.937680 9.588758
## 6 0.9 0.5 5451.059 6.914933 9.705719
## 7 0.8 0.3 5558.518 6.982759 9.713485
## 8 0.8 0.1 5641.057 7.034412 9.769917
## 9 0.8 0.4 5670.898 7.052994 9.790525
## 10 0.9 0.6 5583.550 6.998465 9.873080
## 11 0.8 0.5 5787.796 7.125317 9.897237
## 12 0.8 0.6 5897.872 7.192755 10.083634
## 13 0.9 0.7 5721.062 7.084120 10.086218
## 14 0.7 0.3 6079.221 7.302499 10.111655
## 15 0.7 0.2 6023.454 7.268928 10.155877
## 16 0.7 0.4 6193.768 7.370977 10.165780
## 17 0.8 0.7 5997.579 7.253299 10.266035
## 18 0.7 0.5 6319.698 7.445532 10.273764
## 19 0.7 0.1 6267.741 7.414862 10.349416
## 20 0.9 0.8 5871.227 7.176489 10.368270
## 21 0.7 0.6 6440.661 7.516451 10.379638
## 22 0.8 0.8 6088.059 7.307806 10.425375
## 23 0.8 0.9 6174.164 7.359303 10.577281
## 24 0.7 0.7 6547.818 7.578720 10.593549
## 25 0.9 0.9 6045.088 7.281970 10.642068
## 26 0.6 0.2 6763.494 7.702525 10.807393
## 27 0.7 0.8 6636.285 7.629746 10.815344
## 28 0.6 0.3 6792.979 7.719296 10.853247
## 29 0.6 0.5 7051.588 7.864861 10.883741
## 30 0.6 0.4 6909.696 7.785330 10.922725
## 31 0.6 0.6 7199.850 7.947112 10.948267
## 32 0.6 0.7 7342.876 8.025659 10.961874
## 33 0.7 0.9 6704.650 7.668945 10.985659
## 34 0.6 0.1 7149.641 7.919353 11.029098
## 35 0.6 0.8 7471.338 8.095557 11.167359
## 36 0.6 0.9 7578.360 8.153333 11.434281
## 37 0.5 0.2 7812.536 8.278346 11.544366
## 38 0.5 0.3 7787.114 8.264867 11.700159
## 39 0.5 0.4 7892.308 8.320503 11.799760
## 40 0.5 0.1 8445.102 8.606965 11.847232
## 41 0.5 0.5 8037.187 8.396526 11.924728
## 42 0.5 0.6 8204.003 8.483215 11.991878
## 43 0.5 0.7 8383.798 8.575668 12.020374
## 44 0.5 0.8 8566.233 8.668471 12.050821
## 45 0.5 0.9 8741.122 8.756512 12.176426
## 46 0.4 0.2 9409.252 9.085004 12.314966
## 47 0.4 0.3 9265.525 9.015350 12.482174
## 48 0.4 0.4 9343.546 9.053228 12.848157
## 49 0.4 0.5 9466.724 9.112708 12.987547
## 50 0.4 0.1 10513.841 9.603472 13.032299
## 51 0.4 0.6 9610.481 9.181638 13.089319
## 52 0.4 0.7 9784.689 9.264481 13.171349
## 53 0.4 0.8 9987.581 9.360040 13.408091
## 54 0.4 0.9 10208.994 9.463222 13.565854
## 55 0.3 0.3 11696.993 10.129423 13.710161
## 56 0.3 0.2 12174.040 10.333916 13.888785
## 57 0.3 0.4 11723.418 10.140858 13.922606
## 58 0.3 0.5 11889.062 10.212248 14.440369
## 59 0.3 0.6 12008.662 10.263486 14.777302
## 60 0.3 0.7 12084.219 10.295723 14.935447
## 61 0.3 0.1 14312.753 11.204930 14.956582
## 62 0.3 0.8 12181.772 10.337198 15.025787
## 63 0.3 0.9 12339.840 10.404048 15.268242
## 64 0.2 0.4 16072.618 11.873832 16.501045
## 65 0.2 0.3 16815.211 12.145034 16.571704
## 66 0.2 0.5 16064.138 11.870699 16.690348
## 67 0.2 0.6 16559.821 12.052451 17.043677
## 68 0.2 0.2 18303.217 12.671012 17.155456
## 69 0.2 0.7 17137.460 12.260856 17.720168
## 70 0.2 0.9 17458.303 12.375096 17.925703
## 71 0.2 0.8 17457.269 12.374730 17.972183
## 72 0.2 0.1 23240.446 14.278080 18.460122
## 73 0.1 0.8 25800.110 15.043828 21.226777
## 74 0.1 0.9 25071.526 14.829891 21.321475
## 75 0.1 0.7 27656.572 15.575670 21.576124
## 76 0.1 0.3 34246.130 17.332185 23.206421
## 77 0.1 0.6 30245.603 16.288410 23.273420
## 78 0.1 0.2 38463.161 18.368346 23.421831
## 79 0.1 0.5 32944.224 16.999542 24.364452
## 80 0.1 0.4 34194.339 17.319075 24.422220
## 81 0.1 0.1 63568.917 23.614023 32.350657Fungsi tersebut menampilkan nilai SSE, RMSE, dan MAPE untuk setiap nilai alpha dan beta dalam rentang 0.1 sampai 0.9. Dengan mengacu pada MAPE, didapatkan parameter optimum adalah α = 0.9 dan β = 0.1.
Pemulusan dengan α dan β Optimum
df_des <- HoltWinters(train.ts, alpha = .9, beta = .1, gamma = F)
datades <- data.frame(train.ts, c(NA, NA, df_des$fitted[,1]))
colnames(datades) = c("train.ts", "train.ts_smt")
datades## train.ts train.ts_smt
## 1 98.48 NA
## 2 107.07 NA
## 3 103.02 115.66000
## 4 104.87 111.73640
## 5 86.53 112.39106
## 6 84.96 93.62303
## 7 88.06 89.55356
## 8 96.47 91.80219
## 9 92.19 100.01616
## 10 86.24 96.28120
## 11 88.91 89.64900
## 12 91.82 91.32227
## 13 97.49 94.15339
## 14 92.05 99.83980
## 15 97.23 94.81136
## 16 93.46 99.18819
## 17 91.97 95.71734
## 18 96.56 93.69199
## 19 105.03 97.87858
## 20 107.65 106.56386
## 21 102.33 109.88814
## 22 96.38 104.75234
## 23 92.72 98.13025
## 24 98.42 93.68712
## 25 97.49 98.79876
## 26 102.59 98.35514
## 27 101.58 103.28191
## 28 99.74 102.71242
## 29 102.71 100.73195
## 30 105.37 103.38493
## 31 98.17 106.22288
## 32 95.96 99.30192
## 33 91.16 96.32005
## 34 80.54 91.23746
## 35 66.15 80.20843
## 36 53.27 64.88927
## 37 48.24 50.71962
## 38 49.76 44.55249
## 39 47.60 45.77245
## 40 59.63 44.11493
## 41 60.30 56.17253
## 42 59.47 58.35276
## 43 47.12 57.92434
## 44 49.20 45.79410
## 45 45.09 46.75961
## 46 46.59 43.00690
## 47 41.65 44.30411
## 48 37.04 39.74896
## 49 33.62 34.90064
## 50 33.75 31.22255
## 51 38.34 31.19921
## 52 45.92 35.97055
## 53 49.10 44.16513
## 54 48.33 48.29073
## 55 41.60 48.01382
## 56 44.70 41.35189
## 57 48.24 43.77702
## 58 46.86 47.60721
## 59 49.44 46.68097
## 60 53.72 49.15866
## 61 52.81 53.66895
## 62 54.01 53.22368
## 63 50.60 54.32992
## 64 49.33 51.03585
## 65 48.32 49.40992
## 66 46.04 48.24023
## 67 50.17 45.87324
## 68 47.23 49.74025
## 69 51.67 47.25503
## 70 54.38 51.39985
## 71 57.40 54.52155
## 72 60.42 57.81078
## 73 64.73 61.09253
## 74 61.64 65.62708
## 75 64.94 62.94070
## 76 68.57 65.82200
## 77 67.04 69.62445
## 78 74.15 68.39509
## 79 68.76 75.18910
## 80 69.80 70.43888
## 81 73.25 70.84236
## 82 65.31 74.20439
## 83 50.93 66.59410
## 84 45.41 51.48130
## 85 53.79 44.45561
## 86 57.22 52.13513
## 87 60.14 56.44772
## 88 63.91 59.83929
## 89 53.50 63.93781
## 90 58.47 54.03926
## 91 58.58 57.92117
## 92 55.10 58.46765
## 93 54.07 55.08721
## 94 54.18 53.73062
## 95 55.17 53.73441
## 96 61.06 54.75499
## 97 51.56 60.72550
## 98 44.76 51.94765
## 99 20.48 44.30298
## 100 18.84 19.54244
## 101 35.49 15.52717
## 102 39.27 31.90730
## 103 40.27 37.60995
## 104 42.61 39.31962
## 105 40.22 41.89272
## 106 35.79 39.84849
## 107 45.34 35.29180
## 108 48.52 44.33547
## 109 52.20 48.47845
## 110 61.50 52.53968
## 111 59.16 62.12224
## 112 63.58 60.70789
## 113 66.32 64.80294
## 114 73.47 67.81499ramal_des <- forecast::forecast(df_des, h=17)
data.frame(ramal_des)## Point.Forecast Lo.80 Hi.80 Lo.95 Hi.95
## 115 75.06014 66.33224 83.78804 61.71197 88.40831
## 116 77.21578 64.93423 89.49734 58.43276 95.99880
## 117 79.37142 63.88955 94.85330 55.69393 103.04892
## 118 81.52707 62.98074 100.07339 53.16291 109.89123
## 119 83.68271 62.12108 105.24433 50.70705 116.65837
## 120 85.83835 61.26733 110.40937 48.26022 123.41648
## 121 87.99399 60.39519 115.59280 45.78526 130.20272
## 122 90.14963 59.48997 120.80930 43.25973 137.03954
## 123 92.30528 58.54240 126.06815 40.66941 143.94114
## 124 94.46092 57.54644 131.37540 38.00509 150.91674
## 125 96.61656 56.49811 136.73501 35.26069 157.97243
## 126 98.77220 55.39480 142.14961 32.43219 165.11222
## 127 100.92784 54.23480 147.62089 29.51699 172.33870
## 128 103.08349 53.01705 153.14992 26.51348 179.65349
## 129 105.23913 51.74094 158.73732 23.42071 187.05755
## 130 107.39477 50.40619 164.38335 20.23826 194.55128
## 131 109.55041 49.01275 170.08808 16.96604 202.13478Plot DES dengan α dan β optimum
plot(ramal_des, xlab = "Periode Waktu", ylab = "Price",
main = "Double Exponential Smoothing α = 0.9 dan β = 0.1")
lines(ts.dt, col = "black")
lines (datades[,2], col = "maroon", lwd = 1)
legend("bottomleft", legend = c("Data aktual", "Pemulusan", "Peramalan"),
lty=1, col = c("black", "maroon", "blue"))
Berdasarkan time series plot DES, plot data training hasil fitted menggunakan metode pemulusan DES dengan α = 0.9 dan β = 0.1 menampilkan plot yang berpola sangat mendekati plot data aktual. Selain itu, hampir seluruh data aktual masuk ke dalam selang kepercayaan 95% dari data ramalan DES. Hal ini mengindikasikan bahwa tingkat peramalan menggunakan metode DES terbilang cukup baik dalam meramalkan data.
Perbandingan Akurasi Data Training dan Testing DES
baris.des <- c("Training", "Testing")
rmse.des.train <- sqrt(df_des$SSE/length(train.ts))
rmse.des.test <- sqrt(mean((test.ts-ramal_des$mean)^2))
mape.des.train <- mean(abs((df_des$fitted[,1]-train.ts)/train.ts)*100)
mape.des.test <- mean(abs((test.ts-ramal_des$mean)/test.ts)*100)
rmse_des <- c(rmse.des.train, rmse.des.test)
mape_des <- c(mape.des.train, mape.des.test)
akurasi.des <- data.frame(baris.des, rmse_des, mape_des)
kableExtra::kable(akurasi.des, align = c(rep('c', times = 3)),
col.names = c("DES", "RMSE", "MAPE"), digits = 3)| DES | RMSE | MAPE |
|---|---|---|
| Training | 6.750 | 9.241 |
| Testing | 12.625 | 11.937 |
Perbandingan DMA (n = 2) dan DES (α = 0.9 dan β = 0.1)
Plot
plot(ramal_des, xlab = "Periode Waktu", ylab = "Price",
main = "Plot Peramalan DMA vs DES")
lines(ts.dt, col = "black", lwd = 1.5)
lines(join.data[,5], col = "green", lwd = 1.5)
legend("bottomleft",
legend = c("Data Aktual", "Ramalan DMA", "Ramalan DES"),
lty = 1, col = c("black", "green", "blue"))
Plot di atas menunjukkan hasil peramalan dengan dua metode yang berbeda, yaitu DMA dan DES. Peramalan dengan metode DMA memiliki pola ramalan yang naik secara ekstrem dan menjauhi pola data aktual, seperti ditunjukkan oleh garis hijau. Sementara itu, pola peramalan dengan metode DES (garis berwarna biru) cenderung menurun dan tidak terlalu jauh dengan data aktual. Secara eksploratif, dapat disimpulkan hasil ramalan dengan menggunakan metode DES lebih baik karena cenderung lebih mendekati pola data testing (aktual) dibandingkan hasil ramalan dengan metode DMA.
Perbandingan Akurasi
baris <- c("DMA (n = 2)", "DES (α = 0.9, β = 0.1)")
rmse_1 <- c(rmse.dma.test, rmse.des.test)
mape_1 <- c(mape.dma.test, mape.des.test)
dmavsdes <- data.frame(baris, rmse_1, mape_1)
kableExtra::kable(dmavsdes, align = c(rep('c', times = 3)),
col.names = c(" ", "RMSE", "MAPE"), digits = 3)| RMSE | MAPE | |
|---|---|---|
| DMA (n = 2) | 35.316 | 33.355 |
| DES (α = 0.9, β = 0.1) | 12.625 | 11.937 |
Berdasarkan tabel nilai keakuratan di atas, dapat dilihat bahwa model DES (α = 0.9, β = 0.1) menunjukkan nilai RMSE dan MAPE yang jauh lebih kecil daripada model DMA (n = 2). Oleh karena itu, model dengan metode DES lebih tepat digunakan untuk memprediksi data yang ada dibandingkan dengan metode DMA.
Cek Kestasioneran
Hasil prediksi akan mendekati data aktual atau tidak melenceng jauh apabila pola data yang digunakan stasioner. Uji stasioneritas dilakukan untuk mengetahui apakah terdapat pola stasioner pada data yang digunakan. Data time series dikatakan stasioner jika rata-rata, varian, dan kovarian pada setiap lag adalah tetap sama pada setiap waktu. Kestasioneran dalam ragam dapat dilihat dengan menggunakan Box-Cox. Sementara itu, terdapat beberapa metode untuk uji stasioneritas dalam rataan, seperti uji ADF (Augmented Dickey-Fuller) dan uji PP (Philip Peron).
Kestasioneran dalam Ragam
Stasioneritas terhadap ragam dapat dilihat dengan plot Box-Cox. Data deret waktu dikatakan stasioner terhadap ragam apabila berfluktuasi dengan ragam konstan. Apabila stasioneritas tersebut tidak dipenuhi maka perlu dilakukan transformasi, yaitu transformasi pangkat (power transformation) atau disebut dengan transformasi Box-Cox. Apabila selang kepercayaan 95% bagi λ memuat satu, maka dikatakan bahwa data deret waktu telah stasioner terhadap ragam.
b <- boxcox(dataa$Terakhir[1:114]~dataa$Tanggal[1:114], lambda = seq(3,-5))
lambda <- b$x[which.max(b$y)]
lambda## [1] 0.6565657Berdasarkan output di atas, didapatkan nilai λ adalah 0.65657 dengan selang kepercayaan 95% dari λ memuat nilai satu. Oleh karena itu, disimpulkan bahwa data sudah stasioner dalam ragam.
Selain itu, untuk melihat apakah data sudah stasioner dalam ragam dapat menggunakan Fligner-Killeen test. Fligner-Killeen test adalah suatu uji untuk melihat kehomogenan ragam ketika data yang digunakan tidak berdistribusi normal atau terdapat pencilan pada data.
Hipotesis:
H0 : Ragam memenuhi asumsi homogenitas
H1 : Ragam tidak memenuhi asumsi homogenitas
data7 <- read_excel("C:\\Users\\ASUS\\Downloads\\PROJECT MPDW\\NEW WTI Crude Oil Data1.xlsx", sheet = 7)
data.ts <- ts(data7, start = c(2012, 1), frequency = 12)
library(car)
fligner.test(Yt ~ Tahun, data = data.ts)##
## Fligner-Killeen test of homogeneity of variances
##
## data: Yt by Tahun
## Fligner-Killeen:med chi-squared = 8.5188, df = 9, p-value = 0.4828Berdasarkan hasil uji Fligner-Killeen dengan menggunakan taraf signifikansi α = 5%, diperoleh nilai p−value = 0.4828. Karena p−value > α = 0.05, H0 tidak ditolak sehingga cukup bukti untuk menyatakan bahwa ragam sudah stasioner.
Pemeriksaan kestasioneran dalam rata-rata dilakukan setelah memeriksa kestasioneran data dalam ragam. Dilakukan prosedur secara eksploratif menggunakan grafik dan uji formal.
Prosedur Eksploratif (Plot ACF & PACF)
par(mfrow=c(1,2))
acf(train.ts, main = "ACF Plot for MA(q)")
pacf(train.ts, main = "PACF Plot for AR(p)")
Terlihat bahwa plot ACF membentuk tails off yang menurun secara bertahap sehingga disimpulkan bahwa data yang digunakan tidak stasioner dalam rata-rata.
Uji Formal (ADF Test)
Untuk memastikan hasil eksplorasi melalui plot ACF sebelumnya,
dilakukan pengujian formal yaitu Uji ADF untuk menguji kestasioneran
data dengan hipotesis sebagai berikut:
H0 : data tidak stasioner
H1 : data stasioner
adf.test(train.ts)##
## Augmented Dickey-Fuller Test
##
## data: train.ts
## Dickey-Fuller = -1.3277, Lag order = 4, p-value = 0.8557
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF, didapatkan p-value = 0.8557 > alpha = 0.05 sehingga H0 tidak ditolak. Artinya, belum cukup bukti untuk menyatakan bahwa data yang digunakan stasioner. Oleh karena itu, dilakukan differencing untuk membuat data menjadi stasioner.
Differencing
Suatu deret waktu yang tidak stasioner harus diubah menjadi data
stasioner dengan melakukan differencing, yaitu menghitung perubahan atau
selisih nilai observasi. Nilai selisih yang diperoleh diperiksa lagi
apakah stasioner atau tidak. Jika belum stasioner, dilakukan
differencing lagi.
Berdasarkan hasil pemeriksaan sebelumnya yang menunjukkan bahwa data
tidak stasioner, maka dilakukan differencing sebanyak satu kali pada
data tersebut.
dt.dif1 <- diff(train.ts, differences = 1)Plot Setelah Differencing
ts.plot(dt.dif1,
main = "Plot Hasil Differencing d = 1")
Bedasarkan plot di atas, terlihat bahwa data sudah stasioner setelah dilakukan differencing sebanyak satu kali
Uji ADF
Dilakukan kembali uji ADF pada data hasil differencing untuk
memeriksa kestasioneran data setelah differencing dengan hipotesis
sebagai berikut:
H0 : data tidak stasioner
H1 : data stasioner
adf.test(dt.dif1)##
## Augmented Dickey-Fuller Test
##
## data: dt.dif1
## Dickey-Fuller = -4.7609, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF, didapatkan p-value = 0.01 < alpha = 0.05 sehingga H0 ditolak. Artinya, data hasil differencing dengan diff = 1 menghasilkan pola data yang stasioner.
Plot ACF & PACF
par(mfrow=c(1,2))
acf(dt.dif1, main = "ACF Plot for MA(q)")
pacf(dt.dif1, main = "PACF Plot for AR(p)")
Plot ACF dan PACF di atas tidak menunjukkan adanya pola tails off sehingga sudah stasioner. Selain itu, plot di atas juga tidak menunjukkan adanya pola cuts off sehingga model tentatif dilihat dengan menggunakan plot EACF. Akan tetapi, apabila hal tersebut diabaikan, terjadi cut off after lag 3 pada plot ACF dan PACF sehingga didapatkan kandidat model ARIMA(0,1,3) dan ARIMA(3,1,0).
Plot EACF for ARMA(p,q)
eacf(dt.dif1)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o o o o
## 1 o o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 o x x o o o o o o o o o o o
## 4 x o x o o o o o o o o o o o
## 5 x x x o o o o o o o o o o o
## 6 x x o o o o o o o o o o o o
## 7 x x o o o o o o o o o o o oBerdasarkan plot EACF di atas, triangles of zero terdapat pada ARIMA (0,1,1)
Berdasarkan plot ACF, PACF, maupun EACF, didapatkan tiga model tentatif, yaitu:
- ARIMA(0,1,3)
- ARIMA(0,1,3)
- ARIMA(3,1,0)
Model ARIMA(p,d,q)
Metode Autoregressive Integrated Moving Average (ARIMA) sering juga disebut metode runtun waktu Box-Jenkins. ARIMA menggunakan nilai masa lalu dan sekarang dari variabel dependen untuk menghasilkan peramalan jangka pendek yang akurat. Model ARIMA sebagai fungsi dari p, d, q, di mana p sebagai orde operator dari AR, d orde differencing, dan q sebagai orde operator dari MA.
Pada tahap ini dilakukan pendugaan koefisien bagi model-model yang telah diidentifikasi sebelumnya
ARIMA(0,1,1)
model1 <- Arima(train.ts, order=c(0,1,1), method="ML")
summary(model1)## Series: train.ts
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## 0.1738
## s.e. 0.0868
##
## sigma^2 = 36.88: log likelihood = -363.68
## AIC=731.36 AICc=731.47 BIC=736.81
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1779347 6.019038 4.613207 -0.9146295 8.398176 0.9817368
## ACF1
## Training set 0.007439296lmtest::coeftest(model1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.173836 0.086835 2.0019 0.04529 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pada model ARIMA(0,1,1) didapatkan nilai AIC sebesar 731.36 dan seluruh parameter signifikan pada taraf nyata 5%.
ARIMA(0,1,3)
model2 <- Arima(train.ts, order=c(0,1,3), method="ML")
summary(model2)## Series: train.ts
## ARIMA(0,1,3)
##
## Coefficients:
## ma1 ma2 ma3
## 0.1677 0.0346 -0.1337
## s.e. 0.0954 0.1041 0.0885
##
## sigma^2 = 36.62: log likelihood = -362.3
## AIC=732.6 AICc=732.97 BIC=743.51
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.209285 5.944427 4.577686 -1.013171 8.292235 0.9741775 0.0096606lmtest::coeftest(model2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.167749 0.095422 1.7580 0.07875 .
## ma2 0.034556 0.104097 0.3320 0.73992
## ma3 -0.133679 0.088546 -1.5097 0.13112
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pada model ARIMA(0,1,3) didapatkan nilai AIC sebesar 732.60 dan tidak ada parameter yang signifikan pada taraf nyata 5%.
ARIMA(3,1,0)
model3 <- Arima(train.ts, order=c(3,1,0), method="ML")
summary(model3)## Series: train.ts
## ARIMA(3,1,0)
##
## Coefficients:
## ar1 ar2 ar3
## 0.1815 0.0054 -0.1843
## s.e. 0.0935 0.0951 0.0932
##
## sigma^2 = 36.18: log likelihood = -361.64
## AIC=731.29 AICc=731.66 BIC=742.2
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.2285709 5.908759 4.552318 -1.057525 8.23548 0.9687791
## ACF1
## Training set -0.007079369lmtest::coeftest(model3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.1815376 0.0934931 1.9417 0.05217 .
## ar2 0.0053613 0.0950847 0.0564 0.95504
## ar3 -0.1842896 0.0931874 -1.9776 0.04797 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pada model ARIMA(3,1,0) didapatkan nilai AIC sebesar 731.29 dan di antara tiga parameter hanya satu yang signifikan pada taraf nyata 5%, yaitu AR(3).
Pemilihan model
Model <- c("ARIMA(0,1,1)", "ARIMA(0,1,3)", "ARIMA(3,1,0)")
AIC <- c(model1$aic, model2$aic, model3$aic)
BIC <- c(model1$bic, model2$bic, model3$bic)
Sig.Coef <- c("1 dari 1", "0 dari 3", "1 dari 3")
tabel <- data.frame(Model, AIC, BIC, Sig.Coef)
kableExtra::kable(tabel, align = c(rep('c', times = 4)),
col.names = c("Model", "AIC", "BIC", "Parameter Signifikan"), digits = 2)| Model | AIC | BIC | Parameter Signifikan |
|---|---|---|---|
| ARIMA(0,1,1) | 731.36 | 736.81 | 1 dari 1 |
| ARIMA(0,1,3) | 732.60 | 743.51 | 0 dari 3 |
| ARIMA(3,1,0) | 731.29 | 742.20 | 1 dari 3 |
Pada tabel di atas, nilai AIC masing-masing model tidak berbeda signifikan. Model yang seluruh parameternya signifikan adalah ARIMA(0, 1, 1). Selain itu, model tersebut memiliki nilai BIC terkecil dibandingkan dengan model tentatif lainnya. Dengan demikian, model terbaik sementara adalah ARIMA(0, 1, 1).
Overfitting
ARIMA(1,1,1)
model1a <- Arima(train.ts, order=c(1,1,1), method="ML")
summary(model1a)## Series: train.ts
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.1284 0.0573
## s.e. 0.3070 0.2998
##
## sigma^2 = 37.15: log likelihood = -363.6
## AIC=733.2 AICc=733.42 BIC=741.38
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1703179 6.014652 4.61108 -0.8715949 8.392826 0.9812841
## ACF1
## Training set -0.000867559lmtest::coeftest(model1a)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.128449 0.306989 0.4184 0.6756
## ma1 0.057345 0.299762 0.1913 0.8483ARIMA(0,1,2)
model1b <- Arima(train.ts, order=c(0,1,2), method="ML")
summary(model1b)## Series: train.ts
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## 0.2084 0.0865
## s.e. 0.1004 0.1119
##
## sigma^2 = 37.01: log likelihood = -363.39
## AIC=732.77 AICc=732.99 BIC=740.95
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1567079 6.002836 4.604632 -0.7854681 8.377972 0.979912
## ACF1
## Training set -0.01426636lmtest::coeftest(model1b)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.208433 0.100398 2.0761 0.03789 *
## ma2 0.086505 0.111916 0.7730 0.43955
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model4c <- c("ARIMA(0,1,1)", "ARIMA(1,1,1)", "ARIMA(0,1,2)")
AIC2 <- c(model1$aic, model1a$aic, model1b$aic)
BIC2 <- c(model1$bic, model1a$bic, model1b$bic)
Sig.Coef2 <- c("1 dari 1", "0 dari 2", "1 dari 2")
tabel2 <- data.frame(Model4c, AIC2, BIC2, Sig.Coef2)
kableExtra::kable(tabel2, align = c(rep('c', times = 4)), col.names = c("Model", "AIC", "BIC", "Parameter Signifikan"), digits = 2)| Model | AIC | BIC | Parameter Signifikan |
|---|---|---|---|
| ARIMA(0,1,1) | 731.36 | 736.81 | 1 dari 1 |
| ARIMA(1,1,1) | 733.20 | 741.38 | 0 dari 2 |
| ARIMA(0,1,2) | 732.77 | 740.95 | 1 dari 2 |
Berdasarkan hasil uji signifikan dua model overfitting di atas, keduanya menunjukkan nilai AIC dan BIC yang lebih besar daripada model ARIMA(0,1,1) yang merupakan model terbaik sementara. Selain itu, terdapat parameter yang tidak signifikan pada kedua model overfitting ini. Dengan demikian, model terbaik yang dipilih untuk analisis selanjutnya adalah ARIMA(0,1,1).
Analisis Sisaan
Uji diagnostik model ARIMA dilakukan untuk melihat kenormalan dan kebebasan sisaan dari model yang dimiliki. Uji Shapiro-Wilk digunakan untuk mengetahui kenormalan sisaan, Uji Box-Ljung digunakan untuk mengetahui kebebasan sisaan, dan Uji t Satu Sampel digunakan untuk mengetahui apakah nilai tengah sisaan sama dengan nol atau tidak. Uji diagnostik dilakukan pada model terbaik yang sudah dipilih, yaitu model ARIMA(0,1,1).
Prosedur Eksploratif
sisaan <- model3$residuals
par(mfrow=c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(c(1:length(sisaan)),sisaan)
acf(sisaan)
pacf(sisaan)
Uji Formal
1. Sisaan Menyebar Normal
Uji dilakukan dengan Shapiro-Wilk Normality Test dengan hipotesis uji sebagai berikut:
H0 : Sisaan menyebar normal
H1 : Sisaan tidak menyebar normal
shapiro.test(sisaan)##
## Shapiro-Wilk normality test
##
## data: sisaan
## W = 0.97328, p-value = 0.02204Berdasarkan Shapiro-Wilk Normality Test, didapatkan p-value (0.02204) < α (0.05) sehingga H0 ditolak. Artinya, tidak cukup bukti untuk menyatakan bahwa sisaan menyebar normal pada taraf nyata 5%.
2. Sisaan Saling Bebas
Uji dilakukan dengan Box-Ljung Test dengan hipotesis uji sebagai berikut:
H0 : Sisaan saling bebas (tidak ada autokorelasi)
H1 : Sisaan tidak saling bebas (ada autokorelasi)
Box.test(sisaan, type = "Ljung") ##
## Box-Ljung test
##
## data: sisaan
## X-squared = 0.0058651, df = 1, p-value = 0.939Berdasarkan Box-Ljung Test, didapatkan p-value (0.939) > α (0.05) sehingga H0 tidak ditolak. Artinya, cukup bukti untuk menyatakan bahwa sisaan saling bebas atau tidak terjadi autokorelasi pada taraf nyata 5%.
3. Nilai Tengah Sisaan Sama dengan Nol
Uji dilakukan dengan One Sample t-test dengan hipotesis uji sebagai berikut:
H0 : Nilai tengah sisaan sama dengan nol
H1 : Nilai tengah sisaan tidak sama dengan nol
t.test(sisaan, mu = 0, conf.level = 0.95) ##
## One Sample t-test
##
## data: sisaan
## t = -0.41152, df = 113, p-value = 0.6815
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.3289844 0.8718426
## sample estimates:
## mean of x
## -0.2285709Berdasarkan t-test, didapatkan p-value (0.6815) > α (0.05) sehingga H0 tidak ditolak. Artinya, cukup bukti untuk menyatakan bahwa nilai tengah sisaan sama dengan nol pada taraf nyata 5%.
Forecasting
ramalan <- forecast(model1, h = 17)
data.ramalan <- ramalan$mean
plot(ramalan)
ts.plot(data$Terakhir, main = "Perbandingan Data Aktual dan\n Hasil Forecast")
lines((forecast(model1, h = 17)$mean), col = "maroon")
perbandingan <- matrix(data = c(test.ts[1:17], data.ramalan[1:17]), nrow = 17, ncol = 2)
colnames(perbandingan) <-c ("Aktual","Hasil Forecast")
perbandingan## Aktual Hasil Forecast
## [1,] 73.95 74.65687
## [2,] 68.50 74.65687
## [3,] 75.03 74.65687
## [4,] 83.57 74.65687
## [5,] 66.18 74.65687
## [6,] 75.21 74.65687
## [7,] 88.15 74.65687
## [8,] 95.72 74.65687
## [9,] 100.28 74.65687
## [10,] 104.69 74.65687
## [11,] 114.67 74.65687
## [12,] 105.76 74.65687
## [13,] 98.62 74.65687
## [14,] 89.55 74.65687
## [15,] 79.49 74.65687
## [16,] 86.53 74.65687
## [17,] 91.40 74.65687Berdasarkan plot dan nilai perbandingan nilai data aktual dan hasil ramalan, terjadi perbedaan nilai yang cukup signifikan. Data aktual cenderung meningkat dengan rentang 40 sampai dengan 63, tetapi data forecasting cenderung naik turun pada rentang 28 sampai dengan 40 sehingga perbandingan kedua data terbilang cukup jauh.
Perbandingan Nilai Akurasi
kableExtra::kable(accuracy(ramalan, test.ts),
align = c(rep('c', times = 8)), digits = 3)| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | -0.178 | 6.019 | 4.613 | -0.915 | 8.398 | 0.982 | 0.007 | NA |
| Test set | 13.420 | 19.072 | 15.224 | 13.179 | 15.856 | 3.240 | 0.761 | 1.871 |
Dari nilai keakuratan di atas, terlihat bahwa nilai error pada data testing jauh lebih besar daripada data training. Diambil contoh nilai MAPE pada data training sebesar 7.7926, sedangkan pada data testing sebesar 47.202. Perbedaan kedua nilai ini sangat signifikan. Hal tersebut mengindikasikan bahwa model melakukan prediksi dengan sangat baik pada data training dibandingkan dengan data testing.
Perbandingan Metode Terbaik
DES.acc <- accuracy(ramal_des$mean, test.ts)
ARIMA.acc <- accuracy(ramalan, test.ts)
best.acc <- rbind(DES.acc, ARIMA.acc[-1,-6])
rownames(best.acc) <- c("DES (α = 0.9, β = 0.2)", "ARIMA(2,1,2)")
kableExtra::kable(best.acc,
align = c(rep('c', times = 7)),
digits = 3)| ME | RMSE | MAE | MPE | MAPE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|
| DES (α = 0.9, β = 0.2) | -4.229 | 12.625 | 10.230 | -6.229 | 11.937 | 0.692 | 1.401 |
| ARIMA(2,1,2) | 13.420 | 19.072 | 15.224 | 13.179 | 15.856 | 0.761 | 1.871 |
Berdasarkan tabel perbandingan nilai akurasi di atas, terlihat bahwa model ARIMA(2,1,2) menghasilkan nilai error yang lebih kecil dibandingkan nilai error yang dihasilkan oleh metode DES, baik RMSE, MAE, MAPE, maupun yang lainnya. Dengan demikian, dapat disimpulkan bahwa model ARIMA(2,1,2) merupakan metode terbaik untuk meramalkan harga minyak mentah WTI bulanan berdasarkan harga sepuluh tahun terakhir.