Forecasting Bitcoin Price based on Time Series Model
Praneetha Durage
5 June 2019
Introduction
The purpose of this analysis is, forecasting bitcoin price based on time series model. Bitcoin is the very first cryptocurrency found in 2009. There are many arguments still prevail in financial market how to value a cryptocurrency as it behaves like a fiat currency and a asset. In this analysis I try to find a reasonable time series model by analysing historical behaviour of the currency. The techniques used here are Auto-Regressive Moving Average (ARIMA) and Generalised Auto-Regressive Conditional Heteroskedasticity (GARCH) in finding the best model that fits bitcoin prices. We use historical daily closed prices of the market for 27/04/2013 to 24/02/2019 published by coinmarketcap.com.
Discriptive Analysis
Before starting the modeling techniques, we visually analyse the dataset in order to identify the possible trend or seasonality, visible auto-correlation structure, any turning points or intervention to the market, whether the volatility of the prices changes over times by using data visualisation methods (graphs, charts, etc.).
# Load required packages
library(TSA)
library(fUnitRoots)
library(forecast)
library(lmtest)
library(fGarch)
library(rugarch)
library(tseries)bitcoin <- read.csv("bitcoin.csv")
# Removing non-numeric characters in the series
bitcoin$Close <- gsub('[^a-zA-Z0-9.]', '', bitcoin$Close)
bitcoin$Close <- (as.numeric(bitcoin$Close))
bitcoin.ts <- ts(as.vector(bitcoin$Close),start=c(2013,117),frequency = 365)
plot(bitcoin.ts, xlim=c(2013,2019), xlab="Time", ylab="Price in USD",
main="Bitcoin Price from 27/4/2013 to 24/02/2019")
Distribution of Bitcoin Prices
Based on Figure 1, there is no visible deterministic trend. The seasonality and auto-correlation structure is not obvious as the volatility of prices is very high after 2017 compared to the rest of the period. There is no visible turning points or intervention to the market from outside forces. We use continuously compounding return series of bitcoin to identify the changing variance of the prices over time or possibility of volatility clustering in Figure 2.
We can see the highest volatility is presented around the prices of 2014 and after till 2017 the volatility of prices are quite steady and beyond 2017 it shows high volatility again.
r.bitcoin=diff(log(bitcoin.ts))*100 # Continuous compound return
plot(r.bitcoin, xlab="Time", ylab="Return in %", xlim=c(2013,2019), main="Bitcoin Return")
Distribution of Bitcoin Return
Even though the auto-correlation structure is visibly not clear, we use Figure 3 below to identify evidence of auto-correlation in order to go for ARIMA models. As per Figure 3, auto-correlation of lag 1 is very strong (0.9976471).
plot(y=bitcoin.ts,x=zlag(bitcoin.ts), ylab="Bitcoin Price of Today (in $)",
xlab="Bitcoin Price of Yesterday (in $)", main="Plot of Autocorrelation for Lag 1")
Serial Correlation
x=zlag(bitcoin.ts)
index=2:length(x)
cor(bitcoin.ts[index],x[index])## [1] 0.9976471ARIMA Model
As it is hard to have a deterministic trend in long-run for financial data, we are not going to look for trend models here and directly jumping in to ARIMA models in order to capture the auto-correlation structure of historical prices. When using ARIMA techniques and maximum likelihood model estimation we need to have a stationary and normally distributed series. First differenced log series of bitcoin resulted continuously compounding return of bitcoin. We see in below chunk, the return series is stationary as the ADF test is significant at 5% level (P-Value < 0.01) rejecting the null hypothesis of presenting a unit root in process characteristic equation. However, Q-Q normal plot for return in Figure 4 with fat tails and significant Shapiro-Wilk normality test confirm the series is not normally distributed.
order = ar(diff(r.bitcoin))$order
adfTest(r.bitcoin, lags = order, title = NULL,description = NULL)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 32
## STATISTIC:
## Dickey-Fuller: -7.4259
## P VALUE:
## 0.01
##
## Description:
## Tue Jan 28 12:25:58 2020 by user: smqqnorm(r.bitcoin, main="Q-Q Normal Plot of Bitcoin Return")
qqline(r.bitcoin)
Q-Q Plot for Bitcoin Return
shapiro.test(r.bitcoin)##
## Shapiro-Wilk normality test
##
## data: r.bitcoin
## W = 0.88509, p-value < 2.2e-16Model Specification
Even though it’s not visibly clear that there is a seasonal trend on data, we use seasonal differenced return series of bitcoin to find out the possibility of having seasonality. Figure 5 below presents ACF and PACF of seasonal differenced series. ACF indicates possibility of one seasonal MA and PACF shows evidence for two seasonal AR terms. However, prevailing packages (either TSA or Forecast) do not provide functions to fit SARIMA models with lags over 350. Therefore, we stick our analysis only on ARIMA models.
To determine the orders of ordinary AR and MA terms we use ACF, PACF, EACF and Bayesian Information Criterion (BIC) table of seasonal differenced bitcoin return series with lags up to 365 (Figure 5).
As per Figure 5, it is possible to have higher order AR and MA terms. Therefore, we consider all possible models up to 6 AR and 6 MA terms
# Getting seasonal difference
diff.r.bitcoin=diff(r.bitcoin, lag=365)
par(mfrow=c(1,2))
acf(diff.r.bitcoin, lag.max = 730, main="Seasonal Diff. Return Series")
pacf(diff.r.bitcoin, lag.max = 730, main="Seasonal Diff. Return Series")
ACF anf PACF of Seasonal Differenced Return Series
eacf(diff.r.bitcoin)## 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 x o o o x o o o
## 2 o o o o o x o o o o x o o o
## 3 x o x o o x o o o o x o o o
## 4 o o o x o x o o o o x o o o
## 5 x x o x x x o o o o x o o o
## 6 x x o x x x o o o o o x o o
## 7 x x o x x x o o o o o x o oAs per the above r output with eacf, we do not get clear cut off for selecting possible models (As eacf is more fuzzy and no clear vertices).
res = armasubsets(y=diff.r.bitcoin,nar=7,nma=7,y.name='test',ar.method='ols')## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found## Reordering variables and trying again:plot(res)
BIC Table
As per BIC table (Figure 6), the possible models are ARIMA(6,1,0), ARIMA(5,1,6) and ARIMA(5,1,7). However, models with higher orders are difficult to estimate as convergence problems arise. Therefore, we restrict our candidate models set to {ARIMA(1,1,0), ARIMA(0,1,1), ARIMA(2,1,1), ARIMA(1,1,2), ARIMA(2,1,2), ARIMA(2,1,3), ARIMA(3,1,2),ARIMA(3,1,3),ARIMA(4,1,3), ARIMA(3,1,4), ARIMA(4,1,4), ARIMA(5,1,4), ARIMA(4,1,5), ARIMA(5,1,5), ARIMA(5,1,6), ARIMA(6,1,5), ARIMA(6,1,6)}.
Model Estimation
When estimating all possible models together we use myCandidate() introduced by Yong Kai, Wong. Below chunk dedicates for model estimation.
source('C:/Users/sampa/Desktop/Utility_functions/TSHandy.r')
modelList <- list(c(1,1,0), c(0,1,1),c(2,1,1),c(1,1,2),c(2,1,2),c(2,1,3),c(3,1,2),c(3,1,3),
c(4,1,3),c(3,1,4),c(4,1,4),c(5,1,4), c(4,1,5),c(5,1,5),c(5,1,6),c(6,1,5),c(6,1,6))
modelEstimation_ML <- myCandidate(log(bitcoin.ts), orderList = modelList, methodType = "ML")
modelEstimation_CSS <- myCandidate(log(bitcoin.ts), orderList = modelList, methodType = "CSS")
modelEstimation_ML$IC## p d q AIC AICc BIC
## 15 5 1 6 -7332.099 -7331.952 -7264.139
## 16 6 1 5 -7331.836 -7331.689 -7263.875
## 17 6 1 6 -7329.465 -7329.293 -7255.841
## 13 4 1 5 -7325.694 -7325.590 -7269.060
## 12 5 1 4 -7325.623 -7325.519 -7268.989
## 14 5 1 5 -7325.438 -7325.313 -7263.140
## 11 4 1 4 -7323.862 -7323.777 -7272.891
## 6 2 1 3 -7321.132 -7321.093 -7287.152
## 7 3 1 2 -7321.132 -7321.092 -7287.152
## 8 3 1 3 -7319.343 -7319.291 -7279.700
## 9 4 1 3 -7302.444 -7302.376 -7257.137
## 10 3 1 4 -7302.287 -7302.219 -7256.980
## 2 0 1 1 -7299.666 -7299.661 -7288.339
## 1 1 1 0 -7299.665 -7299.659 -7288.338
## 3 2 1 1 -7296.181 -7296.162 -7273.527
## 4 1 1 2 -7296.156 -7296.138 -7273.503
## 5 2 1 2 -7294.316 -7294.287 -7265.999By comparing significance of the estimates, AIC and BIC values we choose ARIMA(5,1,6) as the best model. In model estimation, we use both Maximum Likelihood Estimation (MLE) and Conditional Some of Square (CSS) methods.
In below chunk, we present the model estimation detail for the best model. The first estimation from MLE and the second from CSS.
m_ml=arima(log(bitcoin.ts), order=c(5,1,6), method="ML")
m_css=arima(log(bitcoin.ts), order=c(5,1,6), method="CSS")
coeftest(m_ml)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.121116 0.059491 -2.0359 0.04176 *
## ar2 0.424654 0.064296 6.6047 3.983e-11 ***
## ar3 -0.525592 0.038573 -13.6259 < 2.2e-16 ***
## ar4 0.087579 0.047489 1.8442 0.06516 .
## ar5 0.852951 0.051592 16.5326 < 2.2e-16 ***
## ma1 0.118409 0.062934 1.8815 0.05991 .
## ma2 -0.429570 0.063943 -6.7180 1.842e-11 ***
## ma3 0.524959 0.043123 12.1734 < 2.2e-16 ***
## ma4 -0.047808 0.051003 -0.9374 0.34857
## ma5 -0.818237 0.053418 -15.3177 < 2.2e-16 ***
## ma6 0.047342 0.026414 1.7923 0.07308 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(m_css)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.250360 0.107598 -2.3268 0.01998 *
## ar2 0.444826 0.189709 2.3448 0.01904 *
## ar3 -0.398113 0.077828 -5.1153 3.132e-07 ***
## ar4 -0.085168 0.118976 -0.7158 0.47409
## ar5 0.581019 0.098896 5.8750 4.228e-09 ***
## ma1 0.245828 0.108768 2.2601 0.02381 *
## ma2 -0.456099 0.194220 -2.3484 0.01886 *
## ma3 0.392388 0.082941 4.7309 2.235e-06 ***
## ma4 0.141236 0.133919 1.0546 0.29159
## ma5 -0.524782 0.101445 -5.1731 2.303e-07 ***
## ma6 0.054730 0.038124 1.4356 0.15112
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Goodness of Fit of the Model
We use residuals generated through the best model to diagnose the possibility of violating the underlying model assumptions. The residuals resulted from the model ARIMA(5,1,6) should be a normally distributed white noise series. In order to check this assumption we use residual.analysis() introduced by Dr.Haydar Demirhan.
source('C:/Users/sampa/Desktop/Utility_functions/residual.analysis.r')
residual.analysis(modelEstimation_ML$model[[15]], std = TRUE,start = 1)##
## Shapiro-Wilk normality test
##
## data: res.model
## W = 0.89292, p-value < 2.2e-16
Residual Analysis for ARIMA(5,1,6)
As per Figure 7, standardised residuals move around zero and report quite large values around 2014. It seems the volatility of residuals are not constant over the time. Histogram for residuals is symmetric across zero, however, Q-Q plot shows flat tails and Shapiro-Wilk normality test confirms the residuals are not normal (p-value < 2.2e-16). ACF, PACF and Ljung-Box test together confirm there is no auto-correlation remains among residuals. In order to check the adequacy of the model we fit ARIMA(6,1,7) and check the significance of the estimates. We can see in ARIMA(6,1,7) many estimates are no longer significant and suggests that ARIMA(5,1,6) is sufficient.
m_616=arima(log(bitcoin.ts), order=c(6,1,7), method="CSS")
coeftest(m_616)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.590284 0.180735 3.2660 0.0010907 **
## ar2 0.028782 0.340432 0.0845 0.9326217
## ar3 -0.099355 0.370611 -0.2681 0.7886350
## ar4 0.239830 0.188222 1.2742 0.2025974
## ar5 0.168719 0.248206 0.6798 0.4966588
## ar6 -0.033442 0.201692 -0.1658 0.8683090
## ma1 -0.598170 0.181388 -3.2977 0.0009747 ***
## ma2 -0.036073 0.342319 -0.1054 0.9160747
## ma3 0.117964 0.371289 0.3177 0.7507006
## ma4 -0.219025 0.188759 -1.1603 0.2459095
## ma5 -0.134728 0.250625 -0.5376 0.5908753
## ma6 0.081483 0.196391 0.4149 0.6782148
## ma7 -0.060941 0.024998 -2.4379 0.0147741 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We noticed the volatility of the residuals swings overtime and provides evidence for heteroskedasticity or volatility clustering. We use McLeod-Li test to find-out the presence of volatility clustering more clearly.
model_residuals = modelEstimation_ML$model[[15]]$residuals
McLeod.Li.test(y=model_residuals, main="MCleod-Li Test for ARIMA(5,1,6) Residuals")
MCleod-Li Test for ARIMA(5,1,6) Residuals
In Figure 8, with all the significant lags it confirms the presence of heteroskedasticity. We can use Generalised Autoregressive Conditional Heteroskedasticity (GARCH) models to capture heteroskedastic errors in a model along with ARIMA concept.
GARCH Model
By using GARCH concept with ARIMA modeling we try to estimate conditional mean and variance of future values based on current and past values. Package rugarch provides sufficient functions to estimate ARIMA-GARCH models.
Specification of the Model
In order to identify the GARCH orders we use ACF, PACF and EACF of squared and absolute residuals resulted from ARIMA(5,1,6) model.
par(mfrow=c(1,2))
abs.res = abs(model_residuals)
acf(abs.res, main="Abs(Residuals) of ARIMA(5,1,6)")
pacf(abs.res,main="Abs(Residuals) of ARMA(5,1,6)")
ACF and PACF of Absolute Residuals of ARIMA(5,1,6)
eacf(abs.res)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x x x x
## 1 x x o o x o o o o x o o o o
## 2 x x o o o o o o o o o o o o
## 3 x x x o o o o o o o o o o o
## 4 x x x x o o o o o o o o o o
## 5 x x x x x o o o o o o o o o
## 6 x o x x x x o o o o o o o o
## 7 x x x x x x x o o o o o o oACF or PACF in Figure 9, does not show clear cut-off to identify the orders of GARCH. As per EACF output vertices suggest the possibility of GARCH (2,2), GARCH(3,2) and GARCH(3,3).
In the same way we check ACF, PACF and EACF of squared residuals of ARIMA(5,1,6) in below chunk.
par(mfrow=c(1,2))
sqr.res = model_residuals^2
acf(sqr.res, main="Squared Residuals of ARIMA(5,1,6)")
pacf(sqr.res,main="Squared Residuals of ARMA(5,1,6)")
ACF and PACF of Squared Residuals of ARIMA(5,1,6)
eacf(sqr.res)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x x x x
## 1 x x x x x x o o o x o o x o
## 2 x x x o o x o o o x o o x o
## 3 o x x o o x o o o o o o x o
## 4 o x x x o o o o o o o o x o
## 5 x x x x x x o o o o o o o o
## 6 x x x x x x o o o o o o o o
## 7 x x x x x x o o o o o o o oEither ACF or PACF in Figure 10 does not provide clear cut-off for GARCH orders. EACF for squared residuals is also too fuzzy and hard to find vertices. Therefore, we decide to go with all possible Garch models up to order 3.
Here is the candidates garch model set we considered. {GARCH(1,1), GARCH(1,2), GARCH(2,1), GARCH(2,2), GARCH(2,3),GARCH(3,2),GARCH(3,3)}
After analysing each garch model with ARIMA(5,1,6) we found GARCH(1,2) is the best model based on significance of the parameters and AIC, BIC values. Below chunk dedicates to see the parameter estimation of ARIMA(5,1,6) with GARCH(1,2).
model1<-ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 2)),
mean.model = list(armaOrder = c(5, 6), include.mean = FALSE),
distribution.model = "norm")
model.fit<-ugarchfit(spec = model1, data = r.bitcoin, out.sample = 100)
model.fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,2)
## Mean Model : ARFIMA(5,0,6)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.33549 0.058463 -5.7384 0.000000
## ar2 0.51919 0.022382 23.1967 0.000000
## ar3 -0.64647 0.028742 -22.4918 0.000000
## ar4 -0.73991 0.031145 -23.7571 0.000000
## ar5 0.22739 0.039901 5.6989 0.000000
## ma1 0.37521 0.059124 6.3461 0.000000
## ma2 -0.50457 0.015234 -33.1219 0.000000
## ma3 0.62575 0.028181 22.2046 0.000000
## ma4 0.79617 0.024179 32.9286 0.000000
## ma5 -0.17736 0.052922 -3.3514 0.000804
## ma6 0.01905 0.017591 1.0829 0.278844
## omega 0.43303 0.093273 4.6426 0.000003
## alpha1 0.18240 0.022230 8.2054 0.000000
## beta1 0.22140 0.061234 3.6157 0.000300
## beta2 0.58838 0.057926 10.1574 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.33549 0.047571 -7.0524 0.000000
## ar2 0.51919 0.059602 8.7109 0.000000
## ar3 -0.64647 0.040332 -16.0285 0.000000
## ar4 -0.73991 0.020923 -35.3631 0.000000
## ar5 0.22739 0.053131 4.2799 0.000019
## ma1 0.37521 0.048095 7.8014 0.000000
## ma2 -0.50457 0.055597 -9.0754 0.000000
## ma3 0.62575 0.036868 16.9727 0.000000
## ma4 0.79617 0.024723 32.2039 0.000000
## ma5 -0.17736 0.048132 -3.6849 0.000229
## ma6 0.01905 0.017251 1.1042 0.269489
## omega 0.43303 0.234997 1.8427 0.065374
## alpha1 0.18240 0.043434 4.1996 0.000027
## beta1 0.22140 0.066304 3.3392 0.000840
## beta2 0.58838 0.073014 8.0585 0.000000
##
## LogLikelihood : -5507.312
##
## Information Criteria
## ------------------------------------
##
## Akaike 5.4434
## Bayes 5.4849
## Shibata 5.4433
## Hannan-Quinn 5.4586
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.111 0.14627
## Lag[2*(p+q)+(p+q)-1][32] 25.185 0.00000
## Lag[4*(p+q)+(p+q)-1][54] 36.105 0.01921
## d.o.f=11
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.458 0.4985
## Lag[2*(p+q)+(p+q)-1][8] 1.288 0.9501
## Lag[4*(p+q)+(p+q)-1][14] 3.616 0.9139
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.9414 0.500 2.000 0.3319
## ARCH Lag[6] 1.0551 1.461 1.711 0.7320
## ARCH Lag[8] 1.2399 2.368 1.583 0.8862
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 2.7063
## Individual Statistics:
## ar1 0.38112
## ar2 0.20921
## ar3 0.42311
## ar4 0.19501
## ar5 0.09263
## ma1 0.47405
## ma2 0.32762
## ma3 0.46353
## ma4 0.21793
## ma5 0.11817
## ma6 0.06654
## omega 0.14004
## alpha1 0.07865
## beta1 0.10529
## beta2 0.11045
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 3.26 3.54 4.07
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.1060 0.2689
## Negative Sign Bias 0.8147 0.4153
## Positive Sign Bias 1.1288 0.2591
## Joint Effect 2.0294 0.5663
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 298.7 3.766e-52
## 2 30 324.6 1.085e-51
## 3 40 333.6 1.862e-48
## 4 50 352.4 1.663e-47
##
##
## Elapsed time : 0.9862592Goodness of Fit of the Model
We use residuals generated from final model to check whether model assumptions are violated. The residuals should be a normally distributed random series with zero mean and constant variance. Again we use residual.analysis() to detect the residuals from ARIMA-GARCH model.
residual.analysis(model.fit,class="ARMA-GARCH",start=2)##
## Shapiro-Wilk normality test
##
## data: res.model
## W = 0.88715, p-value < 2.2e-16
Residual Analysis for ARIMA-GARCH Model
mean(model.fit@fit$residuals)## [1] 0.1695804As per Figure 11, the residuals of final model do not look good. The mean value of residual is 0.17 and histogram is not symmetric and shows very large negative and positive residuals. Q-Q plot shows fat tails and Shapiro-Wilk normality test confirms residuals are not normally distributed. ACF are PACF are quite alright. Ljung-Box test is good till lag 10. Based on all these residual analysis it looks like ARIMA(5,1,6) model with GARCH(1,2) does not have good fit. Therefore, we are not going to use this model for future forecasting.
Forecasting
Earlier we noticed ARIMA+GARCH model doesn’t have sufficient goodness of fit. However, here we compare fitted values from both models with original series.
fitted_ARIMA=exp(log(bitcoin.ts)-model_residuals)
plot(fitted_ARIMA, col=2, ylab="Bitcoin Value (in $)",
main="Arima(5,1,6) Fitted Values and Original Series")
lines(bitcoin.ts)
ARIMA(5,1,6) Fitted Values with Original Series
Figure 12 gives us the fitted series from ARIMA (5,1,6) (in red colour) and original series (in black colour). It’s very hard to identify two series separately as the fitted values are near identical to original values.
par(mfrow=c(1,2))
fitted_GARCH=exp(diffinv(x=model.fit@fit$fitted.values/100, xi=log(bitcoin.ts[1])))
fitted_GARCH=ts(fitted_GARCH, start=c(2013,117),frequency = 365)
plot(fitted_GARCH, col=2, ylab="Bitcoin Price (in $)",
main="ARIMA+GARCH with Original", type="l", ylim=c(130,20000))
lines(bitcoin.ts)
plot(fitted_GARCH, main="ARIMA+GARCH", ylab="Bitcoin Price (in $)", col=2)
Comparing ARIMA-GARCH Fit with Original Series
Figure 13 depicts fitted values (in red) from ARIMA-GARCH versus original bitcoin prices. In contrast to ARIMA model, fitted values from ARIMA-GARCH are extremely deviated from original values. It seems variation we noticed in original series around 2014 is over-captured by ARIMA_GARCH model.
As the ARIMA(5,1,6) looks good with fitted values, we use that model for forecasting future bitcoin prices. Figure 14 depicts the forecasted log(bitcoin) values with possible risk window.
model=Arima(log(bitcoin.ts),c(5,1,6))
plot(forecast(model,h=10), ylab="Log(Bitcoin Price)", xlab="Time")
Log(Bitcoin) Forecast from ARIMA(5,1,6)
bit_foreast <- read.csv("Bit_forecast.csv")# Getting observed out-of sample values
source('C:/Users/sampa/Desktop/Utility_functions/MASE.txt')# Source for MASE()
#Getting predicted values for next 10 days
pred_values=predict(model, n.ahead=10,newxreg = NULL, se.fit = TRUE)
#Calculating MASE for fitted series
MASE_fitted=MASE(bitcoin.ts,fitted_ARIMA)
MASE_fitted## $MASE
## MASE
## 1 1.003302#Calculating MASE for forecast series
MASE_forecast=MASE(bit_foreast$Close,exp(pred_values$pred))
MASE_forecast## $MASE
## MASE
## 1 1.348925Discussion and Conclusion
Our objective of this analysis was to find a reasonable time series model to predict future bitcoin prices. Here we considered two possible model structures namely ARIMA and ARIMA-GARCH. Even though the residuals resulted from ARIMA model do not have a constant variance over time, the fitted values from the model go close with original values. On the other hand, ARIMA-GARCH model captures the volatility structure of the data, however, fitted values are far away from original values. Therefore, we use ARIMA(5,1,6) model to forecast bitcoin prices for next 10 days. In order to gauge the accuracy of the forecast we used Mean Absolute Scaled Error (MASE). We noticed MASE over the fitted value is 0.997 and this value for forecasted data is 1.349. If we can achieve below 1 MASE value, the accuracy of our model is very high. However, achieving such a value is challenging.
The price of bitcoin is mainly driven by profit expectation (speculation) rather based on market fundamentals (intrinsic value). We believe ARIMA-GARCH model gives us the intrinsic value of bitcoin and ARIMA model gives the speculated values. ARIMA-GARCH values are more conservative and mainly driven by market fundamentals while ARIMA values are more vulnerable to changes.