Introduction
Exploratory Data Analysis
Models for Bitcoin Prices
- ARIMA Model
- ARFIMA model
Frequency and Trend Analysis
- Frequency Analysis
- Trend Analysis
Volatility Analysis
Conclusion
References

Introduction

Bitcoin is a type of decentralized cryptocurrency in which encryption techniques are used to govern the creation of units of currency and verify the transfer funds between multiple parties. It has gained popularity worldwide, motivating possible studies about how its price changes. Because Bitcoin is a decentralized market there is price volatility.

In this project, we aim to investigate the price and volatility of Bitcoin from January 1, 2013 to February 12, 2021, and answer the following questions:

Can we fit a model for Bitcoin prices and possibly predict future prices?
Are there any seasonal changes or trends for Bitcoin prices?
Can we model the volatility of Bitcoin prices?

We use historical Bitcoin price data on Investing.com for this project\(^{[1]}\). By observing this dataset and testing various statistical models, a better understanding can be gained on models that may be used to predict Bitcoin prices and variability.

Exploratory Data Analysis

Bitcoin closing price data from Jan 1 2013 to Feb 12 2021 was examined. In total there were 2,965 observations.

The time series plot above highlights notable events over the observed time period. The plot summary above shows that the Bitcoin price remains relatively flat until around the year 2017, when volatility and price sharply increase. The reason for this price increase is still disputed, however it is important to note that the popularity of Bitcoin increased substantially around this time period.

Given there is evidence of non-stationarity, non-stable variance and mean, we take the log and first difference of the time series. This is performed because the log is a known variance stabilizing transformation, and the first-difference helps stabilize the mean of a time series by reducing the trend in the closing price over the level of the time series. This results in the following time series plot, autocorrelation plot, and histogram:

data_ts_2013_log <- log(data_ts_2013)
data_ts_2013_log %>% diff() %>% ggtsdisplay(main="Log and First Difference of Bitcoin Closing Price ", plot.type = "hist", lag.max=20)

data_ts_2013_logdiff <- log(data_ts_2013) %>% diff()

The transformed time series, apart from a few time periods, appears to be stationary. The time periods when the time series is not stationary (e.g., beginning of 2014) are largely due to the high Bitcoin price volatility, even after logarithm transformation of the time series is performed.

Models for Bitcoin Prices

In this section we will fit various time series models to the Bitcoin data and identify which, if any, performs best.

ARIMA Model

We will first try fit data to ARIMA(p,1,q) model\(^{[3]}\), under the assumption that differencing has eliminated the trend of data. We will vary p between 0 and 5, and q between 0 and 6. The AIC table of the fitted 42 models are below:

aic_table <- function(data,P,Q){
table <- matrix(NA,(P+1),(Q+1))
for(p in 0:P) {
for(q in 0:Q) {
table[p+1,q+1] <- arima(data,order=c(p,1,q))$aic
}
}
dimnames(table) <- list(paste("AR",0:P, sep=""),paste("MA",0:Q,sep=""))
table
}
Low_aic_table <- aic_table(data_ts_2013_log,5,6)
require(knitr)
kable(Low_aic_table,digits=2)

	MA0	MA1	MA2	MA3	MA4	MA5	MA6
AR0	-8361.31	-8360.00	-8447.10	-8446.63	-8549.46	-8571.02	-8603.03
AR1	-8359.72	-8381.26	-8445.44	-8488.85	-8557.27	-8582.16	-8607.91
AR2	-8477.99	-8475.99	-8548.00	-8590.22	-8605.75	-8606.70	-8611.46
AR3	-8475.99	-8473.99	-8600.12	-8598.61	-8603.95	-8607.44	-8609.48
AR4	-8560.77	-8578.11	-8599.16	-8598.97	-8608.29	-8615.27	-8613.29
AR5	-8591.05	-8591.74	-8589.92	-8595.23	-8617.48	-8613.19	-8611.29

Since ARIMA(5,1,4) gives the lowest AIC value, we will fit an ARMIA(5,1,4).

model_arima514 <- arima(data_ts_2013_log, order=c(5,1,4))
model_arima514

## 
## Call:
## arima(x = data_ts_2013_log, order = c(5, 1, 4))
## 
## Coefficients:
##          ar1      ar2     ar3      ar4     ar5      ma1     ma2      ma3
##       0.3609  -1.2055  0.5988  -0.6193  0.1576  -0.3631  1.0424  -0.5047
## s.e.  0.0575   0.0671  0.0691   0.0530  0.0252   0.0559  0.0627   0.0510
##          ma4
##       0.6100
## s.e.  0.0415
## 
## sigma^2 estimated as 0.003176:  log likelihood = 4318.74,  aic = -8617.48

