Statistical Analysis of Realized Volatility of Bitcoin Price using Heterogeneous Autoregressive and Generalized Autoregressive Conditional Heteroskedasticity Models

2.1. Bitcoin

Bitcoin is a peer-to-peer (P2P) electronic cash system introduced in the well-known paper of Nakamoto [12]. The P2P mechanism allows an ownership transfer from one party to another without a third-party intervention (financial institution). Payments can be made over the internet without any control or cost of a central authority for the 1^st time. Individuals who want to own bitcoins can either run a program on their own computer that implements the bitcoin protocol or create an account on a website that runs bitcoin for its users. The bitcoins are saved in a file called a wallet, which the user may secure and backup. These programs connect to each other over the internet forming P2P networks, making the system resistant to a central attack by M. Crosby et al. [13].

For now, bitcoins are generated through a process of mining. Any member operates as a miner using their computer knowledge to maintain the network. Mining is a computationally process that requires miners to find a solution to a mathematical problem to create a new block into the blockchain. Miners resolve this issue using the proof-of-work concept. This algorithm involves recurrently difficult mathematical problems until getting to a solution. The first miner to find a solution broadcast it to the network to verify it. Once verified, the block is added to the blockchain. Every 10 min on average is found a new answer and a bitcoin is created. The bitcoin protocol is designed to generate a new bitcoin gradually. The difficulty of solving problems is adjusted every 2 weeks at the rate of six blocks per hour. The size of the reward was initially 50 (genesis block) and it is halved every 4 years, this implies that the number of bitcoins in circulation will never exceed 21 million. Once the last bitcoin is generated, miners will instead be rewarded with transaction fees Lo [13].

2.2. GARCH Model

The GARCH model was first introduced by Bollerslev [7]. Back in these days, the concept of realized volatility modeling was not even introduced. At that period, the daily volatility was calculated as the squared daily return without taking into consideration any subintervals. The GARCH model is a conditional volatility model that allows the conditional variance to depend on the previous lags. It is based on the ARCH model by Engle [14], who used it to show that the conditional volatility is affected by volatility clustering. An autoregressive conditionally heteroskedastic (ARCH) model is a time series model with econometric applications that consider the variance of the current error term as a function of the variance of the error conditions of the previous time periods. One of the disadvantages of the ARCH model is that it responds slowly to large, unusual shocks. Thus, the need of an improvement of this model was crucial. Assuming an autoregressive moving average (ARMA) model for the error variance, then the model is a GARCH model. GARCH models were designed to deal with the problem of volatility clustering, which is the phenomenon where large changes in prices tend to cluster together, as A. J. M. Karim and N. M. Ahmed [15]; Botan et al. (2020) [16]; and R. F. Engle [14].

Before describing the GARCH model, the ARCH specification has to be introduced. The following return process has to be specified:

Where, μ_t is a drift term that is explained by the structural model and z_t is an independent shock with zero mean and unit variance, signifying that ε_t is normally distributed ε_t ~ Z(0, σ_t). The conditional variance in (1) can be transformed into a time-varying by specifying the ARCH (q) process:

Where c is a constant and a_i is the coefficient for the past squared shocks (ε²_t). Then the GARCH (p,q) model is derived by adding p lagged conditional variances, with orders p ≥ 1 and q ≥ 1:

Where β_j is the coefficients for the past conditional variances, p is the past squared error terms, and q is the past estimated volatility terms. When q = 0, then the above equation (3) reduces to an autoregressive conditional heteroskedastic (ARCH) model. Given a distribution of ε_t in equation (1) and setting p = q =1, then the GARCH (1, 1) is derived:

For which the condition c ≥ 0, α₁ ≥ 0, and β₁ ≥ 0 should stand for every positive value of σ_i. Since the GARCH model is non-linear, it cannot be estimated by an OLS regression like the HAR model. Thus, the Gaussian maximum likelihood (GMLE) method should be used for parameter estimation. When assuming normally distributed errors and starting from some parameter vector θ and a time series of size T (y_1,y₂… y_T), the GMLE method calculates the probability density for this specific sample by taking the product over all the marginal conditional probability densities of the observed data. In general, the GARCH model is using the returns to forecast volatility, and it depicts that today’s return consists of yesterday’s return plus some volatility part and this volatility is what we need. This model is also using a rolling regression method to forecast volatility, by moving one day ahead and leaving 1 day behind for every forecast, which means that the data window size remains stable.

2.3. Testing GARCH Effects (Test of Heteroscedasticity)

