Univariate and Multivariate GARCH Models Applied to Bitcoin Futures Option Pricing

: In this paper, the Heston–Nandi futures option pricing model is applied to Bitcoin futures options. The model prices are compared to market prices to give an indication of the pricing performance. In addition, a multivariate Bitcoin futures option pricing methodology based on a multivatiate GARCH model is developed. The empirical results show that a symmetric model is a better ﬁt when applied to Bitcoin futures returns, and also produces more accurate option prices compared to market prices for two out of three expiry dates considered. Author Contributions: Conceptualization, P.J.V. and E.M.; methodology, P.J.V.; software, P.J.V.; vali-dation, E.M.; formal analysis, P.J.V.; investigation, P.J.V.; data curation, P.J.V.; writing—original draft preparation, P.J.V.; writing—review and editing, E.M.; visualization, P.J.V.; supervision, E.M.; project administration, E.M.


Introduction
Cryptocurrencies, and especially Bitcoin, have received a lot of attention in the financial modelling literature in recent years. Financial modelling researchers and practitioners are faced with the problem of price discovery when derivatives on a new asset class are introduced. The focus of this paper is price discovery in the Bitcoin futures option market. Abraham (2020) explains that a valuation model for Bitcoin futures options can provide insight into a central-bank free currency. We consider vanilla options (univariate) and spread options (multivariate).
Modelling the historical returns of an asset as a univariate generalised autoregressive conditional heterskedasticty (GARCH) process is often used as a basis for price discovery in illiquid derivative markets. The model by Heston and Nandi (2000) is often used because it has a convenient closed-form solution. However, this is usually applied to spot price dynamics. This model was extended to futures options on commodities by Li (2019a). The Chicago Mercantile Exchange (CME) was the first established exchange to launch Bitcoin futures options in the first quarter of 2020 (Bharadwaj 2021). Therefore, Bitcoin futures options are actively traded, and model prices can be compared to market prices to give an indication of pricing performance.
An important factor to consider is the ability to model joint dynamics for the pricing of multivariate derivatives, when pricing derivatives on a new asset class. According to Alexander and Heck (2020), crypto-asset futures are exposed to significant basis risk. Therefore, spread options on Bitcoin futures are ideal for hedging basis risk. Spread options on Bitcoin futures do not actively trade. In this study, we consider a modelling approach for price discovery in the Bitcoin futures spread option market. The approach is based on work by Rombouts and Stentoft (2011), who derived the risk-neutral dynamics of the spot price processes for a general class of multivariate heteroskedasticity models. In this study, the model by Rombouts and Stentoft (2011) is extended to multivariate futures options.
The rest of this paper is structured as follows: Section 2 reviews the recent and relevant literature, Section 3 focuses on the theoretical framework (both univariate and mulitvariate options on Bitcoin futures), Section 4 presents the empirical results, and finally, Section 5 provides concluding remarks.

