Bitcoin Volatility and Intrinsic Time Using Double Subordinated Levy Processes

We propose a doubly subordinated Levy process, NDIG, to model the time series properties of the cryptocurrency bitcoin. NDIG captures the skew and fat-tailed properties of bitcoin prices and gives rise to an arbitrage free, option pricing model. In this framework we derive two bitcoin volatility measures. The first combines NDIG option pricing with the Cboe VIX model to compute an implied volatility; the second uses the volatility of the unit time increment of the NDIG model. Both are compared to a volatility based upon historical standard deviation. With appropriate linear scaling, the NDIG process perfectly captures observed, in-sample, volatility.


Introduction
Cryptocurrencies have emerged as an alternative asset class in their own right.By now, institutional as well as private investors have taken a closer look at this asset class.Many view digital currencies, particularly the leading cryptocurrency, bitcoin, as an alternative to gold for providing long-term protection against inflation.In contrast, a major argument against investing in cryptocurrencies is their high and often erratic volatility when compared to conventional speculative assets.
Data on bitcoin volatility, in terms of the Bitcoin Volatility Index (BVI), have been available since 2010. 1 The BVI index is based on historical volatility, specifically, the standard deviation of past daily log-returns.Two versions of this index are provided, based on 30-and 60-day windows of past observations.The problem with historical volatility is that it is only meaningful if the returns are approximately independent and identically distributed (iid).In the presence of temporal dependence in the form of volatility clustering, a phenomenon common to virtually all speculative assets, including cryptocurrencies, historical volatility fails to adequately capture the current risk of an investment in a speculative asset.It simply provides an average of past risk.
Implied volatility measures attempt to tackle this shortcoming by capturing the current market view on future risk as implied by observed transaction prices of option contracts.Assuming an option-pricing method and given the price at which an option has been traded, one can back out the expected volatility of the underlying asset that is impled by the observed price.Several approaches have been developed to price options contracts, including continuum models such as the Black and Scholes (1973) and Merton (1973) (BSM) formula, as well as discrete models based on binomial or trinomial trees, Monte Carlo simulation, discrete stochastic volatility, and finite-differencing methods.Alexander and Imeraj (2021) use bitcoin options data from the Deribit exchange to create a term structure for bitcoin implied volatility indices using the same variance swap fair-value formula (geometric variance swap) that the Chicago Board Options Exchange (Cboe) uses to construct the VIX index (CBOE, 2019) that reflects the volatility of the S&P 500 stock index.Furthermore, Alexander and Imeraj point out that bitcoin prices appear to fluctuate excessively, so that this approach tends to underestimate the fair value of geometric variance swaps.In a recent report, Kim et al. (2019) proposed a volatility index based on the cryptocurrency index CRIX. 2It proxies the next month's mean annualized volatility.Venter and Maré (2020) used the symmetric GARCH option pricing model to construct various implied volatility indices, including the CRIX index.More recently, an alternative volatility index, BitVol (T3I-Pty-Ltd., 2019), has been introduced.It reflects the expected 30-day implied volatility derived from tradable bitcoin option prices using the BSM model.
The BSM model rests upon the assumption that the returns of the underlying asset follow a normal distribution.There is, however, overwhelming empirical evidence that return distributions exhibit heavy tails and asymmetry.Osterrieder (2016) demonstrated that these properties hold for bitcoin returns.The omnipresence of extreme return observations indicates that a heavy-tailed distribution should provide a more realistic description for the behavior of bitcoin returns.Heavy tails and asymmetry contradict the normality assumption and present a challenge to BSM option pricing.The Student's t-distribution would be an obvious candidate for heavy-tailed data.However, in the context of option pricing, use of the t-distribution in the BSM setting results in a divergent integral (Cassidy et al., 2010).Recently, Näf et al. (2019) proposed a mean-variance heterogeneous tails mixture distribution for modeling financial asset returns.They showed that the new model captures, along with the obligatory leptokurtosis, different tail behavior among the assets.
The method of subordination 3 offers an alternative approach to addressing non-normality.By adopting the concept of intrinsic or operational time rather than calendar time, the method of subordination allows the variance of the normal distribution to change over time and to capture heavy-tailed phenomena.The subordination approach leads to a generalization of the classical asset pricing model. 4Shirvani et al. (2021) introduced and analyzed various multiple subordination processes and demonstrated that the resulting dynamics lead to more realistic asset pricing models.The idea of multiple subordination is motivated both by the need to better capture the heavy-tailedness that conventional Lévy-subordinated models fail to explain and the desire to incorporate the views of investors in asset and option pricing models.
In this paper, we develop a bitcoin volatility index based on the Cboe variance swap fair-value formula but with log-prices following a double subordinated process.We briefly introduce the idea of double subordination to model the dynamics of Lévy processes.We show our model implies a heavy-tailed distribution that provides a better description of bitcoin log-prices than the Student's t model.
Rather than using Monte Carlo option pricing in non-Gaussian settings, as put forth in Allen et al. (2011), we employ the arbitrage theorem and mean-correction martingale measure (MCMM) of Yao et al. (2011) to price options when the returns of the underlying assets are doubly subordinated.Specifically, we price European contingent claims when bitcoin returns are assumed to be driven by a double subordinated Lévy process.We determine the equivalent martingale measure to price options using the MCMM approach and demonstrate that our proposed pricing model is arbitrage free.Based on generated bitcoin option data, we derive implied volatilities using the Cboe approach and the double subordination approach.
The remainder of the paper is organized as follows.Section 2 describes the double subordinated process and applies it to model bitcoin log-prices.Section 3 details the option pricing framework under double subordination and applies it to model bitcoin vanilla European call and put options.In Section 5 we assess bitcoin volatility using the Cboe approach, which is based upon option pricing, and our intrinsic time approach based upon spot pricing.Both methods are compared with historical volatililty.The paper concludes with a discussion of our results.