The availability of ARCH/GARCH effects may give serious model misspecification if they are ignored. Logically, ignoring ARCH effects will give the identification of ARMA models that are over-parameterized. In addition, as in heteroscedasticity, estimation assuming its absence will result in inappropriate standard errors of parameter estimates, which are typically smaller than what they should be. Therefore, it is important to check the presence of GARCH effects in time series modeling according to McLeod and Li [17]; Asraa et al. [18]; and Azhy et al. [19].

Two ways of testing GARCH effects are used. First is to check the Ljung-Box portmanteau Q statistics of α²_t.

McLeod and Li show that the sample autocorrelations of have α²_t asymptotic variance n⁻¹ and that portmanteau statistics calculated from their distribution is asymptotically Chi-square if α²_t are independent. Since the sample autocorrelations of a are also pertinent to the identification of a GARCH model for α²_t.

The second is the process of checking for conditional heteroscedasticity which is to utilize the Lagrange multiplier test of Engle. Think about the following model of regression for

Where y_t is the error term, m is a pre-determined positive integer, and n is the total number of data points in the series. Using the coefficient of determination from (5), Engle demonstrates that, under the null hypothesis H₀: α₁ = α₂ =… = α_m = 0, the variance of nR² is approximately distributed according to a chi-square distribution with m degrees of freedom by Weiss [20]; Nakamoto [12].

2.4. Identification of a GARCH Model

If the Ljung-Box statistics and the LaGrange multiplier (LM) test are significant, then conditional heteroscedasticity of α²_t is present, and we need to identify an appropriate GARCH model for α²_t. However, since the GARCH (1,1) model has been shown to be appropriate in many empirical studies, we may employ the GARCH (1,1) model at the beginning of the analysis. As the model is estimated, diagnostic checking procedures may be followed to see if the GARCH (1,1) model is okay, or if the orders of the GARCH model should be increased or decreased. Instead of using this trial-and-error approach, we may use the following procedure for the definition of a GARCH model for the {α²_t} series by Lon-Mu [21]; N. M. Ahmed and A. J. M. Karim [22], Dyhrberg [23].

2.5. Ljung-Box Q-Statistic

Adding to the visual inspection of the plotted autocorrelation, the Ljung-Box Q-Statistic is used for diagnostic checking by Box and Jenkins [24]. The Ljung-Box Q-Statistic is defined by equation (6)

Where n is the number of observations, K is the largest degree of freedom used, and r_j is the sample association function at the jth degree of freedom of a relevant time series a_t, for example. The statistical r_j for is then calculated as:

The Q-statistic was suggested for testing ARIMA and ARMA models; both the test statistics are determined by the calculation of the sample autocorrelation function for the residuals ε^{^}_t from those models. The similar test statistic based on different calculations using the autocorrelation function will be high benefit for small sample applicability, it is defined as Weiss [20]; M. S. Lo. [25].

where k is the number of lags considered in the test, and k is defined by:

k = K – m,

where K: Number of lags used in the test.

m: Number of parameters estimated in the mean and variance equations of the GARCH model, and r^*_j is:

2.6. Likelihood Function of GARCH Models

By defining α = [α₀, α₁,…, α_m, B₁,…, B_r,η]’, the log likelihood functions of α may be derived under the Normality assumption of ε_t. If ε_t is assumed to follow a normal distribution. However, practically, there is substantial evidence showing that this assumption may not all the time be satisfactory by Lon-Mu [21].

For the GARCH (1.1) model, the joint density of the observations α₁…α_T can be calculated as the product of the conditional densities, conditioning on the last observations from M. S. Lo. [25].

fα₁,….,α_T (α₁,…,α_T)

Easy to say, the marginal population will decrease as for the ARIMA (1.1) model. For k = 2.,T, the probability of α_k, given the values α₁…α_(k−1), is

Moreover, the conditional likelihood function given α_t and is σ²_t:

Where σ^*2_t=α₀+α₁α²_i-1+B₁ are obtained recursively. We substitute by σ²_t its expected value:

Using the logarithm and ignoring the constant term, we find that the log likelihood function is:

Where α=(α₁,…α_T)^- and σ²_t=(σ²₁,…σ²_T).

2.7. Model Checking of GARCH (r,m)

For a GARCH model, the standardized errors ε^{^}_t=α^{^}_t/σ^{^}_t are independent and identically distributed random mistakes that are associated with either a standard normal or a non-normal distribution, such as the standardized student-t distribution. As a result, one can assess the effectiveness of a fitted GARCH model by inspecting the series {ε^{^}_t}. Specifically, the sample autocorrelations and the Ljung-Box Q statistics of can be utilized to assess the effectiveness of the mean (primordial) equation and the sample autocorrelations and the Ljung-Box Q statistics of can be utilized to assess the validity of the volatility (secondary) equation A.J.M. Karim, and N.M. Ahmed [15]; Lon-Mu [21].

