Forecasting Bitcoin Prices with using Univariate GARCH model (version 1)

Introduction

Bitcoin is a cryptocurrency, a form of electronic cash. It is considered as one of most valuable currency and its popularity has increased tremendously in recent years. Besides being very popular, it’s has high volatility nature. Former, many studies have been focused on its impact on the market and other variance of the prices. Later it draws attention to deal with its vulnerabilities, which is considered as the utmost economic issues. In this study I have aim to focus on the model selection using ARIMA and used GARCH models to deal with the volatility nature of the series.

Data Preprocessing

d <- read.csv("Bitcoin_Historical_Price.csv")
head(d)

##      Date  Close
## 1 27/4/13 134.21
## 2 28/4/13 144.54
## 3 29/4/13    139
## 4 30/4/13 116.99
## 5  1/5/13 105.21
## 6  2/5/13  97.75

#Checking for obvious errors
which(is.na(d))

## integer(0)

d$Close <- as.numeric(as.character(gsub(",","",d$Close)))

#Converting the data set into time series object
tsObj<- ts(as.vector(d$Close), start=2013,frequency = 365)
plot.ts(tsObj)
title("Time Series plot of Bitcoin Prices ", sub = "(2013-2018)",
      cex.main = 1.5,   font.main= 4, col.main= "blue",
      cex.sub = 0.75, font.sub = 3, col.sub = "red")

Initial Analysis from time series plot:

It is crucial to derive the characteristics of time series to further deal the time dependecy. The major five characteristics of Bitcoin series as mentioned below:

Trend: From the above time series plots, we can observe different trends at different time intervals
Change in variance: Change in variance is obvious from the plot.
Seasonality: We cannot judge the seasonality of the series at this point but at glance look at series, we can assume that this series doesn’t have any seasonal patterns.
Intervention: We can see there is a steep increase in prices till the end of 2017. After there is a radical increase in prices which eventually increase from 5000 to 18000. This can also interpret as high volatile nature of the market at time in 2018.
Behavior: From time series plot one can depicts the autoregressive nature of the series but we can also expect a moving average especially from 2014 to 2017

ACF and PACF plot of bitcoin series:

The ACF plot clearly states the high correlation among successive points. It also shows a strong evidence of an existence of a trend as expected from the time series plot. The PACF plot shows one significant correlation on the plot. Unit root test proves the nature of non-stationarity of the series.

par(mfrow=c(1,2))
acf(tsObj, main=" ACF ")
pacf(tsObj, main=" PACF ")

ar(diff(tsObj))

## 
## Call:
## ar(x = diff(tsObj))
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.0543  -0.0472   0.0040  -0.0265   0.1460  -0.0157  -0.0185   0.0894  
##       9       10       11       12       13       14       15       16  
##  0.0520   0.1259   0.0110  -0.0603  -0.1135  -0.0991  -0.0260  -0.1369  
##      17       18       19       20       21       22       23       24  
##  0.0834  -0.0005   0.1371   0.1285   0.0286   0.0624  -0.0329  -0.0438  
##      25       26       27       28       29       30       31       32  
##  0.0205   0.0320  -0.1346  -0.0425  -0.0571  -0.0286  -0.0425  -0.0183  
##      33  
##  0.1225  
## 
## Order selected 33  sigma^2 estimated as  45811

adfTest(tsObj, lags=33, title=NULL,description=NULL)

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 33
##   STATISTIC:
##     Dickey-Fuller: -1.2548
##   P VALUE:
##     0.2164 
## 
## Description:
##  Mon Aug 05 19:54:58 2019 by user: dmnai

Ensuring the stationarity of the series:

Transformation:

Natural Logarithmic transformation is one of the best approach to look for stationarity of the sereis.

#Log Transformation
data.trans.log<-log(tsObj)

ar(data.trans.log)

## 
## Call:
## ar(x = data.trans.log)
## 
## Coefficients:
##      1  
## 0.9988  
## 
## Order selected 1  sigma^2 estimated as  0.004576

adfTest(data.trans.log, lags =1, title =NULL,description =NULL)

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: 1.3998
##   P VALUE:
##     0.9588 
## 
## Description:
##  Mon Aug 05 19:54:58 2019 by user: dmnai

The ACF and PACF plot are similar to the original series and we can still suspect the non-stationarity of the series and this is also supported by Unit root test.

par(mfrow=c(1,2))
acf(data.trans.log, ci.type='ma', main=" ACF of transformed data")