\(\sigma^2\) of fitted model is very close to 0, however all model coefficients appear to be statistical significant. A 95% confidence interval of the coefficients are below.

t(confint(model_arima514))

##              ar1       ar2       ar3        ar4       ar5        ma1       ma2
## 2.5 %  0.2482452 -1.337038 0.4634485 -0.7232660 0.1081041 -0.4727261 0.9195433
## 97.5 % 0.4735249 -1.074037 0.7341669 -0.5153401 0.2070309 -0.2534787 1.1653514
##               ma3       ma4
## 2.5 %  -0.6046087 0.5286688
## 97.5 % -0.4048694 0.6913478

None of the confidence intervals include 0. Next we verify the model is both causal and invertible. This can be easily done in R by checking the inverse roots. Because we are checking the inverse roots and not the “normal” roots, we want the inverse roots to lie within the unit circle to confirm the model is both causal and invertible.

autoplot(model_arima514, main = "Plotting the ARIMA(5,1,4) characteristic roots")

All the roots lie within the unit circle implying the model is both causal and invertible, however 2 of the AR coefficients and 2 of the MA coefficients are near the boundaries of the unit circle. Hence, in addition to fitting and checking the residuals of the ARIMA(5,1,4) model, we will also fit the smaller ARIMA(3,1,2) model.

The residuals diagnostic plots of the ARIMA(5,1,4) are shown below.

par(mfrow=c(1,2))
acf(model_arima514$residuals, main = "ARIMA(5,1,4) Autocorelation Plot")
qqnorm(model_arima514$residuals, main = "ARIMA(5,1,4) Q Q Plot")
qqline(model_arima514$residuals)

From the autocorrelation plot we observe that there does not appear to exist any serial correlation amongst the residuals, however the qqnorm shows heavy tails on both side, which violates the assumption that noises are normally distributed. We could conclude that ARIMA(5,1,4) does not meet all the assumptions of the arima model for the log-differenced Bitcoin data set.

Next, we will repeat this for the ARIMA(3,1,2) model. The model fit is as follows:

model_arima312 <- arima(data_ts_2013_log, order=c(3,1,2))
model_arima312

## 
## Call:
## arima(x = data_ts_2013_log, order = c(3, 1, 2))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2
##       0.8899  -0.7749  0.2099  -0.8949  0.6060
## s.e.  0.0490   0.0369  0.0210   0.0484  0.0411
## 
## sigma^2 estimated as 0.003203:  log likelihood = 4306.06,  aic = -8600.12

95% confidence intervals of the coefficients is shown below.

t(confint(model_arima312))

##              ar1        ar2       ar3        ma1       ma2
## 2.5 %  0.7938006 -0.8472752 0.1688024 -0.9898003 0.5255391
## 97.5 % 0.9860130 -0.7025577 0.2509432 -0.7999712 0.6864679

None of the confidence intervals of the coefficients include 0. Next, check if the ARIMA(3,1,2) model is both causal and invertible.

autoplot(model_arima312, main = "Plotting the ARIMA(3,1,2) characteristic roots")

All the inverse roots lie within the unit circle implying the model is both invertible and casual. Contrary to the ARIMA(5,1,4) model, none of the inverse roots lie on the boundaries of the unit circle. Next, we can check its’ residuals.

par(mfrow=c(1,2))
acf(model_arima312$residuals, main = "ARIMA(3,1,2) Autocorelation Plot")
qqnorm(model_arima312$residuals, main = "ARIMA(3,1,2) Q Q Plot")
qqline(model_arima312$residuals)

From the ACF plot we observe that autocorrelations are high at the beginning and decrease to be within the confidence interval. The qqnorm plot shows heavy tails on both side, which violates the assumption that noises are normally distributed.

The ARIMA(3,1,2) model appears to perform almost identical to the ARIMA(5,1,4) model in regards to the autocorrelation and qqplots. We can choose the simpler ARIMA(3,1,2) model as a reasonable model. Given the heavy tails observed in the Q-Q plots of the residuals of both models a time series model that can handle non-normal residuals may be needed.

ARFIMA model

First, in order to determine whether long-term effects are obvious, we will compare the ACF plot of the Bitcoin time series and the random walk time series.

library(arfima)
acf(log(data_ts_2013), 100, main="ACF plot of the data")

acf(cumsum(rnorm(1000)), 100, main="ACF plot of random walk") # compare to ACF of random walk

The plots show that the ACF of Bitcoin prices decays slower than the random walk when increasing lags. This suggest that there are possible long-memory effects. In order to take the possible long-memory effects into consideration, an ARFIMA model \(^{[4]}\) can be fitted. The summary is shown below. R examples for the ARFIMA model was found in the text “Time Series Analysis and Its Applications” by Shumway and Stoffer. \(^{[8]}\)