2.8. Forecasting the GARCH (1,1) Model

Forecasts of a GARCH model can be found using methods similar to those of an ARMA model. Consider the GARCH (1, 1) model in assume that the forecast origin is n. For a one-step-ahead forecast, we have:

Where α²_n and σ²_n are known at t = n, therefore, the one-step ahead forecast is:

For multi-step ahead forecasts, we use α²_t=σ²_tε²_t:

When t = n+1, the equation becomes:

Since E(ε (n+1)−1-|F_n) = 0, the two-step-ahead volatility forecast at the forecast origin n satisfies the equation:

In general, we have:

This outcome is identical to the result of an ARMA (1, 1) model with an AR polynomial of degree 1 − (α_1+B_1) B. By repeatedly changing the values in (20), the forward forecast can be written as:

Therefore:

Provided that α₁ + B₁ < 1

As a result, the multi-step-ahead predictions of volatility made by a GARCH (1,1) model match the unconditional variance of as the horizon for predictions increases to nothing if Var (a_t) is present Lon-Mu [21]; Cont [26]; and Nader et al. [27].

2.9. The HAR-RV Model

The idea of realized variance is based on these assumptions and Andersen et al. [28] provide us with an explanation in further detail on how efficient an estimator of volatility the realized volatility is, and moreover, how it outperforms traditional GARCH-type models. Equation (23) presents the model for estimating a daily Realized Variance (RV), where r is high-frequency intraday log-returns as described in equation (24).

where r is Log returns (or continuously compounded returns) are approximately equal to normal price returns, but hold significant benefits in simplicity in multi-period returns by Ruppert and Matteson [29]. Log returns are defined:

Based on the Heterogeneous Market Hypothesis (HMH), Corsi [8] proposes the HARRV as a model that will utilize three realized volatility components in an autoregressive manner, which all represent some time-dependent market component for the model. The following equations 25 and 26 consider the RV over the complementing horizons. They are quite simply the average of the daily RV, so for a weekly RV, we simply extend the model as following:

Moreover, the same definition for monthly volatility, but over 22 daily periods:

The added sum of these three volatilities can be regarded as an additive cascade of volatilities, each representing different components of market volatility. From this, it gets an almost long memory AR type of character (with lags one, five, and 22), but not strictly Corsi [8].

By expanding the expected values and utilizing straightforward recursive substitution, the volatility model will be given by a three-step cascade and has a form of something similar to three AR processes Corsi [8]):

Now given (27), all variables are directly observable and available in the data set. The parameters will be able to be estimated through a simple Ordinary Least Squares estimation (OLS) using the Newey-West covariance matrix estimator. However, due to possible serial correlation, a Newey-West (NW) covariance correction will be applied, since the effects of covariance and autocorrelation must be considered in the estimation.

3.1. Data Description

In this paper, daily observations are used of the bitcoin price, the sample period is June 31, 2017–January 31, 2022, which obtained from the Kaggle website [30], the researcher use R-language to obtain results. Fig. 1 below shows the time series plot of the series during the sample period. Since 2017, the bitcoin price has become more volatile. On October 13^th, 2017, bitcoin price breaks the $5,000 for the 1^st time, on November 28^th, 2017, the $10,000, and on December 18^th, 2017, hits all-time high just below $20,000. The bitcoin price time series can be observed in Fig. 1, non-linear trend and non-stationarity are the first geometrical properties that are shown in Fig. 1.

Fig. 1. Time series plots of the variable.

To build an appropriate model, the series that are used in analysis must be stationary; therefore, it should check the unit-root structure of the data. Although the above graph gives a rough idea about the stationarity structure of the series, we have applied the Augmented Dickey-Fuller test to the series to test unit roots. Table 1 exhibits the results from the ADF test applied to levels, first differences of the series.

TABLE 1: Unit root results of the lag variable

The ADF test results indicate that the variable is non-stationary by not rejecting the null hypothesis of unit-root at the level, but it is stationary after first differencing. Therefore, the researcher uses differenced series in its analysis. Fig. 2 below presents a time series plot of the differenced series.

Fig. 2. Time series plot of the differenced variable.

After achieving the stationarity condition of the series, we should fit the mean equation model as it is shown below:

From Table 2, P-value of the constant and lagged variable are less than the statistical significance (0.05), then one can say that the model is significant.

TABLE 2: The fit of the mean equation model

