Crude Oil Price Forecasting based on ARIMA Model
Untuk memperoleh minyak bumi, perlu dilakukan proses pengeboran. Minyak bumi yang ditemukan biasanya akan bercampur dengan gas alam. Minyak bumi yang telah dipisahkan dari gas alam berbentuk cairan kental hitam dan berbau disebut minyak mentah (crude oil). Pada proyek kali ini, akan dianalisis data rata-rata harga minyak mentah (crude oil) di negara India (dalam Rupee). Dengan mengambil data bulanan harga minyak mentah (crude oil) selama 239 bulan dari bulan Oktober 2000 sampai bulan Agustus 2020 yang diambil dari website Kaggle, dapat dilakukan proses prediksi harga minyak mentah (crude oil) untuk periode berikutnya dengan menggunakan metode ARIMA (Autoregressive Integrated Moving Average).
library yang digunakan:
library(forecast)
library(ggplot2)
library(tseries)
library(lmtest)
library(dplyr)
library(lubridate)Data Preparation
# Read Data
crude <- read.csv("GoldUP.csv")
head(crude)## Date Gold_Price Crude_Oil Interest_Rate USD_INR Sensex CPI
## 1 01-10-2000 4538 1455.51 8.0 46.31830 3711.02 37.23
## 2 01-11-2000 4483 1512.47 8.0 46.78361 3997.99 37.31
## 3 01-12-2000 4541 1178.11 8.0 46.74586 3972.12 36.98
## 4 01-01-2001 4466 1208.18 8.0 46.53603 4326.72 36.90
## 5 01-02-2001 4370 1267.18 7.5 46.51459 4247.04 36.73
## 6 01-03-2001 4269 1166.45 7.0 46.60935 3604.38 36.90
## USD_Index
## 1 116.65
## 2 115.24
## 3 109.56
## 4 110.52
## 5 112.01
## 6 117.37Untuk proyek ini, kita hanya akan berfokus pada variabel Crude_Oil.
# Data Cleaning
crude_clean <- crude %>%
select(Date, Crude_Oil) %>%
mutate(Date = dmy(Date))
head(crude_clean)## Date Crude_Oil
## 1 2000-10-01 1455.51
## 2 2000-11-01 1512.47
## 3 2000-12-01 1178.11
## 4 2001-01-01 1208.18
## 5 2001-02-01 1267.18
## 6 2001-03-01 1166.45Time Series
Agar dapat dianalisis, variabel Crude_oil diubah menjadi data berjenis time series dengan function ts(). Kemudian, kita menggunakan frequency=12 karena data yang kita gunakan merupakan data bulanan.
crude_ts <- ts(data = crude_clean$Crude_Oil, start = c(2000,10), frequency = 12)
crude_ts## Jan Feb Mar Apr May Jun Jul Aug Sep
## 2000
## 2001 1208.18 1267.18 1166.45 1203.33 1292.57 1267.73 1169.12 1216.80 1201.16
## 2002 925.63 972.86 1152.20 1244.02 1257.83 1199.06 1255.67 1301.25 1369.92
## 2003 1473.85 1569.55 1446.38 1210.97 1226.99 1304.45 1321.85 1363.37 1232.46
## 2004 1425.74 1418.36 1515.65 1480.99 1699.29 1617.64 1744.59 1950.09 1917.68
## 2005 1880.23 1957.60 2225.57 2215.05 2080.20 2348.85 2454.20 2700.01 2709.21
## 2006 2773.13 2646.56 2709.89 3055.21 3118.70 3145.13 3365.75 3341.97 2865.03
## 2007 2373.05 2541.99 2667.95 2742.35 2658.01 2779.56 2974.64 2862.54 3096.50
## 2008 3570.71 3710.34 4109.92 4353.72 5165.85 5631.69 5689.98 4919.52 4540.91
## 2009 2142.12 2059.46 2391.75 2517.58 2822.25 3303.55 3137.00 3461.09 3310.80
## 2010 3541.88 3463.32 3607.99 3745.92 3461.17 3479.77 3495.72 3531.14 3504.32
## 2011 4206.66 4450.50 4888.11 5159.40 4851.80 4747.57 4793.30 4550.04 4811.93
## 2012 5484.92 5540.26 5927.55 5888.49 5655.94 5083.60 5372.19 5848.75 5798.48
## 2013 5708.32 5784.94 5575.84 5375.04 5468.66 5817.69 6292.38 6836.67 6926.83
## 2014 6344.00 6529.28 6343.00 6329.60 6273.13 6471.05 6320.51 6092.56 5835.69
## 2015 2927.20 3398.51 3299.23 3611.03 3988.78 3915.30 3458.51 2973.16 3064.43
## 2016 2004.00 2117.58 2503.95 2708.63 3072.75 3208.66 2966.28 3004.16 3006.05
## 2017 3649.89 3647.03 3355.09 3365.04 3213.83 2975.01 3071.23 3194.46 3413.19
## 2018 4215.16 4085.16 4171.72 4516.93 4959.75 4879.75 4992.51 4942.53 5448.55
## 2019 4003.08 4352.94 4432.04 4761.33 4664.02 4149.67 4230.22 4102.97 4282.87
## 2020 4395.91 3811.78 2392.98 1603.02 2298.55 2987.46 3156.01 3243.75
## Oct Nov Dec
## 2000 1455.51 1512.47 1178.11
## 2001 995.46 897.07 887.42
## 2002 1331.61 1184.18 1342.56
## 2003 1316.78 1325.54 1366.36
## 2004 2146.34 1900.58 1717.00
## 2005 2608.22 2516.77 2576.49
## 2006 2633.08 2607.99 2722.26
## 2007 3238.90 3602.45 3530.67
## 2008 3535.80 2644.71 2010.86
## 2009 3461.10 3611.43 3491.65
## 2010 3629.80 3793.46 4065.66
## 2011 4916.61 5343.29 5486.72
## 2012 5476.43 5536.01 5526.62
## 2013 6497.77 6435.78 6534.28
## 2014 5280.84 4748.53 3806.55
## 2015 3055.95 2847.43 2435.40
## 2016 3290.46 3056.29 3572.84
## 2017 3574.53 3887.81 3930.99
## 2018 5648.69 4476.09 3822.07
## 2019 4069.14 4314.32 4509.77
## 2020Visualisasi Data
# Visualisasi data time series
ts.plot(crude_ts,ylab="Value",main="Plot Data Harga Minyak Mentah
(sebelum differencing)",col="orange",lwd=2)
legend("bottomright",c("Value"),cex=0.8,lty=5,text.font=2,col=c("orange"))
Insight: Dapat kita lihat dari plot yang terbentuk bahwa data crude_ts tidaklah memiliki pola dan tidak stasioner. Untuk itu, dilakukan proses differencing untuk mendapatkan data yang stasioner baik dari rata-rata maupun variansnya.
Plot ACF
# plot ACF
ggAcf(crude_ts)+ ggtitle("Plot ACF sebelum differencing")
Plot PACF
# Plot PACF
ggPacf(crude_ts)+ ggtitle("Plot PACF sebelum differencing")
Uji Stasioneritas Varians
# Uji Stasioneritas Varians
lambda=BoxCox.lambda(crude_ts)
lambda## [1] -0.1426602Dari output yang dihasilkan dapat dilihat nilai lambanya jauh dari 1, yaitu -0,1427. Dengan begitu perlunya dilakukan transformasi pada data time series yang telah dibentuk agar memiliki varians yang homogen dan stasioner.
# Transformasi
tf=((crude_ts^lambda)-1.1)/lambda
# Uji Stasioneritas Varians
lambda1=BoxCox.lambda(tf)
lambda1## [1] 1.006074Setelah dilakukan transformasi, output yang dihasilkan menjadi mendekati 1, yaitu 1,0061. Hal ini menunjukan bahwa nilai varians dari data yang digunakan telah homogen.
Uji Stasioneritas Rata-Rata
Hipotesis:
\(H_0\): Data time series tidak stasioner
\(H_1\): Data time series stasioner
Taraf signifikan:
α=5%
Kriteria uji: H0 ditolak, apabila nilai p-value < α, terima dalam hal lainnya.
adf.test(tf)##
## Augmented Dickey-Fuller Test
##
## data: tf
## Dickey-Fuller = -1.5399, Lag order = 6, p-value = 0.7693
## alternative hypothesis: stationaryKesimpulan: Dengan taraf signifikansi 5% dapat dinyatakan bahwa data yang digunakan tidaklah stasioner rata-ratanya sehingga perlu dilakukan proses differencing
# dilakukan 1 kali differencing
df.crude=diff(tf, 1)
adf.test(df.crude)##
## Augmented Dickey-Fuller Test
##
## data: df.crude
## Dickey-Fuller = -6.5566, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationarySetelah dilakukan proses differencing dilanjutkan dengan pengujian Augmented Dickey-Fuller. Dengan hipotesis dan taraf uji yang sama seperti sebelumnya dapat disimpulkan bahwa data time series telah stasioner rata-ratanya, ditandai dengan nilai p-value sebesar 0,01.
Visualisasi Data (setelah differencing)
# Visualisasi data setelah differencing
ts.plot(df.crude,ylab="Value",main="Plot Data Harga Minyak Mentah (setelah differencing)",col="brown",lwd=2)
legend("bottomright",c("Value"),cex=0.8,lty=5,text.font=2,col=c("brown"))
Plot ACF (setelah differencing)
ggAcf(df.crude)+ ggtitle("Plot ACF setelah differencing")
Plot PACF (setelah differencing)
ggPacf(df.crude)+ ggtitle("Plot PACF setelah differencing")
Insight:Setelah proses differencing dilakukan, terbukti bahwa data menjadi stationer seperti yang tampak pada plot di atas. Dapat dilihat pula bahwa plot ACF dan PACF yang dihasilkan sudah berpola cut off, dan keduanya menunjukan adanya cut off after lag 1 dan 2. Untuk itu, nilai p dan q yang diperkirakan adalah antara 1 dan 2.
ARIMA Model
Dengan menggunakan function arima(), kita dapat membuat model ARIMA dengan nilai p, d, dan q mengacu pada plot ACF dan PACF.
model1 <- arima(crude_ts, order = c(1,1,0))
model2 <- arima(crude_ts, order = c(0,1,1))
model3 <- arima(crude_ts, order = c(1,1,1))
model4 <- arima(crude_ts, order = c(2,1,1))
model5 <- arima(crude_ts, order = c(2,1,2))
model6 <- arima(crude_ts, order = c(0,1,2))list(coeftest(model1),
coeftest(model2),
coeftest(model3),
coeftest(model4),
coeftest(model5),
coeftest(model6))## [[1]]
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.379928 0.059745 6.3592 2.028e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.373143 0.053257 7.0064 2.445e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[3]]
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.20755 0.13107 1.5836 0.11329
## ma1 0.20921 0.12535 1.6690 0.09512 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[4]]
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.098924 0.193415 5.6817 1.334e-08 ***
## ar2 -0.383232 0.072637 -5.2760 1.320e-07 ***
## ma1 -0.691377 0.203154 -3.4032 0.000666 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[5]]
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.10120 0.80439 -0.1258 0.8999
## ar2 -0.16481 0.22383 -0.7363 0.4615
## ma1 0.52629 0.80014 0.6577 0.5107
## ma2 0.30568 0.22166 1.3790 0.1679
##
##
## [[6]]
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.436341 0.066569 6.5547 5.576e-11 ***
## ma2 0.133817 0.068678 1.9485 0.05136 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dari hasil pengujian signifikansi koefisien di atas, terlihat bahwa model yang memiliki hasil pengujian signifikansi koefisien signifikan yaitu model1, model2, model4, dan model6.
list(AIC(model1),AIC(model2),AIC(model4),AIC(model6))## [[1]]
## [1] 3364.016
##
## [[2]]
## [1] 3363.702
##
## [[3]]
## [1] 3360.508
##
## [[4]]
## [1] 3362.044Dari hasil list AIC dari model1, model2, model4, dan model6, terlihat bahwa yang memiliki nilai AIC yang paling kecil adalah model4. Oleh karena itu, model4 akan kita pakai untuk melakukan forecasting.
Untuk melihat ukuran kesalahan dari model4, dapat dilihat dengan menggunakan function accuracy()
accuracy(model4)## ME RMSE MAE MPE MAPE MASE
## Training set 6.603373 276.3122 205.7611 -0.03551758 6.695095 0.9569887
## ACF1
## Training set 0.001276973Uji Asumsi
Normalitas
Hipotesis:
\(H_0\): Residual model 4 berdistribusi normal
\(H_1\): Residual model 4 tidak berdistribusi normal
Taraf signifikan:
α=5%
Statistik uji:
Kolmogorov-Smirnov test
Kriteria uji:
H0 ditolak, apabila nilai p-value < α, terima dalam hal lainnya.
# Uji Normalitas
r2=residuals(model4)
n2=length(r2)
mean2=mean(r2)
sd2=sd(r2)
res2=rnorm(n2,mean2,sd2)
cek.normalitas=ks.test(r2,res2)
cek.normalitas##
## Two-sample Kolmogorov-Smirnov test
##
## data: r2 and res2
## D = 0.066946, p-value = 0.6578
## alternative hypothesis: two-sidedKesimpulan:
Dengan taraf signifikansi 5% dapat dinyatakan bahwa residual dari model berdistribusi normal.
Autokorelasi
Hipotesis:
\(H_0\): Tidak terjadi autokorelasi pada residual model 4
\(H_1\): Terjadi autokorelasi pada residual model 4
Taraf signifikan:
α=5%
Statistik uji:
Box-Ljung test
Kriteria uji:
H0 ditolak, apabila nilai p-value < α, terima dalam hal lainnya.
cek.WNA=Box.test(r2,lag=1,type=c("Ljung-Box")) ; cek.WNA##
## Box-Ljung test
##
## data: r2
## X-squared = 0.00039464, df = 1, p-value = 0.9842Kesimpulan:
Dengan taraf signifikansi 5% dapat dinyatakan bahwa tidak terjadi autokorelasi pada residual model 4.
Heteroskedastisitas
Hipotesis:
\(H_0\): Residual bersifat homoskedastik
\(H_1\): Residual bersifat heteroskedastik
Taraf signifikan:
α=5%
Statistik uji:
Box-Ljung test
Kriteria uji:
H0 ditolak, apabila nilai p-value < α, terima dalam hal lainnya.
h2=r2^2
cek.heteros=Box.test(h2,lag=1,type=c("Ljung-Box")) ;
cek.heteros##
## Box-Ljung test
##
## data: h2
## X-squared = 2.5347, df = 1, p-value = 0.1114Kesimpulan:
Dengan taraf signifikansi 5% dapat dinyatakan bahwa residual dari model 4 bersifat homoskedastik.
Forecasting
# Forecast
peramalan <- predict(model4,n.ahead=12)
peramalan$pred## Jan Feb Mar Apr May Jun Jul Aug
## 2020
## 2021 3145.844 3138.795 3136.913 3137.546 3138.964 3140.278 3141.180 3141.667
## Sep Oct Nov Dec
## 2020 3248.355 3219.790 3186.635 3161.147
## 2021# Visualisasi
crude_ts %>%
autoplot()+
autolayer(peramalan$pred, series = "forecast")
Kesimpulan
- Model untuk memprediksi rata-rata harga crude oil di India (dalam satuan rupee) untuk bulan September 2020 dan 11 bulan selanjutnya menggunakan data 239 bulan sebelumnya adalah ARIMA(2,1,1).
- Dengan menggunakan model tersebut, didapatkan nilai AIC sebesar 3360.508 dan nilai error MAPE sebesar 6.7%.
- Seluruh asumsi terpenuhi, mulai dari asumsi normalitas, asumsi non-autokorelasi, hingga asumsi homoskedastisitas.