summary(data_ts_2013.fd <- arfima(log(data_ts_2013)))

## Note: only one starting point.  Only one mode can be found -- this is now the default behavior.
## Beginning the fits with 1 starting values.

## 
## Call:
##  
## arfima(z = log(data_ts_2013))
## 
## 
## Mode 1 Coefficients:
##                Estimate  Std. Error Th. Std. Err.     z-value Pr(>|z|)    
## d.f         4.99756e-01 6.32218e-07   1.43178e-02 7.90480e+05  < 2e-16 ***
## Fitted mean 3.61576e+00 2.95939e+00            NA 1.22179e+00  0.22179    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## sigma^2 estimated as 0.0134651; Log-likelihood = 6382.29; AIC = -12758.6; BIC = -12740.6
## 
## Numerical Correlations of Coefficients:
##             d.f  Fitted mean
## d.f         1.00 0.00       
## Fitted mean 0.00 1.00       
## 
## Theoretical Correlations of Coefficients:
##     d.f 
## d.f 1.00
## 
## Expected Fisher Information Matrix of Coefficients:
##     d.f 
## d.f 1.65

As shown from the summary, the estimate of \(d\) for the \(ARFIMA(0,d,0)\) model is about 0.5. The residuals for this model can also be plotted to show the correlations. We observe that the residuals are correlated, implying that the residuals are not white noises. Because the residuals are not white noises we can deduce that the ARFIMA model does not perform well using this bitcoin dataset.

innov = resid(data_ts_2013.fd)  
plot.ts(innov[[1]],ylab="Residuals")

acf(innov[[1]], main="ACF plot of residuals")

Frequency and Trend Analysis

Frequency Analysis

Using frequency analysis We can try to identify cycles for Bitcoin prices. The smoothed periodogram can be used to determine whether there are any cyclic patterns for Bitcoin prices. The default kernel smoothing will be used for this analysis.\(^{[5]}\)

a = spectrum(data_ts_2013, spans=c(30, 30), plot=TRUE, main="Smoothed Periodogram")

From the plot above, we see that the error bar covers any possible peak of the estimated spectral density. Therefore, there is no obvious seasonal behavior for the Bitcoin prices. There is still variability for the prices however, there are no cycles at predictable lengths. This could be contributed to the high volatility of Bitcoin data.

Since the Bitcoin data doesn’t show enough sign of seasonality, we don’t need to fit SARIMA models to this time series.

Trend Analysis

Using trend analysis we can try to identify the general trend for Bitcoin prices from 2013 to 2021. Is the increasing trend for Bitcoin prices statistically significant? We use the logarithm of prices, and fit a linear regression model with ARMA errors\(^{[6]}\). We have previously found that ARMA(3,2) model may be a reasonable model for the prices since the AIC of an ARMA(3,2) model is significantly smaller than the AIC of an ARMA(0,0) model.

The model with ARMA errors can be written as:

\[Y_t=\beta_0+\beta_1D_t+\eta_t\]

where \(D_t\) is days and \(\eta_t\) is an ARMA(3,2) error.

First, we will need to fit the ARMA(3,2) model with no trends. Then we can fit a linear regression model with ARMA errors and compare the results.

m0 = arima(data_ts_2013_log, order=c(3, 0, 2))
m0

## 
## Call:
## arima(x = data_ts_2013_log, order = c(3, 0, 2))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2  intercept
##       1.1926  -0.8782  0.6854  -0.2025  0.4982     7.2810
## s.e.  0.0621   0.0817  0.0370   0.0750  0.0440     2.9606
## 
## sigma^2 estimated as 0.003263:  log likelihood = 4276.23,  aic = -8538.46

days = time(data_ts_2013_log)
m1 = arima(data_ts_2013_log, order=c(3, 0, 2), xreg=days)
m1

## 
## Call:
## arima(x = data_ts_2013_log, order = c(3, 0, 2), xreg = days)
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2   intercept    days
##       1.1942  -0.8799  0.6831  -0.2077  0.4968  -1406.5755  0.7009
## s.e.  0.0622   0.0817  0.0368   0.0751  0.0440    349.6405  0.1733
## 
## sigma^2 estimated as 0.003251:  log likelihood = 4282.84,  aic = -8549.69

The difference of log likelihoods can now be calculated to determine whether the trend is statistically significant.

cat("Difference of log likelihoods:",m1$loglik-m0$loglik)

## Difference of log likelihoods: 6.614837

Since 6.61 >> 1.92, the general increasing trend for Bitcoin prices is significant.