Literature Review
Research focusing on GARCH models applied to Bitcoin (and other crypto-assets) and Bitcoin derivative pricing is well documented in the literature. In a recent study, Fassas et al. (2020) made use of a vector error correction model to investigate the price discovery process in the Bitcoin market. Their empirical results indicate that volume traded in the futures market is more important than the volume traded in the decentralised spot market when incorporating new information about the value Bitcoin. In addition, Fassas et al. (2020) consider the volatility transmission between the Bitcoin spot and futures market by using multivariate GARCH models (BEKK and dynamic-conditional-correlation). There is evidence of cross-market effects when the variability of Bitcoin spot and futures returns are considered.
It is important to consider a reasonable forecast of Bitcoin returns when trading Bitcoin derivatives. In a study focusing on the use of information on the US-China trade war to forecast Bitcoin returns, Plakandaras et al. (2021) made use of ordinary least squares regression, least absolute shrinkage and selection operator techniques, and support vector regression. The authors also controlled for explanatory variables, which include: financial indices (including the volatility index), exchange rates, commodity prices, political uncertainty indices, and Bitcoin characteristics. Their empirical results indicate that Bitcoin returns are not affected by trade-related uncertainties.
In a recent paper, Shahzad et al. (2019) made use of the bivariate cross-quantilogram to determine whether Bitcoin exhibits safe haven properties (during extreme market conditions) for equity investments. The authors extend the work by Baur and Lucey (2010) to incorporate weak and strong safe haven assets. Furthermore, the safe haven properties of Bitcoin were also compared to those of gold and the general commodity index. Shahzad et al. (2019) conclude that Bitcoin, gold, and the general commodity index can be considered (at best) a weak safe haven asset in some cases. Fang et al. (2019) made use of the GARCH-MIDAS model to investigate how the long-run volatility of Bitcoin, global equities, bonds, and commodities evolves with global economic policy uncertainty. Their empirical results indicate that global economic policy uncertainty is significant for all variables, except bonds. Furthermore, Fang et al. (2019) also considered the impact of global economic policy uncertainty on the correlation between Bitcoin and global equities, commodities, and bonds. Based on the empirical results, the authors argue that Bitcoin can act as a hedge under specific economic uncertainty conditions.
In a study highlighting the role of Bitcoin futures, Chen and So (2020) focused on the relationship between Bitcoin spot and futures prices, with the focus on hedging. Chen and So (2020) tested the hedge performance of the naive method, ordinary least squares, and a dynamic hedge based on a bivariate BEKK-GJR-GARCH model. The results show that the hedge based on the bivariate GARCH model is the most reliable. Furthermore, Chen and So (2020) also show that the structure of Bitcoin volatility is significantly different after the introduction of Bitcoin futures.
In a recent study,  applied symmetric and asymmetric GARCH option pricing models to Bitcoin and CRIX (Cryptocurrency Index). The model Bitcoin option prices were compared to market option prices, which shows that the GARCH option pricing model produces reasonable price discovery. Furthermore, the implied volatility surfaces generated using symmetric and asymmetric GARCH models were compared. This comparison indicates that there is not a significant difference, implying that the symmetric model is a better choice as it is more efficient. Jalan et al. (2020) focused on the pricing and risk of Bitcoin options. Jalan et al. (2020) compared the option prices obtained from classical option pricing models, i.e., the Black-Scholes-Merton model and the Heston-Nandi model. Furthermore, Jalan et al. (2020) also compared the risk (the Greeks) of Bitcoin options to those of traditional commodity options. Jalan et al. (2020) conclude that the classical models produce prices that are slightly different compared to the market. Their empirical results also indicate that Bitcoin deltas are more stable over time compared to traditional commodities. This implies that investors in Bitcoin options are protected from undue Bitcoin price changes.
In another recent study, Siu and Elliot (2021) made use of the self-exciting threshold autoregressive model (to incorporate regime switching) in conjunction with GARCH (Heston-Nandi) to model Bitcoin return dynamics, for the pricing of Bitcoin options. According to Siu and Elliot (2021), conditional heteroskedasticity has a significant impact on Bitcoin option prices. However, the impact of self-exciting threshold autoregressive terms seems to be marginal.
Limited research has focused on the GARCH option pricing model applied to futures options. Li (2019a) extended the Heston-Nandi model to futures options. The overall purpose was the pricing of crude oil futures options. Li (2019a) concludes that optionimplied filtering is superior when compared to futures-based filtering when pricing crude oil futures options. Li (2019b) also applied the model to natural gas futures. However, this approach has not been applied to cryptocurrencies.
The application of GARCH models to multivariate option pricing is also well documented in the literature. Duan and Pliska (2004) developed an option valuation theory for cointegrated assets; the model was used for the pricing of spread options with equity underlying assets. In a recent study, Mahringer and Prokopczuk (2015) applied the model by Duan and Pliska (2004) to the pricing of crack spread options (futures returns were modelled). Mahringer and Prokopczuk (2015) compared this to univariate modelling of the crack spread. Their empirical results show that the univariate approach is superior for the pricing of crack spread options. Rombouts and Stentoft (2011) derived the risk-neutral dynamics (of the spot price processes) for a general class of multivariate heteroskedasticity models. In addition, a feasible way to price multivariate options is also provided. Rombouts and Stentoft (2011) applied the models to options on equity indices. Their empirical results indicate that correlation risk and non-Gaussian features are important factors to consider when pricing multivariate options. In this paper, the framework by Rombouts and Stentoft (2011) is extended to futures options and applied to Bitcoin futures spread options. The theoretical framework is considered in the next section.