pacf(data.trans.log,  main="PACF of transformed data")

Differencing as my next approach to achieve the stationarity of the series.

#Differencing
p<-diff(data.trans.log) 
ar(p)

## 
## Call:
## ar(x = p)
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
## -0.0012  -0.0140   0.0069   0.0226   0.0383   0.0607  -0.0215  -0.0010  
##       9       10       11       12       13       14       15       16  
## -0.0042   0.0495   0.0563  -0.0066   0.0058   0.0052   0.0024  -0.0111  
##      17       18       19       20       21       22       23       24  
##  0.0688   0.0064  -0.0117   0.0552  -0.0230   0.0292  -0.0450  -0.0203  
##      25       26       27       28       29       30       31  
## -0.0139   0.0343   0.0198  -0.0436  -0.0376  -0.0472   0.0386  
## 
## Order selected 31  sigma^2 estimated as  0.001853

At the first difference of the transformed series, we can observe the plots of ACF and PACF shows a bit difference to the previous steep decreasing pattern. As we can ensure the assumption of stationarity with Unit-Root test.

par(mfrow=c(1,2))
acf(p, ci.type='ma', main="ACF of 1st differnce")
pacf(p, main="PACF of 1st differnce")

With reference to the Dickey-Fuller Test, p-value is less than the 0.05 and we can reject the null hypothesis stating the non-stationarity. Hence , we can proceed further for model selection .

adfTest(p, lags =31, title =NULL,description =NULL) # Significance achieved  but acf and Pacf has lot of significant lags

## Warning in adfTest(p, lags = 31, title = NULL, description = NULL): p-value
## smaller than printed p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 31
##   STATISTIC:
##     Dickey-Fuller: -7.7432
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Aug 05 19:54:58 2019 by user: dmnai

Model Selection:

We usually select the order of ARIMA from ACF and PACF. ACF and PACF plot of stationary series shows many significant lags. In order to a clear picture of the p and q, EACF plot is an alternative and the best possible candidate models can be ARIMA (1,1,1) and ARIMA (0,1,1).

#EACF for candidate Models
eacf(p) # ARIMA (1,1,1) ARIMA(0,1,1)

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

Parameter Estimation:

As part of model selection criteria, the candidate model parameters has been tested with help of coeftest() with different methods like conditional sum of squares and maximum likelihood.

The below code chuck results shows that AR(1), MA(1) are highly significant for both candidate models i.e.., ARIMA(1,1,1) and ARIMA(0,1,1).

#ARIMA (1,1,1) and pararmeter estimation by minimizing conditional sum-of-squares 
model_111_css = arima(tsObj,order=c(1,1,1),method='CSS')
coeftest(model_111_css)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.725127   0.081640 -8.8821 < 2.2e-16 ***
## ma1  0.801183   0.070923 11.2966 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA (1,1,1) and pararmeter estimation by Maximum Likelihood
model_111_ml = arima(tsObj,order=c(1,1,1),method='ML')
coeftest(model_111_ml)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.717682   0.082668 -8.6815 < 2.2e-16 ***
## ma1  0.794653   0.072230 11.0018 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA (0,1,1) and pararmeter estimation by Maximum Likelihood
model_011_ml = arima(tsObj,order=c(0,1,1),method='ML')
coeftest(model_011_ml)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 0.080522   0.022266  3.6163 0.0002988 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA (1,1,1) has been selected as the best possible model for the series which corresponds to lowest AIC score.

sc.AIC = AIC(model_011_ml,model_111_ml)
sort.score <- function(x, score = c("bic", "aic")){
  if (score == "aic"){
    x[with(x, order(AIC)),]
  } else if (score == "bic") {
    x[with(x, order(BIC)),]
  } else {
    warning('score = "x" only accepts valid arguments ("aic","bic")')
  }
}

sort.score(sc.AIC, score ="aic")

##              df      AIC
## model_111_ml  3 29217.15
## model_011_ml  2 29227.94

Overfitting:

Overfitting is a good practice to test the model selected from candidate models. As such the selected model is ARIMA (1,1,1). Models ARIMA(1,1,2) ,ARIMA(2,1,1) are used to test overfitting and they produced NA’s which eventually represents the selected model is good enough to fit the series.

#Overfitting ARIMA(1,1,1) the ARIMA(2,1,1) and ARIMA(1,1,2).

#ARIMA(1,1,2)
#Maximum likelihood
model_112_ml = arima(tsObj,order=c(1,1,2),method='ML')
coeftest(model_112_ml)