Double Subordination Model for Bitcoin Log-Prices
The double subordination framework adopted here involves a Lévy subordinator process.A stochastic process X = (X(t), t ≥ 0) defined on a stochastic basis (Ω, F, F, P) is said to be a Lévy process if the following conditions hold.
• X is a process with independent increments; for any partition 0 = t • X is a process with stationary increments; for any 0 ≤ s < t, the increment X (t)−X (s) has the same distribution as • X is continuous in probability; for every ε > 0 and t ≥ 0, there exists h ε,t > 0, such that In dynamic asset-pricing theory, a risky financial asset is defined by its price dynamics S t , t ∈ [0, τ ], where τ < ∞ is the time horizon; that is, τ is the maturity date of a financial contract.A Lévy process T = (T (t), t ≥ 0, T (0) = 0) with non-decreasing trajectories (i.e., non-decreasing sample paths) is called a Lévy subordinator.Since T (0) = 0, the trajectories of T take only non-negative values.In the BSM option pricing model, the price dynamics of the underlying asset are defined by where the log-process is and B = (B t , t ≥ 0) is a standard Brownian motion.To allow for non-normality of asset returns, Mandelbrot and Taylor (1967) and Clark (1973) suggested the use of a subordinated Brownian motion, where the price process S (ss) t and the log-price process are defined as where T = (T (t), t ≥ 0, T (0) = 0) is a Lévy subordinator.Process (4) describes the well-known single subordinated log-price process.Various studies have demonstrated that single subordinated log-price models commonly fail to capture the heavy-tailedness observed in financial return data.For example, Lundtofte and Wilhelmsson (2013) and Shirvani et al. (2021) showed that the normal-inverse Gaussian (NIG) distribution, which is a single subordinated Lévy process, fails to explain the equity premium puzzle.This is partly due to the fact that the tails produced by NIG are not heavy enough.However, Shirvani et al. (2021) demonstrated that the high value for the risk aversion coefficient, which gives rise to the equity premium puzzle, is compatible with a return process driven by a double subordinator model.Shirvani et al. (2021) defined and investigated the properties of various multiple subordinated log-return processes designed to model leptokurtic asset returns.They showed that multiple subordinated log-return processes can imply heavier tails than single subordinated models and that they are capable of capturing skewness and kurtosis.Therefore, a double subordination framework may be a more appropriate candidate for modeling the rather extreme behavior of bitcoin.
To apply double subordination to modeling the bitcoin price process, let S t denote the price process with the dynamics (5) ) where the triplet members {B s , T (s), U (s)}, s ≥ 0, are independent processes generating the stochastic basis (Ω, F, F = (F t , t ≥ 0) , P) representing the real world.Here, {B s , s ≥ 0} is a standard Brownian motion, while {T (s), s ≥ 0, T (0) = 0} and {U (s), s ≥ 0, U (0) = 0} are Lévy subordinators.B t , T (t) and U (t) are F t -adapted processes, whose trajectories are right-continuous with left limits.Shirvani et al. (2021) referred to T (U (t)), t ≥ 0 as the double subordinator process; hence, the process modeled by ( 6) is a double subordinated log-price process.
Consider the case where the subordinators T (t) and U (t) are inverse Gaussian (IG) Lévy processes; that is, T (1) ∼ IG(λ T , µ T ) having the probability density function (pdf), Similarly U (1) ∼ IG (λ U , µ U ).In this case, Shirvani et al. (2021) referred to T (U (t)) as the double inverse Gaussian subordinator and to X t as a normal double inverse Gaussian (NDIG) log-price process.The characteristic function (chf) of X 1 is given by The NDIG model 6 and 8 has eight parameters, namely µ 3 , σ 3 , γ, ρ, µ τ , λ τ , µ T and λ T , which can make fitting the model to data a challenging task.However, it is worth noting that that only six of these parameters are identifiable within the model.1To set the remaining two parameters, we can consider the expectation: As the processes U and T are independent inverse Gaussian (IG) processes, we can uniquely identify γ and ρ by requiring: By utilizing equation 10, the set of model parameters that can be identified becomes µ 3 , σ 3 , γ, ρ, λ T , and λ U .
The central moments of the NDIG can serve as a source of data for fitting these six parameters.The moment generating function (MGF) M X 1 (w) for X 1 , which generates the moments of its probability distribution, can be obtained by evaluating E e iwX 1 for w ∈ R.This can be obtained from equation 8 by setting w = iv.The MGF is written in terms of the cumulant generating function K(w): Using the identity we can express M Xt (w) as M Xt (w) = exp (K X 1 (w)t), where K X 1 (w) can be written using equation 11 as: Now, by finding the first four central derivatives of K X 1 (w), we can have the first four centered moments of X 1 .