Theoretical Framework
In this paper, we test the pricing performance of the Heston-Nandi futures option pricing model when applied to Bitcoin. Furthermore, we extend the work by Rombouts and Stentoft (2011) to multivariate futures options to price spread options on Bitcoin futures. The section is divided into three subsections. The first part focuses on the Heston-Nandi futures option pricing model, the second considers the multivariate GARCH option pricing framework, and finally, the focus of the third subsection is the constant conditional correlation (CCC) and dynamic conditional correlation (DCC) GARCH (multivariate) models.

Heston-Nandi Futures Option Pricing Model
The model by Heston and Nandi (2000) is widely used in the literature for the pricing of vanilla options. This model was extended to futures options by Li (2019a), who applied the model to crude oil futures. The futures dynamics under the real-world measure P are given by (Li 2019b): where F t,T is the futures price at time t with expiry T, λ is the unit risk premium, and z t is a standard normal random variable. The conditional variance takes the following form: In this study, the parameters λ, α 0 , α 1 , γ, and β 1 are calibrated to historical Bitcoin futures returns (under the real-world measure) using maximum-likelihood estimation. When estimating the symmetric Heston-Nandi (HN) model, the asymmetry parameter γ takes a value of zero. For the pricing of futures options, risk-neutral dynamics are required. According to Li (2019a), the risk-neutral futures price process in the HN framework is given by The risk-neutral conditional variance takes the following form, Given the risk-neutral dynamics, a closed-form formula for a European call option on a futures contract can be obtained, see Li (2019a) for more detail. The parameters are estimated using maximum likelihood estimation; the log-likelihood function is given by (Wang et al. 2017): where N is the estimation sample size. The multivariate GARCH futures option pricing model is outlined in the next section.

Multivariate GARCH Futures Option Pricing Model
We assume the following futures return dynamics ln F j,t,T j F j,t−1,T j = R j,T j under the real-world measure P: where F j,t,T j is the futures price of asset j with expiry T j at time t; µ j,T j is the conditional mean of asset j; f (·) denotes the cumulant generating function; c j is a vector of zeros except for position j, which takes a value of one. Furthermore, we assume a multivariate heteroskedastic process; therefore, where z t is identically and independently distributed with mean zero and covariance matrix equal to the identity matrix under the real-world measure P. In addition, H t is an n × n matrix of full rank, more specifically, where Σ t is the conditional covariance matrix, driven by a multivariate GARCH process. In order to obtain the risk-neutral dynamics (Q measure) required for pricing, we use the following Radon-Nikodym derivative, where F t is the information set available at time t, and ν i is an N dimensional vector sequence. Rombouts and Stentoft (2011) prove that Equation (3) is in fact a Radon-Nikodym derivative. Furthermore, it can also be shown (using the tower property) that as required.
Proposition 1. The risk-neutral measure Q defined by the Radon-Nikodym derivative in Equation (3) is an equivalent martingale if and only if for j = 1, . . . , n.
Proof. Clark (2014) explains that under the T j -forward measure, the futures price process is driftless. However, when assuming constant interest rates, the risk-neutral (Q) and T-forward measures are equivalent. Therefore, Hence, by making use of the Radon-Nikodym derivative, Using the driftless property of the futures price under the risk-neutral measure, it follows that which completes the proof.
To derive the risk-neutral dynamics, the following lemma from Rombouts and Stentoft (2011) is required: Lemma 1. Under the risk-neutral measure, the conditional moment generating function takes the following form, Proof. See Rombouts and Stentoft (2011).
Using the lemma above, the following expression for the risk-neutral cumulant generating function is obtained Hence, for any choice of ν t , the risk-neutral dynamics can be obtained by substituting Equations (4) and (5) into the mean Equation (2). The risk-neutral futures log-returns are given by: R * j,T j = f * (−c j ) + * j,T j , * denotes that the variables are considered under the risk-neutral measure. In this paper, we assume a multivariate Gaussian distribution. Rombouts and Stentoft (2011) show that the conditional cumulant-generating function is given by: where u is an arbitrary vector, if the multivariate Gaussion distribution is assumed. By substituting Equation (6) into Equation (4), it is easily shown that ν t takes the following form: In addition, the risk-neutral cumulant-generating function is given by: We assume the following mean model, where diagΣ t is a diagonal matrix of conditional variances, and λ is the unit risk premium.
This suggests that This implies that the risk-neutral dynamics are given by: The multivariate GARCH models are outlined in the next subsection.