Fig. 3 shows that there is a prolonged period of high volatility from day 1 to the end of 2017, and also, there exists a prolonged period of low volatility from the end of 2018, to the beginning of 2019. This suggests that the residual or error term is conditionally heteroskedastic and it can be represented by the ARCH and GARCH models.

Fig. 3. The residuals of the mean equation.

In the series, hypothesis of ARCH has an effect or not on the mean equation. For this purpose, a heteroskedasticity test ARCH have been used to test the following hypothesis below, and its results are shown in Table 3:

TABLE 3: The heteroskedasticity test ARCH

Table 3 presents P-value of the heteroskedasticity ARCH test, which is less than the statistical level (0.05), then the hypothesis can be rejected; in another word, there exists an ARCH effect in the series. The researcher achieved two main assumptions of using the GARCH model, which are the stationary of the series and the effect of ARCH in the mean equation model, and then GARCH model can be used to forecast the volatility. The GARCH (1,1) model has been run, the results of its fit are shown in Table 4 below:

TABLE 4: The fit of the GARCH (1,1) model

From the above table, it is obvious that the estimators of the variance equation are significant, depending on P-value, which is less than 0.05. The residuals of the GARCH (1,1) model should be tested to find out that the model is suffered from serial correlation of residuals or not for the hypothesis below:

From Table 5, P-value of the Q-statistic test for the 36 laggs are greater than the statistical level (0.05), then we cannot reject the hypothesis, that is mean there is no serial correlation of residuals.

TABLE 5: Testing of ACF and PACF

Another test is the heteroskedasticity test to figure out that the postulated model is adequate or not for the following hypothesis below, and the results are shown in Table 6 below:

TABLE 6: The test of heteroskedasticity for the ARCH effect

Table 6 shows the test for heteroskedasticity and it is clear that P-value of the test is greater than the statistical level (0.05) then cannot be reject the hypothesis, which means that there is no heteroskedasticity.

The final test is a normality test to figure out that the residual of the postulated model is normally distributed or not form the following hypothesis below, and the results are shown in Table 7 below:

TABLE 7: A normality test for the residual of the model

Table 7 shows the test for residual normality and it is clear that P-value of the test is less than the statistical level (0.05) then one can be reject the hypothesis, which means that the residual are not normal and this is a good result because residual not normal, which indicate to good model. Moving now to the HAR-RV model, first, the heteroskedasticity test has to be done for the residuals. Using heteroskedasticity test, as shown in Table 3, that there is indeed heteroskedasticity was found. After conducting this test for the HAR-RV model, the realized volatility of yesterday (RV1) was significant at a level of 5%, the realized volatility of past week (RV5) was no again significant, and the realized volatility of past month (RV22) was no again significant at a level of 5%. All the coefficients were positive and the F-test for the model is very large, which is significant. The coefficient for RV1 is more than RV5 and RV22. Table 8 below provides an overview of the coefficients of the estimation.

TABLE 8: The forecasted values

Table 8 reports the results of the HAR model, where the in-sample forecasting results show consistently that the lagged RV has a strong and persistent positive relationship with future realized variance log (RVt: t+h), especially in the first forecasting horizons h = 1. To detect the best model between GARCH and HAR models, one can use the criteria below such shown in Table 9.

TABLE 9: Comparison of postulated models

Table 9 represents the comparison between models, according to all criteria; the HAR model is the best model than the GARCH model, which means that can be used HAR model to forecast the price of bitcoin cryptocurrency. In addition, the full-sample forecasting assumes the realized variance time series to be stable, so the researcher implements the rolling window method to allow the parameters to change over time, and then more reasonable comparisons can be obtained. The adaptive method mimics an investor who updates the forecasting model based on the most recent information. The window size T of the adaptive HAR models employed here is 90 days, that is, models are estimated using past 90-day samples, such as shown in Table 10. Moreover, the model is re-estimated every day. After the re-estimation of each day, the out-of-sample forecasts are performed in horizons h = 1; 5; 22, spontaneously. The parameters of the daily aggregated realized variances are evolving systematically, which justifies the adaptive forecasting method.

TABLE 10: he forecasted values

Table 10 reports about the forecast values. These forecasts are obtained by first estimating the parameters of the models on the full sample and then performing a series of static one-step-ahead forecasts. Fig. 4 reports the results for out-of-sample forecasts of the realized volatility in which the model is re-estimated daily.

Fig. 4. (a) Observed and forecasted realized volatility (RV), (b) forecasted RV only.

Fig. 4 illustrates that the upward and downward risk estimators serve complementary roles in forecasting over time. Specifically, the upward risk coefficients showed a gradual increase from 2017 until early 2021. After that point, they exhibited a stronger upward trend, indicating growing exposure to potential positive risks or gains in the forecasted variable.