## Warning in sqrt(diag(se)): NaNs produced

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.038775         NA      NA       NA
## ma1  0.039536         NA      NA       NA
## ma2 -0.037663         NA      NA       NA

#CSS
model_112_CSS = arima(tsObj,order=c(1,1,2),method='CSS')
coeftest(model_112_CSS)

## Warning in sqrt(diag(se)): NaNs produced

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.038687         NA      NA       NA
## ma1  0.039662         NA      NA       NA
## ma2 -0.037697         NA      NA       NA

#ARIMA(2,1,1)
model_211_ml = arima(tsObj,order=c(2,1,1),method='ML')
coeftest(model_211_ml)

## Warning in sqrt(diag(se)): NaNs produced

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.039405         NA      NA       NA
## ar2 -0.032123         NA      NA       NA
## ma1  0.038631         NA      NA       NA

#CSS
model_211_CSS = arima(tsObj,order=c(2,1,1),method='CSS')
coeftest(model_211_CSS)

## Warning in sqrt(diag(se)): NaNs produced

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.039319         NA      NA       NA
## ar2 -0.032081         NA      NA       NA
## ma1  0.038889         NA      NA       NA

Residual Analysis:

Residuals which is more we care about during model fitting. As of they tells about how extend we had a good fit for the model. The below code chunks tell about the residuals in terms of:

Shapiro-Wilk normality test: This hypothesis states that residuals are highly correlated and are not normally distributed.

Ljung-Box test: This test visualizes the p values of the residuals which are all below the significance level

Time series plot of Standardized Residuals and QQplot of standardized suspects the volatility of the series.

residual.analysis <- function(model, std = TRUE,start = 2, class = c("ARIMA","GARCH","ARMA-GARCH")[1]){
  # If you have an output from arima() function use class = "ARIMA"
  # If you have an output from garch() function use class = "GARCH"
  # If you have an output from ugarchfit() function use class = "ARMA-GARCH"
  library(TSA)
  library(FitAR)
  if (class == "ARIMA"){
    if (std == TRUE){
      res.model = rstandard(model)
    }else{
      res.model = residuals(model)
    }
  }else if (class == "GARCH"){
    res.model = model$residuals[start:model$n.used]
  }else if (class == "ARMA-GARCH"){
      res.model = model@fit$residuals
  }else {
    stop("The argument 'class' must be either 'ARIMA' or 'GARCH' ")
  }
  par(mfrow=c(3,2))
  plot(res.model,type='o',ylab='Standardised residuals', main="Time series plot of standardised residuals")
  abline(h=0)
  hist(res.model,main="Histogram of standardised residuals")
  acf(res.model,main="ACF of standardised residuals")
  pacf(res.model,main="PACF of standardised residuals")
  qqnorm(res.model,main="QQ plot of standardised residuals")
  qqline(res.model, col = 2)
  print(shapiro.test(res.model))
  k=0
  LBQPlot(res.model, lag.max = 30, StartLag = k + 1, k = 0, SquaredQ = FALSE)
}

residual.analysis(model_111_ml, std = TRUE,start = 1)

## Loading required package: lattice

## Loading required package: leaps

## Loading required package: ltsa

## Loading required package: bestglm

## 
## Attaching package: 'FitAR'

## The following object is masked from 'package:forecast':
## 
##     BoxCox

## 
##  Shapiro-Wilk normality test
## 
## data:  res.model
## W = 0.49007, p-value < 2.2e-16

## Warning in (ra^2)/(n - (1:lag.max)): longer object length is not a multiple
## of shorter object length

Dealing with Conditional Heteroscedaticity:

To deal with volatile nature of bitcoin series, I have transformed the series as returns of the bitcoin.

#Volatility 
r.bitcoin = diff(log(tsObj))*100
plot(r.bitcoin)
title("Plot of returns of Bitcoin", sub = "(2013-2018)",
      cex.main = 1.5,   font.main= 4, col.main= "blue",
      cex.sub = 0.75, font.sub = 3, col.sub = "red")

McLeod Li test confirms the volatile nature as almost at all lags the p-values fall below the significance levels.

McLeod.Li.test(y=r.bitcoin,main="McLeod-Li test statistics for Daily return series")

In order to get an order of GARCH , we further transform the return series into absolute values and squared return values.

abs = abs(r.bitcoin)
sqr = r.bitcoin^2

GARCH Model specification:

From ACF and PACF we see many lags are significant. Hence, we plot EACF to get the candidate models.