Multivariate GARCH Models
According to Francq and Zakoian (2019), the CCC-GARCH model is formulated as follows: where R is the constant correlation matrix of j,T j (estimated using historical data). The diagonal matrix D t takes the following form The conditional variance of each asset is assumed to be consistent with Equation (1). An obvious shortcoming of the CCC-GARCH model is the assumption of constant correlation. To address this problem, Engle (2002) extended the model to incorporate dynamic conditional correlation (DCC). The DCC-GARCH model is formulated as follows: Q t is modelled using an autoregressive process: where v t = j,T j / h j,t ,Q is the unconditional covariance matrix of j,T j , and to ensure stationarity and positive definiteness, θ 1 + θ 2 < 1 and θ 1 , θ 2 > 0. The log-likelihood function (up to a constant) of both the CCC-GARCH and DCC-GARCH models is given by: In this study, we consider futures prices on Bitcoin with different expiry dates. Therefore, highly correlated asset price processes are expected (the same underlying ones). Given the volatility process, risk-neutral sample paths of Bitcoin futures prices can be simulated. The price of a Bitcoin futures spread option at time t, that expires at time T, is given by where T 1 , T 2 ≥ T, T 1 = T 2 , s K is the spread, and DF(t, T) is a discount factor used to discount a cashflow from time T to t (in this paper, we use the US 3 Month Treasury Bill rate as a proxy for the risk-free rate). It is clear from the above that when s K = 0, it is an exchange option. The empirical results are considered in the next section.