Volatility Analysis

To further investigate the Bitcoin price from 2013 and onward, we will attempt to model the timeseries as conditionally heteroskedastic, meaning with non-constant variance throughout time. This is easily motivated by the plot of the growth rate \(y_t = \nabla log(x_t)\) where \(x_t\) is price at time t, as the data presents different variance at different points in time.

This motivates the use of an auto-regressive conditionally heteroskedastic model, or ARCH (similarly GARCH for generalized-ARCH). GARCH models\(^{[9]}\) model the return as a function of the variance at time t and the volatility as a function of both past returns and volatilities. The simplest GARCH model, the GARCH(1,1), models the return and volatility as \[r_t=\sigma_t\epsilon_t \\ \sigma^2_t = \alpha_0+\alpha_1r^2_{t-1}+\beta_1\sigma^2_{t-1}\] where \(\epsilon_t \sim iid. \mathcal{N}(0,1)\). The first term of the GARCH(p,q) model accounts for the order of lags of squared returns in the \(\sigma_t\) equation, and the second term accounts for the order of lags of the volatility. The GARCH(p,q) model retains the return model \(r_t=\sigma_t\epsilon_t\) but allows for greater lags in the volatility model, extending to \[\sigma_t = \alpha_0 + \sum_{j=1}^p\alpha_jr^2_{t-j}+\sum^q_{j=1}\beta_j\sigma^2_{t-j}\]

Building on the ARIMA(3,1,2) model, a low order model we investigated with relatively low AIC, we include a GARCH term in the model to account for the autocorrelation observed in the residuals. We examine the GARCH(1,1), GARCH(1,2), GARCH(2,1), and GARCH(0,1) models by fitting them to the residuals of the ARIMA(3,1,2) fit. We use the R package “tseries” and its “garch” function\(^{[7]}\) to do this.

library(tseries)
data.arima <- arima(data_ts_2013_log,order=c(3,1,2))
data_ts_2013_logdiff.g11 <- tseries::garch(data.arima$residuals, order = c(1,1))
data_ts_2013_logdiff.g12 <- tseries::garch(data.arima$residuals,order=c(1,2),maxiter=3000)
data_ts_2013_logdiff.g21 <- tseries::garch(data.arima$residuals,order=c(2,1))
data_ts_2013_logdiff.g01 <- tseries::garch(data.arima$residuals,order=c(0,1))

We achieve convergence on all models except the GARCH(1,2). We now use the parameters estimated from the GARCH(1,1) model to predict the variance of the time series as it had the highest log likelihood among the four models.

sigPred <- predict(data_ts_2013_logdiff.g11,data_ts_2013_logdiff)

plot(data_ts_2013_logdiff,type='l')
lines(sigPred[,1]^2,col="red", lwd=2)

Above is a plot of our time series with the predicted variance superimposed in red. By applying the GARCH(1,1) model to the residual of the ARIMA(3,1,2) fit, we see that the predicted variance matches rather well with the periods of high volatility in the data. Investigating the ACF of the squared residuals, we see that the GARCH(1,1) fit’s residuals have no significant autocorrelation. We can be confident that the GARCH(1,1) model fit to the ARIMA(3,1,2) model’s residuals provide a satisfactory explanation for the correlation in the data, as what is left appears to be uncorrelated noise.

len <- length(data_ts_2013_logdiff.g11$residuals)
garch_res2 <- data_ts_2013_logdiff.g11$residuals[7:len]^2
garch_res2_acf <- acf(garch_res2,plot=FALSE)
plot(garch_res2_acf,main="ACF of GARCH(1,1) fit squared residuals")

arima_res2 <- data.arima$residuals^2
arima_res2_acf <- acf(arima_res2,plot=FALSE)
plot(arima_res2_acf,main="ACF of ARIMA(3,1,2) fit squared residuals")

Conclusion

In this project we analyzed the Bitcoin prices from 2013 to 2021. After using all the time series analysis above, we have drawn the following conclusions:

ARIMA(3,1,2) is a reasonable model for Bitcoin prices on a log scale, it can probably be used to predict the approximate range of future Bitcoin prices.
As can be seen from the volatility analysis, Bitcoin prices are very volatile and difficult to predict. GARCH(1,1) model is reasonable for this time series considering the volatility.
The price of Bitcoin doesn’t show signs of apparent seasonal or cyclic changes. However, the increasing trend for Bitcoin price is statistically significant.

Bitcoin price is affected by many factors, like government policies, market situation and other unpredictable events like COVID-19, so it is hard to predict. Our models are based on basic time series techniques, so they might not fully account for the variability of Bitcoin price. Our models can be improved by making use of other relevant data.