Obtained Models from EACF are ARMA (1,2), ARMA (2,2), ARMA (1,3). As GARCH order is GARCH (max (p, q), p) the candidate models are as following GARCH (2,1), GARCH (3,1), GRACH (2,2)

par(mfrow=c(1,2))
acf(abs, ci.type="ma",main=" ACF for abs. returns")
pacf(abs, main=" PACF plot for abs.returns")

eacf(abs) # ARMA (1,2),ARMA (2,2),ARMA (1,3) # GARCH(2,1),GARCH(3,1), GRACH (2,2)

## 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 x o o o x o  o  o  o 
## 2 x x o o o x o o o x 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 x o o o o o  o  o  o 
## 6 x x x x x o o o o o o  o  o  o 
## 7 x x x x x o o o o o o  o  o  o

From the squared returns ACF and PACF plot, it is not that clear to derive the order of p and q. Hence, I approach EACF and the order of ARMA are ARMA (2,3), ARMA (3,3), ARMA (2,4). Thus, GARCH candidate models would be GARCH (3,2) GARCH (3,3) GARCH (4,2)

par(mfrow=c(1,2))
acf(sqr, ci.type="ma",main="ACF  for sqr. return")
pacf(sqr, main="PACF for sqr. return")

eacf(sqr) # ARMA(2,3) ARMA(3,3) ARMA(2,4) # GARCH(3,2) GARCH(3,3) GARCH(4,2)

## 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 x o x o  o  x  o 
## 2 x x x o o x o o o x o  o  x  o 
## 3 x x x o o x o o o x o  o  x  o 
## 4 x x x o o o o o o o o  o  x  o 
## 5 x x x o 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 x o o o o  o  o  o

Finally, the total candidate models from above are GARCH (2,1), GARCH (3,1), GRACH (2,2), GARCH (3,2), GARCH (3,3), GARCH (4,2)

MODEL ESTIMATION:

GARCH (2,1):

All model parameters are highly significant and Box- Ljung test confirms that the squared residuals are normally distributed.

# GARCH(2,1)
m.21 = garch(r.bitcoin,order=c(2,1),trace =FALSE)
summary(m.21)

## 
## Call:
## garch(x = r.bitcoin, order = c(2, 1), trace = FALSE)
## 
## Model:
## GARCH(2,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -8.39977 -0.36387  0.05552  0.50986  4.73252 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0   0.56551     0.05010   11.287  < 2e-16 ***
## a1   0.17891     0.01157   15.461  < 2e-16 ***
## b1   0.22436     0.04621    4.855  1.2e-06 ***
## b2   0.57941     0.04145   13.978  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 5021.7, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.061595, df = 1, p-value = 0.804

GARCH (2,2):

This model can be interpreted as an overfit model of GARCH(2,1) and p values from residual tests confirms that residuals are highly correlated. Thus this model is not consider to be a good fit.

# GRACH (2,2)
m.22 = garch(r.bitcoin,order=c(2,2),trace =FALSE)
summary(m.22)

## 
## Call:
## garch(x = r.bitcoin, order = c(2, 2), trace = FALSE)
## 
## Model:
## GARCH(2,2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.41911 -0.29952  0.04549  0.42410  5.88351 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 1.515e+01   1.137e+00   13.325  < 2e-16 ***
## a1 1.504e-01   1.391e-02   10.814  < 2e-16 ***
## a2 1.091e-01   2.413e-02    4.522 6.12e-06 ***
## b1 6.463e-14   1.059e-01    0.000    1.000    
## b2 9.075e-03   5.462e-02    0.166    0.868    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 2964.9, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 6.4863, df = 1, p-value = 0.01087

GARCH (3,1):

This model can be interpreted as an overfit model of GARCH(2,1) and GARCH (2,2). This model may not be consider to be a good fit.

# GARCH(3,1)
m.31 = garch(r.bitcoin,order=c(3,1),trace =FALSE)
summary(m.31)

## 
## Call:
## garch(x = r.bitcoin, order = c(3, 1), trace = FALSE)
## 
## Model:
## GARCH(3,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.84647 -0.29858  0.04541  0.42561  5.57910 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 1.514e+01   6.883e-01   22.000   <2e-16 ***
## a1 2.010e-01   1.704e-02   11.796   <2e-16 ***
## b1 2.583e-02   1.401e-02    1.843   0.0653 .  
## b2 2.381e-02   1.231e-02    1.934   0.0531 .  
## b3 2.342e-13   2.707e-02    0.000   1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 3284.3, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 2.1099, df = 1, p-value = 0.1463