Bitcoin Option Pricing under Double Subordination
To price European contingent claims, we assume that X t follows a NDIG log-price process.We need to derive an equivalent martingale measure (EMM) Q of P on (Ω, F, F = (F t , t ≥ 0) , P), such that the discounted price process, e −rt S t with r ≥ 0 denoting the riskless rate, is a martingale. 7To do so, we use the MCMM, since the proposed pricing model is arbitrage-free under MCMM, and estimate the parameters specifying the process.We then add the appropriate drift term to the process, such that the discounted price process becomes a martingale.Yao et al. (2011) constructed a martingale measure using MCMM for the geometric Lévy process model and showed that this measure is an EMM if there is a continuous Gaussian part in the Lévy process.In the case that X t is a pure jump Lévy process, they pointed out that this measure cannot be equivalent to a physical probability.However, pricing European options under this measure is still arbitrage free.
Let C be a European call option having an underlying risky asset S with a price and logprice process as in ( 5) and ( 6), respectively.Let B be a riskless asset with price b t = e rt , t ≥ 0, where r ≥ 0 is the riskless rate.Then, the price of C is where τ > 0 denotes the maturity, K > 0 is the strike price, and is the price dynamic of S on Q (an equivalent martingale measure of P).Using MCMM to find the EMM, describes the dynamics of S on Q, where M Xt (v) is the MGF of X t and K Xt (w) = ln M Xt (w) is the cumulant-generating function of X t .The chf of the log-price, ln where Carr and Madan (1998) showed how to use the fast Fourier transform (FFT) to value options in the case in which the chf of the log-price of the underlying asset is known analytically.They consider the modified option price c a (S 0 , r, k, τ ) = e ak C(S 0 , r, k, τ ), where k = ln(K), which, for a range of values a > 0, guarantees that c a (S 0 , r, k, τ ) is square integrable over k ∈ (−∞, ∞).The access point for applying the FFT is their derived relationship Numerical solution of this integral involves two fundamental concerns, an "optimum" value for a and control over the error produced by truncating the integral (18) over a finite domain [0, v max ].Addressing the first of these two concerns, Carr and Madan note that a positive value for a guarantees square-integrability of c a (S 0 , r, K, τ ) over the negative k axis but aggravates square-integrability over the positive k axis.They derive a sufficient condition, < ∞ which, combined with the analytic expression for the chf, can be used to determine an upper bound on a. Addressing the second of these concerns, they develop a bound on the error involved in truncating the integration range of (18), where ϵ is a desired level of truncation error and √ A/v 2 is an upper bound on the magnitude of the integrand in (18).
Constructing the implied volatility surface for call prices requires estimating (18) over a discrete mesh of (k, T ) values.Put options can then be valued assuming that put-call parity holds.Implementing a FFT requires numerical discretization of (18) into the form which the FFT can solve in O(M ln 2 (M )) operations.Writing (18) as discretization of Ψ(k) over a finite range [0, v max ] using the left hand rectangle rule 8 gives where which is identical to (19) with the identifications Ẑp = Ψ(k p ), Z j = exp(i kv j )h(vj)∆v, M = N , and 2π/M = ∆v∆k.This last identification gives the familiar FFT tradeoff between the span covered in the "space" domain and in the "frequency" domain, In our computations in section 4.1, we used N = 2 10 = 1024.Carr and Madan (1998) introduced the parameter a to ensure that the call pricing function 15 is square integrable as K → 0 (i.e., as k = lnK → −∞).They note that a sufficient condition for square integrability is provided by the requirement that ϕ ln(S From 8 and 11, note that ψ(iw) ≤ Kψ(w) for w ∈ R. Hence, 24 and To ensure that the cumulant generating function remains real-valued, requirement 25 can be reduced to positive argument requirements for the square root evaluations in 13 for where a ∈ (0, ∞).From 11, under the assumption that γ = 0, this can be further reduced to the following requirements: Solving the latter equation for h(w) yields Combining 27 with the first equation in 26 gives The quadratic equation has two roots, and Equation 28 is satisfied for w = 1 + a when Since a is a non-zero value, it follows that By utilizing 28, equation 30 can be rewritten as 4 Numerical computation In this section, we illustrate the double subordinated method of sections 2 and 3 and fit the NDIG model to the log-return Bitcoin time series.Figure 1 shows the daily bitcoin price and log-return, r t time series covering the period from July 19, 2010 to July 28, 2023.The return series consists of 4,026 observations, to which we fit the NDIG model specified by ( 5)-( 8).
The NDIG model is characterized by six parameters, θ = (µ 3 , λ U , µ U , λ T , µ T , σ 3 ).To further simplify the parameter set µ 3 , σ 3 , γ, ρ, λ U , and λ T , we assume that the subordinator U (t) is used to model the intrinsic time of the return process, while the subordinator T (t) is used to model the return skewness and heavy-tailed behavior.It is reasonable in this model to demand that there be no γU (t) term in 6, i.e., γ = 0.With this requirement, the first four moments of X become: As there is no analytical expression for the probability density function, we employ the method of moments and the empirical chf to estimate the parameters of the model.We follow Paulson et al. (1975) and Yu (2003) and use the fact that the pdf is the Fourier transform of the chf.So, the five parameters µ 3 , σ 3 , ρ, λ U , and λ T can be estimated via minimization where r t denotes the observed daily bitcoin return time series.The inclusion of the term ∆CF relies on the one-to-one correspondence between the cumulative distribution function and the characteristic function (since the probability density function is the inverse Fourier transform of the characteristic function).The integral for ∆CF can be estimated using the method described by Yu (2003).
The parameters are estimated by the method of moments and the empirical chf fitting, which is a method as efficient as likelihood-based methods (Yu, 2003).
The resulting parameter estimates are reported in Table 1.
Table 1 NDIG parameter estimates for the bitcoin log returns.For comparison, we also fit the normal and the Student's t distributions to the Bitcoin daily log-return values and estimated the distribution parameters using maximum likelihood.The empirical density, the fitted NDIG model, and the best-fit Student's t (df = 1.60), and normal (μ = 0.0032, σ = 0.0551) distribution densities are compared in Figure 2. The results confirm, as was noted in the introduction, that the normal distribution is not well suited for modeling asset returns.The graphs indicate that both the Student's t and NDIG models match the tails of empirical density rather well, while NDIG matches the skewness in the empirical distribution more evenly than does Student's t.The low estimate of df = 1.60 for the Student's t distribution implies that even second moments do not exist.Thus, the fat tails of the Student's t distribution would cause divergence of the BSM integral needed to evaluate price options.

NDIG Model for European Call Option Pricing
In this section, the NDIG model is used to price a European call option C with the underlying risky asset S being the return of Bitcoin.The dynamics of S under Q is given by equation We apply the NDIG Lévy model to the pricing of European plain vanilla bitcoin options.Let C be a European call option where the underlying risky asset, S, follows the log-price process given in (6).The dynamics of S on Q are given by ( 16) and the chf of the log-price by ( 17).We evaluate the integral (18) using FFT for a range of strike levels and maturity horizons with the parameter estimates reported in Table 1.When using the global parameter values provided in Table 1, the upper bound of a in 32 is ≲ 21.93.
To ensure that option prices are within acceptable limits, a more stringent determination of the value of σ ∈ (0, ∞) is required.Specifically, for a European call option, the upper and lower bounds on the price are given by max(S − Ke −r f T, , 0) ≤ C(S t , T, K) ≤ S T .2 As noted above, an upper bound value for a and hence choice of a suitable value in performing the integrations via FFT is very dependent on the data set (the price process).When a > 1 the call option surface was unstable.We follow Schoutens (2003) and set a = 0.15.
Figures 3 (a) and (b) show the resulting prices for call and put options plotted against maturity, τ , and strike price K. Figures (c) and (d) show the implied volatility surfaces, i.e., the market's view on future volatility, against time to maturity and relative moneyness, defined as K/S, where S is the asset price obtained by the BSM model.Figure 3-(c) investigate the behavior of implied volatility concerning moneyness and time to maturity for both call and put options.It is observed that as moneyness increases, moving towards out-of-the-money values, the implied volatility also increases.Conversely, as moneyness decreases, a distinct pattern emerges, referred to as the "smirk" in the implied volatility curve.Initially, as moneyness decreases and options move in-the-money, implied volatility rises, but subsequently, it starts to decrease for further in-the-money options, eventually approaching near-zero values.As we extend the time to maturity (denoted as τ ) of the options, we notice a noteworthy trend.The sensitivity of implied volatility to changes in moneyness diminishes for call options, and this effect appears to converge to a constant level irrespective of the moneyness.In other words, as the options approach their expiration dates, their implied volatility becomes less responsive to changes in moneyness.

Bitcoin Volatility Measures
Understanding the volatility of speculative assets is critical for investment decisions.Given that bitcoin is, at least by some, considered to be a potential alternative to fiat money, volatility characteristics are of especial concern.It is, therefore, of interest to understand and adequately model the process that governs the volatility of bitcoin.
Model parameters reflecting the volatility of asset prices are known to vary stochastically and exhibit clustering.To allow for these phenomena in option pricing, Hull and White (1987) and Heston (1993) suitably randomize the volatility parameter in the Black-Scholes model, where the volatility process is governed by Brownian motion.An alternative strategy was adopted by Carr et al. (2003) based upon ideas by Geman et al. (2001) and Clark (1973).Clark conjectured that price processes are controlled by a random clock, which is a cumulative measure of economic activity, and used transaction volume as a proxy for this measure.Hurst et al. (1997) considered various subordinated processes to model the leptokurtic characteristics of stock-index returns.In the option pricing literature, Carr and Wu (2004) extended the approach in Carr et al. (2003) by providing an efficient way to allow for correlation between the stock price process and random time changes.Klingler et al. (2013) introduced two new six-parameter processes based on time changing tempered stable distributions and developed an option pricing model based on these processes.Shirvani et al. (2021) introduced the volatility intrinsic time or volatility subordinator model to reflect the heavy-tail phenomena present in asset returns.They studied the question of whether the VIX index is a volatility index that adequately reflects intrinsic time.They showed that this index fails to properly capture the intrinsic time for the SPDR S&P 500.Apparently, the VIX index, as a measure of time change, does not reflect all the information needed to correctly capture the skewness and the fat-tailedness of the S&P 500 index.A model with a suitable volatility subordinator should adequately account for such empirical phenomena.
In this section, we apply three views to analyze bitcoin volatility.The first measures historical realized volatility through sample standard deviation.As noted in the Introduction, this is the method employed in the Bitcoin Volatility Index.The second measures implied (future) volatility in the risk-neutral, derivatives world, Q, reflecting the views of the option traders.The third measures implied volatility in the real world, P, reflecting the views of spot traders.We employ an intrinsic time formulation in ths latter case.In the following, we empirically compare these three modeling strategies using daily bitcoin data from January 1, 2015 to April 8, 2021.
Historical Volatility. Figure 7 plots the historical volatility series based on measuring average standard deviation on a 1008-day rolling window of bitcoin returns.The historical volatility varies substantially over the sample period, assuming values between 60.69 and 159.62 percent.
Implied Volatility in Q. Option traders commonly use implied volatility as a proxy for future volatility.The VIX, the index reflecting volatility expectations for the S&P 500 stock index, is based on S&P 500 index options estimated on a real-time basis by the Cboe.The VIX can be regarded as an efficient volatility forecast for the S&P 500 index, provided that option markets are efficient.
The value of the VIX index is Subscript "1" denotes near-term options while "2" denotes next term options.As SPX options expire on Fridays, depending on the trading day, near-term options have 23 to 30 days until expiration, while next-term options have 31 to 37 days until expiration.The weights W j (j = 1, 2) express these expiration times in a normalized form, accurate to the minute, where M j represents the number of minutes until the settlement of the j-th term options, while M 30 represents the total number of minutes in a 30-day period.As a result, W 1 and W 2 are constrained to the range 0 ≤ W 1 , W 2 ≤ 1, with the additional requirement that W 1 + W 2 = 1.The equation below, provided by Demeterfi et al. (1999), can be used to calculate the near-term and next-term volatilities: Here, T j is the time to expiration measured in fractions of a year, r f is the annualized riskfree rate, F j is the forward index level derived from index option prices, K 0j is the first strike price below F j , K i is the strike price of the ith out-of-the-money option (a call if K i > K 0 , a put if K i < K 0 , and both call and put if K i = K 0 ), ∆K i is the interval between strike prices, and Q K i is the midpoint of the bid-ask spread for each option with strike price K i .For more information on the computation of each of these parameters, refer to the Cboe white paper by CBOE (2019) on VIX.
To calculate a "VIX-like" volatility known as BVIX, the NDIG model is used to compute prices for European put and call options that have between 23 and 37 days to expiration.To compare with the historical volatility, BVIX values are calculated using historical data as follows: for each moving window of length 1008 Bitcoin log-return data, 1.The NDIG model is fitted to the 1008 Bitcoin log-return data (4 years), and the model parameters are estimated using the minimization 33.The resulting parameter values are displayed in Fig. 4. The results depicted unveil discernible patterns among the parameters under investigation.Particularly noteworthy is the conspicuous observation that λ T displays the least temporal variation, while ρ showcases the highest level of temporal variability.Furthermore, an examination of the estimated values for µ 3 and σ 3 indicates a declining trend leading up to April 2017, followed by a subsequent abrupt change, with relatively minor fluctuations thereafter.Additionally, it is worth noting that the parameter ρ exhibits the smallest magnitude in its estimated value, accompanied by a considerable degree of dispersion or deviation from the mean.These observations are graphically represented in Figure 4, which presents a comparative visualization of the first four empirical moments derived from the time series depicted in Figure 4 .These moments are contrasted with those computed from the parameters derived via Equation 33.This analysis offers valuable insights into the dynamics and characteristics of the examined parameters.
Fig. 5 compares the first four empirical moments computed from the time series in Fig. 1 and the moments computed from the fitted parameters using 32. 2. The dynamics of the Bitcoin price S Q T in the risk-neutral world and subsequently, the characteristic function (CF) of the log price ln S Q T are determined using 17.
3. Call option prices with appropriate expiration dates are computed using the FFT formulation 18, and put option prices are computed using put-call parity based on the portfolio spot price S Q T for the last date of the moving window.In the context of option price calculations based on 18, the numerical values for the parameters a, a max , and v max are determined through the utilization of conditions outlined in ?? and Equation ??, respectively.Choosing too small values for a results in option prices that exceed the maximum threshold, while choosing a too large produces prices that go below the minimum threshold (and maybe negative).As we discussed, a stricter determination for the value of a ∈ [0, a max ] is based on the requirement that option prices remain within allowed bounds.Fig. 6 shows the computed values for the upper bound a max over period 05/08/2014 through 07/28/2023.For the computations in the remainder of this section, we utilize the smallest possible value, a = 0.40.(e) The range of strike prices considered in the BVIX computation is from 0.75S t to 1.5S, as NDIG computed option prices do not go to zero.
As mentioned in step 4d, since the underlying asset portfolio is not traded, there is no market sentiment setting bid-ask prices for the call and put options used in the BVIX calculation.Therefore, the BVIX volatility directly reflects the NDIG option price computations, and there is no risk premium or market price of risk that is typically captured by the VIX.
Figure 7 shows the values for BVIX.As observed from this figure, BVIX values are higher than STD values.The fact that implied volatility computed from options based upon a risky asset is higher than realized volatility is commonly observed and is referred to as the "volatility risk premium."An explanation for the volatility risk premium is that option traders allow for the probability of a significant market move.To be compensated for the inherent volatility risk that cannot be hedged away, option sellers charge a premium.Option traders typically view implied volatility as the most reliable approach to volatility forecasting (Jorion, 1995).Figure 7 and a Pearson correlation coefficient of ρ = 0.48 between BVIX and historical volatility indicates that implied volatility provides only partial information about future empirical volatility.
To account for the different scales, the historical and BVIX volatility time series can be normalized separately by subtracting the series mean and dividing by the series standard deviation.These normalized plots, as shown in Fig. 8, demonstrate close agreement between the two, with the BVIX volatility displaying a slight daily fluctuation that is not present in the historical volatility.This emphasizes the fact that the sole source of variance in this numerical example is the daily returns of the underlying risky asset, the Bitcoin data, which is directly quantified by the historical volatility.Since all option prices in this example are calculated using the Carr-Madan formulation and the NDIG model, no additional volatility is introduced to these option prices through, for instance, trader sentiment (bid-ask pricing).Therefore, the BVIX volatility formulation only captures the original asset return volatility.Implied Volatility in P: Intrinsic Time.Mandelbrot and Taylor (1967) coined the term "intrinsic time" in finance to describe how market information is not continuously available, but instead arrives in discrete events occurring at varying intervals.Market orders for assets, for instance, are a prime example of such events, which can differ significantly in both timing and frequency throughout the trading day.These events provide information on asset value (price), and their occurrence, magnitude, and sign characterize intrinsic time.Price change information may also differ in its level of informativeness, with larger changes providing more information than smaller ones, and consecutive price changes of the same sign providing more information than those with opposite signs.Researchers have explored the concept of intrinsic time and applied it to financial time series in various works, including Guillaume et al. (1997), Tsang (2010), andAloud et al. (2012).
According to the intrinsic time perspective, no information is available between events, and therefore no time has elapsed.An analogy to conceptualize intrinsic time is to imagine an analog clock's seconds-hand that ticks intermittently and randomly, varying in volume with each tick.Time only moves forward when the hand advances, and the tick's volume represents the amount of information conveyed.
By utilizing the double subordinated method, a novel measure of volatility can be derived under the P measure.From (5), we view X t as the process governing the intrinsic time of the double subordinated price process.From (32) with γ = 0, the standard deviation (i.e. the volatility) of the unit increment of X t , when X t follows a NDIG log-price process, is given by Vol where λ T , λ U , and σ 3 are model parameters as defined in Table 1.We can refer to the volatility Var(X 1 ) as the NDIG volatility, which takes into account both the Brownian motion (Gaussian) and Lévy subordinator components of the model.However, it is important to note that this measure differs from other measures of volatility in the literature.
To derive the intrinsic time for the bitcoin log-price process, we compute rolling NDIG parameter estimates from the recent six years of daily data, using overlapping windows with a length of 1008 days.We use (37) to derive the annualized volatility implied by the NDIG model.The volatility computed by the NDIG intrinsic time model (NDIG IT) is also shown in Figure 7. Examination of all three volatility plots reveals substantial similarity between the patterns of NDIG intrinsic time and historical volatility, confirming that (37) captures realistic bitcoin dispersion patterns.We have addressed above the observed differences between the BVIX implied volatility and historical volatility.We note that the NDIG IT volatility predictions are consistently larger than historical, but otherwise demonstrate very strong synchronicity.Normalized 9 versions of the three series are plotted in Figure 8.The normalized NDIG IT and STD time series are now virtually identical, indicating an almost perfect linear correlation between the two.The BVIX variant is less synchronous but still co-moves to a certain degree, with a strong pattern of anticipating cycles of high volatility.Cointegration analysis is one approach to assess whether or not there are significant deviations in co-movement among time series.To apply cointegration, we first check whether each volatility series is integrated.The results of both the augmented Dickey-Fuller and KPSS tests (Kwiatkowski et al., 1992) shown in Table 2 support the hypothesis that each series is integrated of order one (has a unit root).Performing a Johansen test (Johansen, 1988) involving all three time series (Table 2) indicates a cointegrating relationship of orders 1 and 2 among the three time series.
Taken together, the results from Figure 8 and Table 2 show that the dynamic NDIG IT model, which is consistent with dynamic asset pricing theory and has predictive capabilities, provides very good agreement when applied "in-sample".The implied volatility model, BVIX, which is also dynamic and predictive, has poorer agreement with the in-sample data.The pattern of BVIX often, though not always, anticipating cycles of high volatility may speak to option traders' heightened sensitivity to potential increases in market volatility.
Finally, we investigate the question of whether there is long-range dependency (LRD) in bitcoin volatility, as has been found in stock volatility (see, for example, Asai et al., 2012).We fit a fractionally integrated autoregressive conditional heteroskedastic (FIGARCH) model to investigate both short-and long-range dependency in the bitcoin variance series.Specifically, we fit an AR(1)-FIGARCH(1,d,1) 10 with Student's t innovations.Fractional integration in the GARCH part, i.e., d > 0, indicates LRD in volatility and non-stationarity.For d < 0 we have anti-persistency and short-range memory.The AR(1)-FIGARCH(1,d,1) model parameters are estimated in a rolling fashion using a window of 252-day window.Figure 9 plots the sequence of d estimates in the conditional heteroskedastic equation for the variance.The estimates fall almost exclusively in the interval (0.8, 1.0), indicating nonstationarity and the presence of LRD in bitcoin volatility.This implies that volatility measures derived from historical data can be used to determine bitcoin volatility-a finding that also follows from Figure 8, showing that the (normalized) historical standard deviation and NDIG intrinsic time move closely together.

Discussion
Dealing with heavy-tailed asset returns remains a challenge in theoretical and empirical finance.Implicit in our approach is that there is a relationship between arrival times of financial news, the "intrinsic time" of finance, and price volatility for an asset, and particularly for bitcoin, whose volatility is especially pronounced.As noted by Clark (1973), "The different evolution of price series on different days is due to the fact that information is available to traders at a varying rate.On days when no new information is available, trading is slow, and the price process evolves slowly.On days when new information violates old expectations, trading is brisk, and the price process evolves much faster."Our analysis is based on the assumption that information (whether unexpected or reinforcing) arrives at a rate that defines an "intrinsic" or "operational" time scale.Whether this intrinsic time is measured in rate (ticks) or volumes of transactions, or another proxy is irrelevant; we take as a stylized fact that the intrinsic time distribution of information arrival is skewed and extremely heavy tailed.(See, for example, ? on the microstructure of the time between successive transactions.) The philosophy behind our approach is therefore as follows.The log-price process is additive Brownian motion (Equation (2)), but the time variable that drives the Brownian motion is not the physical time t as suggested by ( 2), but is an intrinsic or operational time.We utilize subordinated Lévy processes to mimic a relationship between intrinsic time and price processes observed in physical time.Unfortunately, as has been demonstrated by ?, a single level of subordination (equation ( 4)) is insufficient to capture both the stylized facts of the operational time scale and the stylized facts (skewness, heavy-tails, clustering) observed in asset returns.In this paper we show that, for a very volatile asset, a doubly subordinated process (equation ( 6)) using a normal inverse Gaussian process at each level of subordination is capable of capturing (Figure 2) the skewness and heavy-tailed behavior of the asset's return process in physical time, and is, presumably also accounting for the

Figure 2
Figure 2 (a) Comparison of normal, Student's t, and NDIG distribution fits to the kernal density of the bitcoin return.(b) Same plot with the y-axis now in natural log scale.

Figure 3
Figure 3 NDIG-based price surfaces for bitcoin (a) call and (b) put options as functions of time to maturity and strike price.(c) Implied volatility surfaces as functions of time to maturity and moneyness.

Figure 5
Figure5Comparison of the first four empirical ("emp") moments computed from the return time series and the moments ("th") computed from the fitted parameters.

Figure 8
Figure 8 Normalized (mean subtraced, standard deviation scaled) versions of the bitcoin volatility series in Figure 7.

Figure 9
Figure 9 Time-varying LRD (d) in the bitcoin volatility log-return process.
Geman et al. suggested the price process must have a jump component; thus the price process can be regarded as Brownian motion subordinated to this random clock.Carr et al. developed this subordinated model using normal inverse Gaussian and variance gamma examples of pure jump Lévy processes.