Empirical Results
In this section, the empirical results are presented and discussed. In this study, daily data 1 from 30 October 2020 to 1 April 2021 were used. The expiry dates of the Bitcoin futures prices are as follows: 30 April 2021, 28 May 2021, and 25 June 2021. The Bitcoin futures prices and returns are plotted in Figures 1 and 2   It is clear from the above that Bitcoin futures prices are trended. Furthermore, the returns show signs of volatility clustering, which is consistent with the typical stylised facts of financial returns (McNeil et al. 2015).
The descriptive statistics of the Bitcoin futures returns are reported in Table 1 below. The descriptive statistics indicate that the conditional expectation of the returns is close to zero, the returns are not normally distributed, and the returns also show signs of leptokurtosis. This is consistent with the stylised facts of financial returns (McNeil et al. 2015).
The estimated parameters (maximum-likelihood) and information criteria of the symmetric and asymmetric HN model applied to Bitcoin futures returns are reported in Tables 2 and 3 below. The AIC indicates that the symmetric time-varying volatility model is a better fitting model. This is consistent with previous findings in the literature (see, e.g., Conrad et al. 2018;Dyhrberg 2016). In addition to the AIC, a likelihood ratio test is also used to compare the estimated symmetric and asymmetric models. The test statistic of the likelihood ratio test for each expiry is reported in Table 4 below: We do not reject the null hypothesis for each expiry (Held and Sabanés Bové 2014), which is also in favor of the symmetric HN model.
The futures option prices of the two models and market prices 2 (scaled by the relevant futures price) are plotted in Figure 3 below. Furthermore, the pricing performance metrics of the two models applied to different futures are outlined in Tables 5-7  It is clear that similar European futures option prices are obtained when the different models are compared. The models produce reasonable prices compared to market prices. The RMSE and MAE of the symmetric HN model are slightly lower compared to the asymmetric HN model.
The performance metrics of the two models are similar. Therefore, to determine whether the predictive accuracy of the two models are the same, the test by Diebold and Mariano (1995) was applied. The test statistic for each expiry is reported in Table 8 below: The null hypothesis of the test is that the symmetric and asymmetric HN model have the same accuracy, and the alternative hypothesis is that the symmetric model outperforms the asymmetric model. The Diebold-Mariano test indicates that the symmetric model outperforms the asymmetric model for options that expire in April and May, but the predictive accuracy of the models is the same for options that expire in June.
Based on the pricing performance of univariate options, the symmetric HN GARCH process is used for the pricing of short-dated spread options on Bitcoin futures. For the CCC-GARCH model, no additional parameters need to be estimated. The additional parameters of the DCC-GARCH model are reported in Table 9 below: The CCC-GARCH and DCC-GARCH models are compared using a likelihood ratio test. The value of the test statistic is −3.8 × 10 −7 , which is insignificant. Hence, we do not reject the null hypothesis, which is in favour of the CCC-GARCH model. Therefore, the CCC-GARCH model is used for the pricing of Bitcoin futures spread options.
As mentioned previously, Bitcoin futures are highly correlated. To illustrate this concept, risk-neutral sample paths are illustrated in Figure 4 below. Spread options on Bitcoin futures do not actively trade and, therefore the model prices cannot be compared to market prices.

Conclusions
Bitcoin futures options were launched in the first quarter of 2020. GARCH modelling of Bitcoin returns and Bitcoin option pricing is well documented in the literature. However, the pricing of Bitcoin futures options in a GARCH framework has not been considered. In addition, a methodology for price discovery of multivariate options on Bitcoin futures has not been developed.
In this study, Bitcoin futures options were priced using the Heston-Nandi futures option pricing model (Li 2019a). The empirical results show that the symmetric Heston-Nandi model is a better-fitting model, which is consistent with previous studies that focused on Bitcoin spot return dynamics. The pricing performance metrics show that the Heston-Nandi model produces reasonable Bitcoin option prices, and that the symmetric Heston-Nandi model also produces more accurate option prices compared to market prices for two out of the three expiry dates considered.
In addition to the pricing of univariate Bitcoin futures options, the work by Rombouts and Stentoft (2011) was also extended to the pricing of multivariate futures options. The model was applied to Bitcoin futures spread options. The model produces reasonable spread option prices. However, spread options on Bitcoin futures do not actively trade. Therefore, the model prices cannot be compared to market prices.
The empirical results show that the symmetric Heston-Nandi model is more accurate in most cases. This implies that the symmetric model is a better choice when pricing exotic options (univariate) and other illiquid derivatives written on Bitcoin futures. Furthermore, the multivariate GARCH analysis showed that the CCC-GARCH model is more appropriate (based on historical data) when pricing multivariate Bitcoin futures options. Hence, the symmetric Heston-Nandi model and CCC-GARCH model can serve as a basis for pricing and risk measurement (quantifying market risk and capital calculations) of Bitcoin futures options.
Areas for future research include the use of skewness and kurtosis in the Heston-Nandi futures model (see e.g., Christoffersen et al. 2006) applied to univariate Bitcoin futures options, and the use of different multivariate GARCH processes (e.g., BEKK) applied to multivariate Bitcoin futures options. Furthermore, the hedge performance of the Heston-Nandi model applied to univariate and multivariate Bitcoin futures options should also be considered.