GARCH (3,2):

This model can be interpreted as an overfitting model and p values from residual tests confirms that residuals are highly correlated. Thus this model is not consider to be a good fit.

# GARCH(3,2) 
m.32 = garch(r.bitcoin,order=c(3,2),trace =FALSE)
summary(m.32)

## 
## Call:
## garch(x = r.bitcoin, order = c(3, 2), trace = FALSE)
## 
## Model:
## GARCH(3,2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.46412 -0.30530  0.04672  0.43204  5.89009 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 1.420e+01   1.308e+00   10.853  < 2e-16 ***
## a1 1.509e-01   1.391e-02   10.845  < 2e-16 ***
## a2 1.131e-01   2.733e-02    4.140 3.48e-05 ***
## b1 6.634e-14   1.191e-01    0.000    1.000    
## b2 6.971e-03   7.610e-02    0.092    0.927    
## b3 2.216e-02   3.665e-02    0.605    0.545    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 2965.6, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 5.7328, df = 1, p-value = 0.01665

GARCH (3,3):

This model can be interpreted as an overfitting model and p values from residual tests confirms that residuals are highly correlated. Thus, this model is not consider to be a good fit.

# GARCH(3,3)
m.33 = garch(r.bitcoin,order=c(3,3),trace =FALSE)
summary(m.33)

## 
## Call:
## garch(x = r.bitcoin, order = c(3, 3), trace = FALSE)
## 
## Model:
## GARCH(3,3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.54365 -0.31092  0.04697  0.43143  6.08994 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 1.325e+01   1.906e+00    6.953 3.57e-12 ***
## a1 1.423e-01   1.275e-02   11.166  < 2e-16 ***
## a2 9.502e-02   2.833e-02    3.354 0.000796 ***
## a3 7.758e-02   3.051e-02    2.543 0.010996 *  
## b1 2.111e-15   1.568e-01    0.000 1.000000    
## b2 1.281e-02   1.508e-01    0.085 0.932309    
## b3 1.904e-02   9.500e-02    0.200 0.841122    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 3101.7, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 7.111, df = 1, p-value = 0.007661

GARCH (4,2):

This model can be interpreted as an overfitting model and p values from residual tests confirms that residuals are highly correlated. Thus, this model is not considered to be a good fit.

# GARCH(4,2)
m.42 = garch(r.bitcoin,order=c(4,2),trace =FALSE)
summary(m.42)

## 
## Call:
## garch(x = r.bitcoin, order = c(4, 2), trace = FALSE)
## 
## Model:
## GARCH(4,2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.55804 -0.30917  0.04787  0.42873  5.81905 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 1.325e+01   1.559e+00    8.500  < 2e-16 ***
## a1 1.561e-01   1.444e-02   10.810  < 2e-16 ***
## a2 1.163e-01   3.916e-02    2.969  0.00299 ** 
## b1 4.709e-14   2.058e-01    0.000  1.00000    
## b2 7.654e-03   1.478e-01    0.052  0.95870    
## b3 2.243e-02   1.013e-01    0.221  0.82476    
## b4 2.664e-02   6.069e-02    0.439  0.66062    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 2906.6, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 4.9389, df = 1, p-value = 0.02626

Model Selection:

Best possible model is selected by AIC scores of the models. From the below sort function, GARCH(2,1) would be the best model for the return series.

sc.AIC = AIC(m.21,m.32,m.33,m.22,m.42,m.31)
sort.score(sc.AIC, score ="aic")

##      df      AIC
## m.21  4 11599.99
## m.33  7 11907.11
## m.42  7 11916.49
## m.32  6 11945.45
## m.22  5 11982.50
## m.31  5 12022.18

Model Fitting:

Model selected :- ARMA of order(1,1) and GARCH(2,1)

We can observe three lags namely at 6 ,11 and 26 are above significance level from ACF of Standardized residuals but they are highly significant which barely contributes to 8% significance and ACF of squared residuals look good enough.

If we notice the QQ plot normality is assumed very bad hence I have tried with other distribution following like t-distribution and skewed t-distribution to get a view of performance and residuals along with Box-Ljung residual test.

#ARMA(1,1)+GARCH(2,1)
model21<-ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(2, 1)), 
                  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), 
                  distribution.model = "norm")
                  
m.21<-ugarchfit(spec=model21,data=r.bitcoin, out.sample = 100)
plot(m.21,which="all")

## 
## please wait...calculating quantiles...

Model Diagnostics

From the below tests , p values are significant which contradicts the qqplot but the goodness fit is very poor. Thus, normality is not satisfied. Hence, I tried with t-distribution and skew t-distribution as shown below.

The goodness of fit increases for t-distribution but the Box-Lung test p-values fall critically low. Hence, we can consider the model with assumption of normal distribution is good enough for the forecasting.

#normality assumption is not satisfied
m.21

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(2,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.096575    0.075659  1.276442 0.201799
## ar1     0.641911    0.591684  1.084887 0.277972
## ma1    -0.607595    0.612502 -0.991989 0.321203
## omega   0.377305    0.113078  3.336669 0.000848
## alpha1  0.122307    0.020557  5.949545 0.000000
## alpha2  0.000000    0.032767  0.000002 0.999998
## beta1   0.867308    0.025108 34.542494 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.096575    0.072574  1.330697 0.183289
## ar1     0.641911    1.045880  0.613751 0.539380
## ma1    -0.607595    1.082943 -0.561059 0.574757
## omega   0.377305    0.409747  0.920824 0.357142
## alpha1  0.122307    0.031353  3.900943 0.000096
## alpha2  0.000000    0.090762  0.000001 0.999999
## beta1   0.867308    0.091526  9.476093 0.000000
## 
## LogLikelihood : -5526.648 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       5.4546
## Bayes        5.4739
## Shibata      5.4545
## Hannan-Quinn 5.4617
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic  p-value
## Lag[1]                      2.665 0.102569
## Lag[2*(p+q)+(p+q)-1][5]     3.197 0.353723
## Lag[4*(p+q)+(p+q)-1][9]    10.705 0.004549
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.570  0.2102
## Lag[2*(p+q)+(p+q)-1][8]      2.874  0.7086
## Lag[4*(p+q)+(p+q)-1][14]     4.959  0.7745
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]     1.472 0.500 2.000  0.2251
## ARCH Lag[6]     1.534 1.461 1.711  0.6017
## ARCH Lag[8]     1.797 2.368 1.583  0.7823
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.5828
## Individual Statistics:              
## mu     0.32516
## ar1    0.04374
## ma1    0.04715
## omega  0.15776
## alpha1 0.09109
## alpha2 0.08410
## beta1  0.12370
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value    prob sig
## Sign Bias          1.80271 0.07158   *
## Negative Sign Bias 0.07369 0.94126    
## Positive Sign Bias 0.75321 0.45141    
## Joint Effect       4.26599 0.23414    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     329.3    1.912e-58
## 2    30     364.0    1.374e-59
## 3    40     372.0    6.434e-56
## 4    50     393.5    2.596e-55
## 
## 
## Elapsed time : 0.342247

t-distribution

A better qqplot but the auto correlation in standardised residuals are significant and Ljung test does support this assumption.

Though this model has better goodness of fit, we can’t consider due to above stated reasons.

model21_t_dis<-ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(2, 1)), 
                  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), 
                  distribution.model = "std")
                  
m.21_t_dis<-ugarchfit(spec=model21_t_dis,data=r.bitcoin, out.sample = 100)
m.21_t_dis

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(2,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu       0.15054    0.042536  3.539002 0.000402
## ar1      0.36245    0.323954  1.118829 0.263213
## ma1     -0.38919    0.320011 -1.216180 0.223916
## omega    0.21465    0.097928  2.191925 0.028385
## alpha1   0.15499    0.032900  4.710994 0.000002
## alpha2   0.00000    0.041697  0.000002 0.999999
## beta1    0.84401    0.027157 31.078554 0.000000
## shape    3.26089    0.178136 18.305605 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu       0.15054    0.045705  3.293606 0.000989
## ar1      0.36245    0.169267  2.141291 0.032251
## ma1     -0.38919    0.171458 -2.269889 0.023214
## omega    0.21465    0.214566  1.000395 0.317120
## alpha1   0.15499    0.036522  4.243793 0.000022
## alpha2   0.00000    0.058694  0.000001 0.999999
## beta1    0.84401    0.052385 16.111582 0.000000
## shape    3.26089    0.200703 16.247313 0.000000
## 
## LogLikelihood : -5241.388 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       5.1744
## Bayes        5.1965
## Shibata      5.1743
## Hannan-Quinn 5.1825
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      17.77 2.498e-05
## Lag[2*(p+q)+(p+q)-1][5]     23.02 0.000e+00
## Lag[4*(p+q)+(p+q)-1][9]     32.47 2.776e-15
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                     0.01425  0.9050
## Lag[2*(p+q)+(p+q)-1][8]    2.39644  0.7939
## Lag[4*(p+q)+(p+q)-1][14]   5.02529  0.7664
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]     1.857 0.500 2.000  0.1730
## ARCH Lag[6]     2.067 1.461 1.711  0.4762
## ARCH Lag[8]     2.515 2.368 1.583  0.6376
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  10.0727
## Individual Statistics:             
## mu     0.5769
## ar1    0.1557
## ma1    0.1613
## omega  0.3921
## alpha1 1.3713
## alpha2 0.8089
## beta1  0.6284
## shape  0.5260
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           1.0088 0.3132    
## Negative Sign Bias  0.6838 0.4942    
## Positive Sign Bias  0.4735 0.6359    
## Joint Effect        2.1322 0.5454    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     45.08    6.688e-04
## 2    30     66.58    8.776e-05
## 3    40     69.20    2.049e-03
## 4    50     84.87    1.119e-03
## 
## 
## Elapsed time : 0.5186179

# skew t-distribution
#much similar to t-distribution qqplot and the auto correlation in standardised residuals are significant and Ljung test does support this assumption.
#Though this model has better goodness of fit than normality assumption, we donot consider due to above stated reasons
model21_t_skw<-ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(2, 1)), 
                  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), 
                  distribution.model = "sstd")
                  
m.21_t_skw<-ugarchfit(spec=model21_t_skw,data=r.bitcoin, out.sample = 100)
m.21_t_skw

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(2,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.073871    0.051159  1.443956 0.148751
## ar1     0.370369    0.305095  1.213945 0.224769
## ma1    -0.399277    0.300891 -1.326983 0.184514
## omega   0.219032    0.095571  2.291827 0.021916
## alpha1  0.157419    0.032875  4.788444 0.000002
## alpha2  0.000000    0.041246  0.000002 0.999999
## beta1   0.841581    0.026498 31.760628 0.000000
## skew    0.935237    0.023789 39.314645 0.000000
## shape   3.285786    0.180643 18.189410 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.073871    0.060214  1.226802 0.219897
## ar1     0.370369    0.164032  2.257908 0.023951
## ma1    -0.399277    0.166132 -2.403368 0.016245
## omega   0.219032    0.203083  1.078532 0.280796
## alpha1  0.157419    0.036315  4.334847 0.000015
## alpha2  0.000000    0.056466  0.000001 0.999999
## beta1   0.841581    0.049357 17.051050 0.000000
## skew    0.935237    0.025366 36.870083 0.000000
## shape   3.285786    0.202799 16.202184 0.000000
## 
## LogLikelihood : -5237.891 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       5.1719
## Bayes        5.1968
## Shibata      5.1719
## Hannan-Quinn 5.1810
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      18.27 1.918e-05
## Lag[2*(p+q)+(p+q)-1][5]     23.71 0.000e+00
## Lag[4*(p+q)+(p+q)-1][9]     33.14 9.992e-16
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                    0.005683  0.9399
## Lag[2*(p+q)+(p+q)-1][8]   2.381750  0.7964
## Lag[4*(p+q)+(p+q)-1][14]  5.014373  0.7677
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]     1.809 0.500 2.000  0.1786
## ARCH Lag[6]     2.036 1.461 1.711  0.4829
## ARCH Lag[8]     2.456 2.368 1.583  0.6494
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  10.7665
## Individual Statistics:             
## mu     0.6091
## ar1    0.1609
## ma1    0.1680
## omega  0.3403
## alpha1 1.4412
## alpha2 0.8586
## beta1  0.6368
## skew   0.1962
## shape  0.5284
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.1 2.32 2.82
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           1.3827 0.1669    
## Negative Sign Bias  0.8844 0.3766    
## Positive Sign Bias  0.3916 0.6954    
## Joint Effect        3.1798 0.3647    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     52.86    4.909e-05
## 2    30     68.83    4.381e-05
## 3    40     77.16    2.609e-04
## 4    50     81.82    2.266e-03
## 
## 
## Elapsed time : 0.99175

Forecasting:

The below code chunk produces the forecast of return series.

forc = ugarchforecast(m.21, data = bitcoin, n.ahead = 10, n.roll =10)
print(forc)

## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 10
## Roll Steps: 10
## Out of Sample: 10
## 
## 0-roll forecast [T0=1975-07-23 10:00:00]:
##         Series Sigma
## T+1  -0.124287 3.657
## T+2  -0.045199 3.690
## T+3   0.005569 3.721
## T+4   0.038157 3.753
## T+5   0.059076 3.783
## T+6   0.072504 3.813
## T+7   0.081123 3.843
## T+8   0.086656 3.872
## T+9   0.090208 3.901
## T+10  0.092488 3.929

plot(forc, which= "all")

## Forecast of bitcoin series from returns:

Returns(t)= diff(log(price(t)))100 Returns(t+1)/100= log(price(t+1) - log((price(t+1-1) Returns(t+1)/100= log(price(t+1)/price(t)) Price(t+1) = price(t) e^(Returns(t+1)/100)

# Forecasting bitcoin price from returns forecast
  p.t_1 = 3810.43
  R_t <-c(-0.124287,-0.045199,0.005569,0.038157,
          0.059076,0.072504,0.081123,0.086656,0.090208,0.092488)
  p_t= 0 
  for (i in 1:10){
    p_t =  p.t_1 *((2.71828)^(R_t[i]/100))
    print(p_t)
    p.t_1=p_t
  }

## [1] 3805.697
## [1] 3803.977
## [1] 3804.189
## [1] 3805.641
## [1] 3807.89
## [1] 3810.652
## [1] 3813.744
## [1] 3817.051
## [1] 3820.495
## [1] 3824.031

Creating a data frame for fitted values and observed (real)values of bitcoin in order to calculate MASE.

#Creating the fitted and observed values for bitcoin

#importing the test set 
library(readxl)
Bitcoin_Prices_Forecasts <- read_excel("Bitcoin_Prices_Forecasts.xlsx")
test<-Bitcoin_Prices_Forecasts

#Creating the a column to the forcasted values
fitted<-c(3805.697,3803.977,3804.189,3805.641,3807.89,
               3810.652,3813.744,3817.051,3820.495,3824.031)

#Creating  the date column for the purporse of fitted values
Date<-c(test$Date)
test$Date <- as.Date(test$Date)

#Creating a new data frame for the fitted values
fitted_data <- data.frame(Date,fitted)
str(fitted)

##  num [1:10] 3806 3804 3804 3806 3808 ...

# Converting both fitted and observed data into 'ts' objects
fitted_data<- ts(as.vector(fitted_data$fitted), start=2019,frequency = 365)
observed<-ts(as.vector(test$`Closing price`), start=2019,frequency = 365)

MASE = function(observed , fitted ){
  # observed: Observed series on the forecast period
  # fitted: Forecast values by your model
  Y.t = observed
  n = length(fitted)
  e.t = Y.t - fitted
  sum = 0 
  for (i in 2:n){
    sum = sum + abs(Y.t[i] - Y.t[i-1] )
  }
  q.t = e.t / (sum/(n-1))
  MASE = data.frame( MASE = mean(abs(q.t)))
  return(list(MASE = MASE))
}

#MASE of test data and fitted data 
MASE(observed,fitted_data)

## $MASE
##       MASE
## 1 1.778466

Conclusion:

In this study, a detailed investigation of bitcoin prediction has been done using ARIMA and GARCH modeling techniques. Towards this end, first dealt with stationarity of the series which is followed by specifying ARIMA model. ARIMA (1,1,1) has been selected within the candidate models by careful examining the parameter selection and overfitting techniques. As, the residuals are highly correlated, GARcH modelling has been used to deal with changing variance. GARCH (2,1) has been selected as the best possible candidate model which is examined with Box-Ljung test.

As financial series have high changing variance, we didn’t get an exact normalized residual, but they are not significant above significance level of 10%. Thus,forecasted the values of the bitcoin prices for next 10 days and MASE has been calculated with observed values which is obtained as 1.78 approx.

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

Forecasting Bitcoin Prices with using Univariate GARCH model (version 1)

Manikanta Naishadu Devabhakthuni

Introduction

Data Preprocessing

Initial Analysis from time series plot:

ACF and PACF plot of bitcoin series:

Ensuring the stationarity of the series:

Transformation:

Model Selection:

Parameter Estimation:

Overfitting:

Residual Analysis:

Dealing with Conditional Heteroscedaticity:

GARCH Model specification:

MODEL ESTIMATION:

GARCH (2,1):

GARCH (2,2):

GARCH (3,1):

GARCH (3,2):

GARCH (3,3):

GARCH (4,2):

Model Selection:

Model Fitting:

Model Diagnostics

t-distribution

Forecasting:

Conclusion: