A first attempt at pricing cryptocurrency options: a direct modelling approach

A first attempt at pricing cryptocurrency options: a

direct modelling approach

J.W. Lutterop

s2226243

A first attempt at pricing cryptocurrency options: a

direct modelling approach

J.W. Lutterop

Abstract

Contents

1 Introduction 2

2 Theoretical background 4

2.1 ARMA-GARCH models . . . 4

2.1.1 Standard GARCH model (sGARCH) . . . 4

2.1.2 Integrated GARCH model (IGARCH) . . . 5

2.1.3 Exponential GARCH model (EGARCH) . . . 5

2.1.4 Glosten-Jagannathan-Runkle GARCH model (GJR-GARCH) . . . 5

2.1.5 Component standard GARCH model (CS-GARCH) . . . 6

2.1.6 The asymmetric power ARCH model (APARCH) . . . 6

2.1.7 The threshold GARCH model (TGARCH) . . . 6

2.1.8 Distribution of the innovations . . . 6

2.1.9 Comparing nested models . . . 8

2.2 Econometric option pricing . . . 8

2.2.1 _{Historical P-dynamics . . . .} 9

2.2.2 The stochastic discount factor . . . 9

2.2.3 _{RN Q-dynamics . . . .} 9

2.2.4 Pricing strategies . . . 10

3 Data source and description 11 4 Model specifications 13 4.1 Selection procedure and resulting models . . . 13

4.2 Discussion and concluding remarks . . . 14

5 Option pricing 16 5.1 Methodology . . . 16

5.1.1 _{Obtaining the Q-dynamics: the GHD and NIGD case . . . .} 16

5.1.2 _{Obtaining the Q-dynamics: the sGED case . . . .} 18

5.2 Implementation and results . . . 20

6 Concluding remarks 22

7 References 23

A Appendix: Tables and figures 25

1

Introduction

The launch of the Bitcoin network by Satoshi Nakamoto, marked the inception of a new type of “currency ”, the so-called cryptocurrency, where crypto refers to the the fact that it relies upon cryptography. The Bitcoin network was introduced as an electronic currency using a purely peer-to-peer system which would allow parties to conduct financial transac-tions without mediation of a third party, i.e. a financial institution (Nakamoto [2008] [41]]). The driving technology behind cryptocurrencies is the so-called blockchain, which guarantees the validity and security of the transactions made 1_{. The concept of a cryptocurrency goes} further back than 2009, electronic currency systems were already mentioned in Chaum [1983] [16]. The practical implementation of the concept, however, had not been realised until 2009, due to the technical difficulties surrounding decentralised transactions such as the validation process of transactions. A typical example of such a technical anomaly is double spending, a potential flaw in the network that enables one to spend a digital token more than once. The introduction of the blockchain technology in 2009 made the practical implementation of decentralised cryptocurrencies possible and sparked the inception of many more cryptocur-rencies in the years to come.

In present days, the term cryptocurrency is well known. Despite the availability of more than 1600 alternative cryptocurrencies, Bitcoin is still by far the largest in terms of market capitalisation, which is roughly twice as high as the market capitalisation of the number two, Ethereum 2_{. Even though cryptocurrencies are well-established and known to the greater} public, this was not the case in the years following the launch of bitcoin. In fact, bitcoin led a relatively anonymous existence in the years subsequent to its creation, it started getting more public attention throughout 2013. Near the end of the year bitcoin had reached a record value of approximately 1200 USD, compared to a value of roughly 10 USD at the beginning of the year. Cryptocurrencies have endured an abundance of controversy in the years there-after. Critics have compared cryptocurrencies to Ponzi schemes 3_{, pyramid schemes}4 _and economic bubbles 5_{. The increased popularity of cryptocurrencies attracted interest from} academic fields other than computer science, including finance and economics (e.g. Dwyer [2015] [25], Ciaian et al. [2016] [21], Cheah and Fry [2015 [15]]), law (e.g Ly [2013] [36], Dion [2013] [24]) and ethics (e.g. Angel and McCabe [2015] [1]). The peak of the cryptocurrency popularity was in the last months of 2017; in November 2017, a surge of interest in the cryp-tocurrency market resulted in steep increases of several crypcryp-tocurrency prices (most notably bitcoin) during the months of November and December, after which the prices plummeted at a similar pace (see figure 2).

The aforementioned turbulence in cryptocurrency prices certainly underlines the topical rel-evance of the cryptocurrency market. Moreover, it caused cryptocurrencies to be widely regarded as investment assets, rather than a means of payment 6_{. Unsurprisingly, markets} for futures and derivatives on cryptocurrencies are emerging gradually (e.g. on BitMEX). Due to these recent developments and the immature state of the current acedemic literature, cryptocurrencies form a fruitful field of research. The aim of this study is to shed light on modelling surrounding cryptocurrencies - by looking at the five biggest (in terms of market cap) cryptocurrencies and a cryptocurrency index, of which further details shall be provided in section 3 - and use that information to take a first step towards option pricing of plain vanilla European call options on cryptocurrencies. Firstly, the time series properties and other characteristics of the prices and log returns on the cryptocurrencies are investigated. Secondly, an algorithm is proposed and applied that selects a suitable ARMA-GARCH type model. Finally, the ARMA-GARCH model specification is used to price plain vanilla Euro-pean call options, which are subsequently compared to option prices with more conventional benchmark pricing methods.

Similar studies are scarce, but do exist in the current literature. However, the majority of the similar literature encompasses the forecasting ability of GARCH or ARIMA type

1_{The mechanical details surrounding cryptocurrency related technology are beyond the scope of this paper}

and shall not be discussed here.

2_{See E.g. www.coinmarketcap.com.}

3

e.g. head of the World Bank Jim Yong Kim, February 2018

4_{E.g. Dan McCrum in the Financial times, November 2015.}

5

E.g. Robert Schiller in the New York Times, March 2014.

6_{Though, famous investor Warren Buffet has cautioned people not to look at cryptocurrencies as an}

models for cryptocurrencies. For instance Catania et al. [2018] [14] test the forecasting abil-ity of future volatilabil-ity of the standard GARCH model (see subsection 2.1.1) against other “more refined” models, though it is generally outperformed by the other tested models. Their findings further indicate that incorporating leverage and time varying skewness into the model specification can improve volatility predictions and therefore could be useful for pricing derivatives. Chu et al. [2017] [20] estimate twelve different GARCH models on seven different cryptocurrencies and compare the models in terms of in-sample fit, forecasting abil-ity and acceptabilabil-ity of value at risk estimates. Their model specifications are flexible in the sense that the distribution of the innovations (see subsection2.1.8 for more details) can be one of eight different distribution choices. However, the GARCH orders chosen are not flexible as they are fixed at (1,1) (see subsection 2.1 for more information on GARCH orders). They compare the in-sample fit in terms of information criteria and conclude that the GARCH models with normal innovations usually have the best fit. Furthermore, they conclude that the IGARCH (see subsection 2.1.2) and the GJR-GARCH (see subsection 2.1.4) models perform best in terms of in sample fit. Similar GARCH modelling oriented studies include e.g. Chen et al. [2017] [17] and Dyhrberg and Haubo [2016] [26]. Option pricing methods have barely been developed in the literature due to the absence of an established market for cryptocurrency options. Chen et al. [2018] [18] develop an option pricing method for CRIX and Bitcoin in continuous time. They develop a stochastic volatility model with correlated jumps and calibrate it using a Markov Chain Monte Carlo method. They conclude that their method produces stable results. Other than the aforementioned paper, the literature is practically non-existent, though I shall not exclude the possibility that other researchers have ventured in this direction.

2

Theoretical background

This section provides a theoretical background of the econometric ingredients used in this study. In subsection 2.1 the discrete time models used, specifically ARMA-GARCH type models, are discussed and subsection 2.2 describes the mathematical tools for (discrete time) econometric option pricing.

2.1

ARMA-GARCH models

The econometric models used in this study to model the log returns of the cryptocurrency prices are ARMA-GARCH type models, where ARMA stands for autoregressive moving average and GARCH stands for generalised autoregressive conditional heteroskedasticity. The ARMA component of the model models the log returns of the cryptocurrencies. ARMA models are nowadays well established in the econometric literature and are of the form

rt= µ + p X j=1 φjrt−j+ q X j=0 θjεt−j (1)

where rtis the log return of the cryptocurrency, p is the autoregressive order, q is the moving average order, the εt’s are assumed to be zero mean white noise error terms and θ0 ≡ 0. The GARCH component models the volatility of the return series. GARCH models are generalisations of the ARCH model proposed by Engle [1982] [28], which is given by

εt= σtzt zt i.i.d. ∼ (0, 1) σ2_t = ω + n X j=1 αjε2t−j = ω + α(L)ε2_t (2)

where the lag polynomial α(L) = Pn

j=1αjLj and L is the lag operator. (1) and (2) can respectively be seen as the conditional mean equation and the conditional variance equation of the log return process. The addition of the ARCH conditional variance equation given above allows one to relax the usually inaccurate assumption of having constant volatility over time.

The models used in this study all have a conditional mean equation of the form (1), but have different conditional variance equations, which are all generalisations of (2) and modifi-cations thereof. Consequently, the rest of the subsection proceeds by describing the volatility models used in this study.

2.1.1

Standard GARCH model (sGARCH)

The standard GARCH model was proposed by Bollerslev [1986] [9] and it generalises the conditional volatility equation 7of (2) by adding autoregressive terms:

σ2t = ω + n X j=1 αjε2t−j+ m X j=1 βjσ2t−j = ω + α(L)ε2t+ β(L)σt2. (3)

where the lag polynomial β(L) =Pm

j=1βjLj. The attractive feature of this generalisation is that the model, in addition to non-constant volatility, also captures volatility clustering, i.e. the tendency that periods of high volatility tend to be followed by periods with high volatility and vice-versa. Clearly, high values of βj are associated with high persistence in volatility, i.e. high volatility clustering.

Note that (3) only makes sense if σ2

t is non-negative almost surely. Bollerslev [1986] [9] im-poses the restrictions that ω ≥ 0, αj≥ 0 ∀1 ≤ j ≤ n and βj ≥ 0 ∀1 ≤ j ≤ m, which clearly ensures this condition. Nelson and Cao [1992] [43] offer a less restrictive set of constraints,

7_{In fact all of the different GARCH models described in this section only differ in the conditional volatility}

partially based on rewriting the GARCH(n,m) specification to an ARCH(∞) specification. Rewriting the conditional volatility equation corresponding to (3) gives

σ2_t = [1 − β(L)]−1[ω + α(L)]ε2_t = ω∗+ ∞ X k=0 ψkε2t−k−1.

Necessary conditions for the non-negativity of σ2

t are that ω∗ and {ψk}∞k=0are non-negative. In order for them to be well defined one needs the following two sufficient conditions. i. The roots of 1 − β(x) = 0 lie outside the unit circle.

ii. 1 − β(x) = 0 and α(x) = 0 have no common roots. The above conditions, however, do not guarantee that σ2

t is finite with probability one, nor that it is strictly stationary. To guarantee this one needs further restrictions, see Nelson and Cao [1992] for further details. A sufficient, but non-necessary condition for strict stationarity, proposed by Bollerslev [1982] [9] is that

α(1) + β(1) = n X j=1 αj+ m X j=1 βj< 1. (4)

2.1.2

Integrated GARCH model (IGARCH)

Integrated GARCH models, first discussed by Engle and Bollerslev [1986] [29], arise when condition (4) is not satisfied, in particular when α(1) + β(1) = 1. As an illustration, the IGARCH(1,1) model is characterised by the fact that α1+ β1 = 1 and consequently the conditional variance equation can be written as

σ2_t = ω + α1ε2t−1+ (1 − α1)σt−12 .

The property that α(1) + β(1) = 1 is commonly referred to as persistent variance.

2.1.3

Exponential GARCH model (EGARCH)

Nelson [1991] [42] proposes a modification of the sGARCH models due to three perceived drawbacks of the sGARCH approach: 1. evidence in the literature of a negative correlation between current returns and future returns volatility, 2. parameter restrictions imposed by GARCH models are often violated by coefficient estimates and may unduly restrict the dy-namics of the conditional volatility process, 3. investigating the persistence in the conditional variance is difficult as different norms for persistence often disagree. The EGARCH model is proposed to meet these objections, its conditional variance equation is given by

ln σ2_t = ω + n X j=1 [αjεt−j+ γj(|εt−j| − E|εt−j|)] + m X j=1 βjln σt−j2 . (5)

In this setup, σt2depends not only on the size (through αj) of {εt−j}nj=1, but also on the sign (through γj). Moreover, since the left hand side of (5) is expressed in natural logarithms, the right hand side is allowed to be negative, hence there are no parameter restrictions on the coefficients. Furthermore, the EGARCH model allows for the leverage effect, the phenomenon where the volatility of an asset is negatively correlated with its returns.

2.1.4

Glosten-Jagannathan-Runkle GARCH model (GJR-GARCH)

The GJR-GARCH model, named after its inventors (Glosten et al. [1993] [32]), has a con-ditional variance equation given by

σ2_t = ω + n X j=1 αjε2t−j+ γj1{εt−j< 0}ε2t−j + m X j=1 βjσ2t−j (6)

where1{·} is the indicator function and ω, the αj’s and βj’s are restricted to be non-negative. The negative sign effect, governed through γj, was motivated by empirical evidence of a negative relationship between conditional expected returns and conditional variance. The GJR-GARCH model is like the EGARCH model an example of an asymmetric modification of the sGARCH model. Here asymmetry refers to the notion that the response of σ2

2.1.5

Component standard GARCH model (CS-GARCH)

The CS-GARCH model, proposed by Engle and Lee [1999] [30] decomposes the conditional variance into a permanent (long-run) and a transitory (short-run) component. Its conditional variance equation is given by

σ2_t = qt+ n X j=1 αj(ε2t−j− qt−j) + m X j=1 βj(σt−j2 − qt−j) qt= ω + ρqt−1+ φ(ε2t−1− σ 2 t−1).

A key characteristic of the CS-GARCH model is that the intercept of the conditional variance equation is no longer constant over time (compared to the sGARCH model).

2.1.6

The asymmetric power ARCH model (APARCH)

The APARCH model, introduced by Ding et al. [1993] [23], was motivated by the obser-vation that the autocorrelations of the absolute log returns are higher than those of the squared log returns. Phrased differently, the autocorrelation functions of the absolute power transformation of the log returns given by |rt|δ typically peaks around δ = 1. This property is called the Taylor effect, named after Taylor [1986] [45]. Consequently, the authors propose a model which enables one to look beyond just conditional variance process:

σδ_t = ω + n X j=1 αj(|εt−j| − γjεt−j)δ+ m X j=1 βjσδt−j

where ω > 0, δ > 0, the αj’s and the βj’s are restricted to be non-negative and |γj| < 1. Note that the APARCH model nests several other GARCH specifications. If γj = 2 and δ = 0 we have the sGARCH model. If γj 6= 0 and δ = 2 we have the GJR-GARCH model and if γj6= 0 and δ = 1 we have the TGARCH model (defined below).

2.1.7

The threshold GARCH model (TGARCH)

Finally, the TGARCH model, introduced by Zakoian [1994] [47], bears similarities to the GJR-GARCH model, but is characterised by its conditional standard deviation given by,

σt= ω + n X j=1 α− j1{εt−j < 0}εt−j+ α+j1{εt−j≥ 0}εt−j + m X j=1 βjσt−j.

An advantage of this approach is that no non-negativity constraints are required, since σ2 t will be non-negative by construction. However, the study of the probabilistic properties of σ2

t becomes rather complicated if σt is not restricted to be non-negative, for specifics see Zakoian [1994] [47].

Since the introduction of the ARCH model, countless variations have been proposed. The list provided above is nowhere near complete, but does list the most well established variations in the current literature and will suffice for the purposes of this study. For a more complete list, the interested reader can consult Bollerslev [2009] [11].

2.1.8

Distribution of the innovations

What has been discussed so far is that the GARCH models used follow from (2) with dif-ferent specifications for σ2

introduced by barndorff-Nielson [1977] [4], which is characterised by fGH(y | λ, α, β, δ, µ) = (pα2_{− β}2_/δ)λ √ 2πKλ(δ p α2_{− β}2₎exp{β(y − µ)} Kλ−1/2  αpδ2_{+ (y − µ)}2₎ (p(δ2_{+ (y − µ)}2_)/α)1/2−λ Kλ(x) = 1 2 Z ∞ 0 zλ−1exp x 2  z +1 z  dz , x > 0 Y = R Θ =(λ, α, β, δ, µ)0 _{∈ R}5_{: δ > 0, α > |β| > 0 .} (7) Its mean and variance - along with the normalisation conditions - are given by

E(y) = µ + δβKλ+1(δ p α2_{− β}2₎ p α2_{− β}2_K λ(δ p α2_{− β}2₎ = 0 Var(y) = δKλ+1(δ p α2_{− β}2₎ p α2_{− β}2_K λ(δ p α2_{− β}2₎+ β2_δ2 p α2_{− β}2 ( Kλ+2(δ p α2_{− β}2₎ Kλ( p α2_{− β}2₎ − K_λ+12 (δpα2_{− β}2₎ K2 λ(δ p α2_{− β}2₎ ) = 1. (8) The generalised hyperbolic distribution is an example of a normal mean-variance mixture. That is, a random variable X follows a generalised hyperbolic distribution if

X= µ + βW +d √W Z Z ∼ N (0, 1)

W ∼ GIG(λ, δ2, α2− δ2_{) independent of Z}

where W follows a generalised inverse Gaussian distribution. The great advantage of the GHD is its flexibility. It’s location and scale are respectively governed through µ and δ, β in-fluences the asymmetry of the distribution (the distribution is symmetric for β = 0), α drives the kurtosis and λ distinguishes several subclasses. For instance, a GHD with λ = −1/2 is also known as the normal inverse Gaussian distribution (NIGD). Other special cases of the GHD are given in table 4.

The Student t-distribution, which is a special case of the GHD, is also a suitable candi-date, but does not meet the unit variance requirement. The standardised version of the Student t-distribution was introduced by Bollerslev [1987][10] and is characterised by

ft(y | ν) = Γ ν+1 2
p(ν − 2)πΓ ν 2
 1 + y 2 ν − 2 −ν+12 Θ = (2, ∞) Y = R.

The condition ν > 2 ensures that the variance of a non-standardised t-distribution is finite, otherwise the standardised t-distribution does not exist.

The generalised error distribution (GED) (Nadarajah [2005] [40]) can also be standardised, its standardised version is characterised by

fGED(y | κ) = κ exp  −1 2    q 2−2/κ Γ(1/κ) Γ(3/κ)y    κ q 2−2/κ Γ(1/κ) Γ(3/κ)2(1+κ)/κΓ(1/κ) Θ = (0, ∞) Y = R.

To all of the symmetric distributions mentioned above, skewness can be introduced due to Fern´andez and Steel [1998] [31] using the following method. Let f (y) be the pdf of a univariate, unimodal distribution which is symmetric around zero. Then by introducing a skewness parameter ξ ∈ (0, ∞), one can define a class of skewed distributions with density

Clearly, (9) amounts to f (y) if ξ = 1. The central moments can then derived be using M (τ ) =R∞

0 s

τ_{2f (y)dy. Specifically, the mean and the variance for the skewed distribution} corresponding to (9) are given by

E(y) = M (1)(ξ − ξ−1) (10)

Var(y) = (M (2) − M (2)2)(ξ2+ ξ−2) + 2M (1)2− M (2). (11) The normal, Student-t and GED distributions have skewed versions which can be standard-ised using moment conditions (10) and (11).

2.1.9

Comparing nested models

Nested (and unnested) models can be compared to eachother using information criteria. These criterion functions favour parsimonious models with a high log likelihood function, i.e. they penalise models with more parameters, which is useful for determining the ARMA-GARCH order. Let θ ∈ Θ be the k-dimensional estimator to be estimated by the maximum likelihood estimator denoted by ˆθMLEand let `(θ) be the likelihood function. A selection of information criteria is listed below8_.

AIC = 2k − 2 ln `(ˆθMLE) (Akaike information criterion) (12) BIC = k ln n − 2 ln `(ˆθMLE) (Bayes information criterion) (13) CAIC = k[ln(n) + 1] − 2 ln `(ˆθMLE) (Consistent AIC) (14)

AICc = AIC +2k(k + 1)

n − k − 1 (Corrected AIC) (15)

HQC = −2 ln `(ˆθMLE) + 2k ln(ln n) (Hanan-Quinn criterion) (16) where n is the total amount of observations. Criteria (12)-(16) are respectively due to Akaike [1974] [2], Schwarz [1978] [44], Bozdogan [1987] [12], Hurvich and Tsai [1989] [35] and Hanman and Quinn [1979] [33]. Note that they only differ in the penalisation of the number of independent parameters. The penalty term of the BIC is larger than that of the AIC 9 _{and increasing in sample size. The CAIC and the AICc are both modifications of} the AIC. Bozdogan [1987][12] proposed the CAIC as an improvement of the AIC to make the criterion asymptotically consistent in estimating the order of the model and to penalise overparameterisation more stringently. The AICc corrects the AIC to make it more suitable in case of a small sample size, or if the number of parameters is large relative to the amount of observations. Certainly, there exists literature on the adequacy of information criteria (e.g. Burnham and Anderson [2004] [13]). However, this falls outside the scope of this study and hence a thorough discussion is not provided here.

2.2

Econometric option pricing

In the area of option pricing one can distinguish two approaches, the continuous time ap-proach and the discrete time apap-proach. The former typically aims at finding closed form option pricing formulas based on absence of arbitrage arguments, whereas the latter ap-proach is focused on econometric tractability as empirical data is observed in discrete time. Since data is only observed in discrete time, the estimation of the parameters of continuous time models can be extremely complex. Sufficiently complex that practitioners usually re-sort to calibration instead of estimation, which entails choosing the parameters such that the model implied prices are as close to the observed prices as possible. Note that there are no quoted prices for cryptocurrency options, rendering the continuous time approach infeasible. As a consequence, this study relies on the discrete time approach.

The discrete time approach has several advantages. For one that one can rely on econo-metric techniques to find a model that fits the historical return data the best, instead of relying on assumptions as to which modelling choice is best. Moreover, one can find op-tion pricing strategies for opop-tions that are not quoted or for highly illiquid opop-tions. The theoretical background provided in this section is a description of the three key ingredients of option pricing according to Bertholon et al. [2008] [7]: the historical P-dynamics of the 8

Note that researchers might use linear transformations of the listed criteria, naturally these do not affect the model selections.

asset returns, the risk-neutral (RN) Q-dynamics of the asset returns and the stochastic dis-count factor (SDF) as well as a description of the three pricing approaches based on these ingredients. Technical details surrounding the implementation shall be deferred to section 5.

2.2.1

Historical P-dynamics

The P-dynamics are based on the observed asset returns rt, where t ∈ T = {0, . . . , T }. The information available at time t is given by r_t= {r0, r1, . . . , rt} and the σ-algebra generated by r_t is given by σ(r_{t). The historical P-dynamics are defined by the joint distribution} of r_t and is characterised by either the conditional density fP_(r

t+1 | rt), the conditional Laplace transform ϕP_{(u | r}

t) = EP[exp(urt+1)|rt] or the conditional log-Laplace transform ψP_{(u | r}

t) = log ϕP(u | rt).

2.2.2

The stochastic discount factor

Bertholon et al. [2008] [7] propose three assumptions based on Hansen and Richard [1987] [34]. These are respectively regarding the existence and uniqueness of a price for any con-tingent claim, linearity and continuity of the concon-tingent claim price and absence of arbitrage (for specific technical details see Bertholon et al. [2008] [7]). Under those three conditions one can guarantee the existence and the uniqueness of the SDF, denoted by Mt,t+1(rt+1), between t and t + 1 for all t ∈ T . Suppose the current date is t and a contingent claim pays Φs(rs) at date s > t. Then the unique price of the contract is given by

Pt[Φs(rs)] = EMt,t+1(rt+1) · · · Ms−1,s(rs)Φs(rs) | rt . The SDF is assumed to be exponential-affine, i.e. it can be written as

Mt,t+1(rt) = exp[αt(rt)rt+1+ βt(rt)].

Furthermore, the relationship between the predetermined short rate ρt+1 and the SDF is given by

E[Mt,t+1(rt+1) | rt] = exp(−ρt+1).

2.2.3

RN Q-dynamics

The historical dynamics of the asset returns rt10 are given by the joint distribution of rT under the measure P, whereas the RN dynamics are given by the joint distribution of rtunder the measure Q. Given the historical conditional density fP_(r

t+1| rt), the RN-dynamics are characterised by

fQ_(r

t+1| rt) = fP(rt+1| rt)dQ(rt+1 | rt) under the measure P, where

dQ_(r

t+1| rt) =

Mt,t+1(rt+1) E[Mt,t+1(rt+1) | rt] = exp(ρt+1)Mt,t+1(rt+1). Similarly, the historical conditional pdf is related to the RN pdf by

fP_(r

t+1| rt) = fQ(rt+1 | rt)dP(rt+1 | rt) under the measure Q where

dP_(r t+1 | rt) = 1 dQ_{(rt+1| r} t) .

Similar to the P-dynamics, it is common practice to characterise the Q-dynamics using the (log) Laplace transform.

10

Bertholon et al. [2008] [7] describe the methods more generally for a vector of risk factors wt. In the case of

2.2.4

Pricing strategies

The combination of the aforementioned mathematical tools/ingredients forms an economet-ric asset peconomet-ricing model (EAPM). Since none of these key tools are observed, they have to be estimated. Bertholon et al. [2008] [7] describe three estimation strategies: the direct modelling approach, the RN constrained direct modelling approach and the back modelling approach. The three methods differ in the way the three mathematical tools are obtained. In all of the three approaches two of the three tools are specified/parameterised and estimated after which the third ingredient follows.

Firstly, in order to form an EAPM, one needs the short rate ρt+1 which can be known, in which case it can be either exogenous or endogenous, or it can be unknown in which case it has to be parameterised. Once the short rate has been specified the EAPM can be formed. In case of the direct modelling approach one chooses a parametric family for the P-dynamics and specifies the stochastic discount factor, after which the Q-dynamics are im-plied. Hence, the Q-dynamics are obtained as a by-product. In case of the RN constrained direct modelling approach one specifies a parametric family for both the P-dynamics and the Q-dynamics and the SDF is obtained as a by-product. Finally, the back modelling approach involves specifying the SDF and the parametric family of the Q-dynamics, after which the P-dynamics are obtained.

3

Data source and description

The dataset used in the present study contains time series of five cryptocurrency prices and one cryptocurrency index. The cryptocurrencies are Bitcoin (BTC), Ethereum (ETH), Ripple (XRP), Bitcoin Cash (BCH) and EOS (EOS), of which the data is from YAHOO Finance. These currencies account for roughly 70 % of the market capitalisation of the cryptocurrency market. The cryptocurrency index, or CRIX (Trimborn and H¨ardle [2016] [46]), is an index comrprising of numerous different cryptocurrencies which can be found on the website of the Humboldt-university of Berlin. The end date of the dataset, i.e. all timeseries, is the 31st of March. The first date of each time series is different depending on the availability of the data. Bitcoin is the oldest cryptocurrency in the dataset and has data available from the sixteenth of July 2010 (2812 observations), whereas the youngest cryptocurrency - Bitcoin Cash - has data available from the first of August 2017 (244 ob-servations). The CRIX, Ethereum, Ripple and EOS time series respectively have 1340, 969, 1166 and 277 observations. The rest of this section focuses on the time series properties of for one the cryptocurrency prices (typically denoted by Pt), but especially the log re-turns on the cryptocurrencies (defined by rt= log(Pt/Pt−1)). The forthcoming analysis is largely inspired by the analysis of financial time series by McNeil et al. [2005] [39], chapter 4. Firstly, consider the time series paths of the cryptocurrencies given in figure 1 in which the turbulence in the cryptocurrency prices mentioned in section 1 is clearly visible. All of the price paths exhibit a spike in prices from the end of 2017 followed by a sharp decline around the beginning of 2018. Note that the increase and decline appear to be less prominent for Bitcoin cash and EOS, this is due to the fact that these are fairly young cryptocurrencies of which there is not even a year of data available. The empirical densities of the cryptocur-rency prices are given in figure 2. The price behaviour of the cryptocurrencies is clearly reflected in the extremely right skewed shapes of the empirical distributions. Again, less so for the younger cryptocurrencies Bitcoin cash and EOS. Table 1 provides the summary statistics of the cryptocurrency prices. The reader can clearly observe the price ranges the cryptocurrencies have attained during their life-span. Where Bitcoin has reached a maximum of almost twenty thousand USD, Ripple never got past three USD. Again the right skewed nature is reflected in the medians relative to the means and the skewness statistics. A final observation is that the excess kurtosis indicates that the price distributions are leptokurtic, i.e. heavy tailed.

Turning to the time series paths of the log returns, depicted in figure 3, once can observe trajectories that fluctuate erratically around zero. Table 2 provides the summary statistics of the log returns. Indeed, the distributions of the log returns appear to be more stable than the distributions of the prices in the sense that they all have a mean and median of approximately zero. The fact that the means and medians of the cryptocurrencies are not far apart could indicate that the distributions are rather symmetric. Though, one should not overlook the fact that the skewness statistics are all non-zero and in two cases negative. Hence, the log return distributions are in some cases left skewed and in other cases right skewed, but never completely symmetrical. Perhaps the most striking observation that can be made when examining table 2 is the high values of excess kurtosis, indicating that the log return distributions are extremely leptokurtic. This underlines the prevalence of extreme observations in the log returns, which is not surprising considering the hype in the end of 2017. Indeed, the trajectories depicted in figure 3 exhibit some extreme spikes. Especially the Bitcoin trajectory, where one daily log return reaches a value of of nearly 1.5.

the serial correlation is negligible for Bitcoin Cash and EOS. The same can be observed when comparing the ACFs of the log returns to the squared log returns, though the differences are far less pronounced than the difference with the ACFs of the absolute log returns. Again this is a stylised fact which is also established in McNeil et al. [2005] [39], and it buttresses the observation of volatility clustering. It is no surprise therefore that the ACFs of the absolute and the squared log returns of Bitcoin Cash and EOS exhibit little serial correlations, since their trajectories barely exhibit any volatility clustering.

Finally, a well-established fact about financial time series is that financial log return series are highly non-normal. Figure 7 depicts a plot of the empirical densities of the log returns, with a fitted normal distribution and a fitted location-scale t-distribution 11_{. Clearly the} normal distribution does not fit well compared to the empirical distribution. The empirical distribution has more mass in the very centre of the distribution than the normal distribu-tion, which is more evenly distributed. The location-scale t-distribution provides a better fit, which can be explained by the fact that the t-distribution accommodates leptokurtosis, whereas the normal distribution does not. The inadequacy of the normal distribution was further established by performing a Shapiro-Wilk test, a Jarque-Bera test and a Kolmogorov-Smirnov test, all of which p-values were equal to zero for all cryptocurrencies, from which can be concluded that the log return distributions are non-normal. Indeed, the QQ-plots in figure 8 indicate that the log return distributions are fairly non-normal. The ‘S-shapes’ of the quantile plots are indicative of the inability of the normal distribution to accomodate leptokurtosis, which is indeed a frequently cited inadequacy of the normal distribution when it comes to modelling financial time series.

In conclusion, in this section some statistical properties of the cryptocurrency prices and especially their log returns have been established. For one the prices of cryptocurrencies have been extreme volatile during their relatively short life span (compared to e.g. stock prices/indices). The log returns all fluctuate erratically around zero and exhibit extreme leptokurtosis and non-symmetry. Furthermore, volatility clustering has been observed in four out of the six trajectories of the log returns, the exception being Bitcoin Cash and EOS. Moreover, the distributions of the log-returns are non-normal. It is therefore imperative that the modelling choices for the log return series accommodate the aforementioned statistical properties. Note that the GARCH specifications given in subsection 2.1 accomodate non-constant volatility and volatility clustering by construction and that some of the innovation distributions provided in subsection 2.1.8 accomodate skewness, leptokurtosis or both. The next section proceeds by discussing the model selection process and the resulting estimated models.

11_{That is if a rv. X follows a location-scale t-distribution with location µ and scale σ, then X} d

4

Model specifications

4.1

Selection procedure and resulting models

In this section the model specifications for CRIX and the five cryptocurrencies described in section 3 are given. All of the models are ARMA(p,q)-GARCH(n, m) models, where the GARCH specifications and innovations distributions ((skew) normal, (skew) Student-t, (skew) GED, NIGD and GHD) are chosen from the models and distributions listed in sub-section 2.1. The model selection process involves two steps, selection of the order (p, q, n, m) and selction of the GARCH/distribution combination.

Step one involves an order selection algorithm based on the AIC (see (12)), which works as follows. Given a GARCH/distribution type, the algorithm starts with the benchmark order (p, q, n, m) = (1, 1, 1, 1). In each step, the algorithm checks if increasing an ARMA-GARCH order by one (in the order n → m → p → q) improves the AIC. If it does, the algorithm continues using the inrcreased order and if it does not, the algorithm advances to the the next step with the previous order, which is now fixed permanently. The algo-rithm stops if all orders (p, q, n, m) have been fixed. This algoalgo-rithm is applied to all the GARCH/distribution types, yielding 56 model specifications. Among these 56 models, the model with the highest log likelihood is chosen.

The reason that the AIC is the criterion used for picking the ARMA-GARCH order is simply that it is the most conventional one. This may appear arbitrary, however, there is for one no a priori reason to believe the information criteria (13)-(16) are more adequate than the AIC, nor that they will provide very different results. On the other hand performing a sensitivity analysis with the other information criteria is rather tedious and time consuming and is unlikely to provide substantial insight. Finally, the log-likelihood is used in the second step as a selction criterion, as opposed to an information criterion. The reason for this is that parsimony is already established in the first step. Hence, model adequacy is preferred over parsimony in the second step. The rest of this section proceeds by describing estimation results of the selected models for all of the cryptocurrencies and the CRIX.

CRIX

The optimal model selected for the CRIX index is an ARMA(3,3)-eGARCH(3,1) model with sGED innovations. Specifically, the estimated model is given by

ˆ rt= −0.0009 (0.0007) + 0.6020 (0.0079)rt−1+ −0.2325(0.0037) rt−2+ 0.6270 (0.0087)rt−3− 0.6368(0.0070)εˆt−1+ ˆεt + 0.2262 (0.0001)εˆt−2− 0.5655(0.0065)εˆt−3 ˆ εt= ˆσtzt, zt i.i.d. ∼ sGED(0.8600 (0.0424), 0.9554(0.0137)) ln ˆσ_t2= −0.0937 (0.0128)− 0.0610(0.0273)εˆt−1+ 0.4594(0.0867)Σt−1+ −0.0119(0.0119)εˆt−2− 0.1235(0.0800)Σt−2+ + 0.1259 (0.0297) ˆ εt−3− 0.0976 (0.0418) Σt−3+ 0.9866 (0.0020) ln ˆσ2_t−1 Σt= ˆεt− E|εt|

where the standard errors are given in parentheses. The parameters of the sGED here are respectively the shape parameter κ and the skewness parameter ξ.

Bitcoin

where the parameters of the GHD are respectively λ, α, β, δ and µ, of which there are no standard errors, since the software used to estimate the model (R) uses a different parame-terisation.

Ethereum

The optimal model selected for Ethereum is an ARMA(3,1)-CS-GARCH(1,2) model with sGED innovations. Specifically, the estimated model is given by

ˆ rt= 0.0024 (0.0002)+ 0.8307(0.0073)ˆrt−1+ 0.0881(0.0038)rˆt−2+ 0.0248(0.0017)rˆt−3+ ˆεt− 0.9192(0.0085)εˆt−1 ˆ εt= ˆσtzt, zt∼ sGED(0.9362 (0.0433), 1.067(0.0140)) ˆ σ_t2= ˆqt+ 0.1952 (0.000) (ˆε2_t−1− ˆqt−1) + 0.5273 (0.0289) ˆ σ_t−12 + 0.2046 (0.0364) ˆ σ_t−22 ˆ qt= 0.0001 (0.000)+ 0.9994(0.0004)qˆt−1+ 0.1172(0.0046)(ˆε 2 t−1− ˆσ2t−1).

Ripple

The optimal model for Ripple is an ARMA(1,1)-CS-GARCH(1,1) with sGED innovations. Specifically, the estimated model is given by

ˆ rt= −0.0015 (0.0000) + 0.1998 (0.0041) ˆ rt−1+ ˆεt− 0.4541 (0.0090) ˆ εt−1 ˆ εt= ˆσtzt, zt∼ sGED(0.8082 (0.019) , 1.0881 (0.0104) ) ˆ σ2t = ˆqt+ 0.1999 (0.0148)(ˆε 2 t−1− ˆqt−1) + 0.2578 (0.0234)(ˆσ 2 t−1− ˆqt−1) ˆ qt= 0.0000 (0.0000)+ 0.9999(0.0003)qˆt−1+ 0.04568(0.0017)(ˆε 2 t−1− ˆσ 2 t−1).

Bitcoin Cash

The optimal model for Bitcoin Cash is an ARMA(2,1)-APARCH(3,1) model with sGED innovations. Specifically, the estimated model is given by

ˆ rt= −0.0020 (0.0052)+ 0.4975(0.4992)ˆrt−1− 0.0771(0.0556)rˆt−2− 0.4930(0.4937)εˆt+ ˆεt−1 ˆ εt= ˆσtzt, zt∼ NIG(0.7820, 0.1321, 0.7488, −0.1283) ˆ σ 3.5000 (1.7264) t = 0.0073 (0.0183)(|ˆεt−1| − 0.6964(1.1615)εˆt−1) 3.5000 (1.7264)_{+ 0.0000} (0.0000)(|ˆεt−2| − 0.4114(1.3893)εˆt−2) 3.5000 (1.7264) + 0.0022 (0.0052)(|ˆεt−2| − 0.9985(0.0039)εˆt−2) 3.5000 (1.7264)_{+ 0.8001} (0.1774)σˆ 3.5000 (1.7264) t−1

where the parameters of the NIGD are specifically α, β, δ and µ, of which the standard errors cannot be provided for the same reason as for the GHD.

EOS

The optimal model for EOS is an ARMA(2,1)-gjrGARCH(2,2) model with sGED innovations ˆ rt= −0.0031 (0.0000) + 0.9137 (0.0034)rˆt−1+ 0.0505(0.0000)ˆrt−2+ ˆεt− 0.9406(0.0054)εˆt−1 ˆ εt= ˆσtzt, zt∼ sGED(0.9251 (0.0865) , 1.1236 (0.0279) ) ˆ σt2= = 0.0012 (0.000)+ 0.1558(0.0007)εˆ 2 t−1− 0.3057 (0.0012)1{ˆεt−1< 0}ˆε 2 t−1+ 0.0009 (0.0000)εˆ 2 t−2 + 0.5417 (0.0027)1{ˆεt−2< 0}ˆε 2 t−2+ 0.1927 (0.0009)ˆσ 2 t−1+ 0.4545 (0.0021)ˆσ 2 t−2

4.2

Discussion and concluding remarks

5

Option pricing

The aim of this section is to outline an option pricing methodology using the ingredients presented in section 2. The specific aim is to price plain vanilla European call options for different levels of moneyness and maturities, though the methods presented here can easily be adjusted to price different contingent claims.

5.1

Methodology

The pricing methodology applied in this study is the methodology used by Chorro et al. [2012] [19] in case of GHD or NIGD innovations and the methodology proposed by Elliott and Madan [1998] [27], the latter of which is based on what one could consider as the Girsanov Theorem applied extended to a discrete time framework. As is mentioned in subsection 2.2, the methodology is a direct modelling approach, hence the first order of business is to specify the P-dynamics and the SDF. The specification of the P-dynamics depends on the outcome of the algorithm described in section 4 and the SDF shall be exponentially affine as in subsection 2.2.2. The rest of this subsection proceeds by describing how to obtain the Q-dynamics in case of the GHD and the NIGD and in case of the sGED.

5.1.1

Obtaining the Q-dynamics: the GHD and NIGD case

In the general case a discrete time economy is considered on the time axis T = {0, . . . , T } where T ∈ N++ ≡ {k ∈ N : k > 0}. The economy consists of a cryptocurrency priced at Pt at time t ∈ T and a risk-free zero coupon bond priced at Bt. The log-return on the cryptocurrency defined on the filtered probability space (Ω, F , P) with F0= {∅, Ω} and Ft= σ(rt) for t ∈ T , is given by

rt≡ log Pt Pt−1

= mt+ σtzt, t ∈ T \{0} (17)

where mt is the deterministic part of (1) conditional on Ft−1, σt arises from a GARCH specification as in subsection 2.1 and zt is a zero-mean, unit-variance random innovation distributed according to one of the specifications in subsection 2.1.8. The bank account fol-lows the dynamics Bt= exp(ρ)Bt−1 for t ∈ T \{0} and B0= 1, where ρ is the risk-free rate expressed on a daily basis. Finally, markets are assumed to be frictionless which encompasses 1. no transaction costs 2. borrowing at the risk free rate and short-selling the cryptocurrency is allowed and 3. the cryptocurrency is perfectly divisible.

Key to obtaining the Q-dynamics is the notion of absence of arbitrage opportunities (AAO). Key to obtaining any reasonable option pricing method is that the market is free of arbitrage 12

. The measure Q is chosen such that Q is an equivalent martingale measure (EMM), i.e. Q is equivalent to P and we have that ∀t ∈ T \{0} that EQt−1[Pt] = exp(ρ)Pt−1, where the subscript t − 1 of the expectation operator indicates that the expectation is conditional on Ft−113. One can prove that the existence of such an EMM and the notion of AAO are equiv-alent, known as the Fundamental Theorem of Asset Pricing, see e.g. Dalang et al. [1990] [22]. Contrary, to the continuous case, this EMM is generally not unique in the discrete case. Consequently, a method of specifying the EMM is required. Chorro et al. [2012] [19] do this by using an exponentially affine stochastic discount factor given by

Mt,t+1= exp(θt+1rt+1+ ξt+1) (18)

where θt+1 and ξt+1 need to be computed explicitly. Note that using the EMM, the AAO price Xtat time t of a European contingent claim paying Ψτ at time τ > t is given by

Xt= exp(−ρ(τ − t))EQt[Ψτ] = EPt[ΨτMt,τ] (19) where Mt,τ =Q

τ −1

k=1Mk,k+1. Now we let τ = t + 1 so the internal consistency constraints that follow from (19) are given by

EPt[exp(ρ)Mt,t+1] = 1 EPt[exp(rt+1)Mt,t+1] = 1.

(20) 12

In this study it is assumed that the reader is familiar with the notion of arbitrage and basic mathematical option pricing. If this is not the case one can consult any educative book on mathematical option pricing or mathematical finance.

13_{Note that in subsection 2.2 I typically condition on r}

t−1. In the GARCH setting used in this study Ft−1and

Let GP

rt+1|Ft : D

P

G→ R++be the conditional moment generating function (conditional MGF) of rt+1conditional on Ftgiven by

GPrt+1|Ft(u) = E

P

t[exp(urt+1)] (21)

where DP

G is the domain which depends on the distribution of rt+1 under the measure P. Now let Φt: ΘGt → R be defined by

Φt(u) = log GPrt+1|Ft(1 + u) GPrt+1|Ft(u) ! (22) where ΘG

t = {u ∈ R : u, 1 + u ∈ DP_G}. Using the notation introduced in (21) and (22) in combination with (20) gives

Φt(θt+1) = ρ

GPrt+1|Ft(θt+1+ 1) = exp(ξt+1).

(23)

Again, when one prices a European contingent claim at time t, paying Ψτ at time τ > t, one can do so by using the discounted expected payoff under the Q-dynamics as in (19). Once the parameters {θk}τk=t+1and {ξ}τk=t+1have been calculated explicitly, one can calculate the expected payoff under the Q-dynamics by shifting the P-dynamics using the Radon-Nikodym derivative given by dQ dP = exp(ρ(τ − t))Mt,τ = τ Y k=t exp(θkrk) GPrk|Fk−1(θk) (24) which is also known as the conditional Esscher transform. Chorro et al. [2012] [19] specifi-cally apply the methodology described above to the case where the innovations are from the GHD. Note that the fact that the GHD with λ = −1

2 is the NIGD, which permits us to use the methodology in case of the NIGD as well.

Assume the same economic settings as in the previous subsection and let ztin (17) follow a GHD. The moment generating function of ztis given by

Gzt(u) = exp(uµ)  _α2_{− β}2 α2_{− (β + u)}2 λ2 _K λ(δpα2− (β + u)2) Kλ(δ p α2_{− β}2₎ , |β + u| < α.

Since the GHD is closed under linear transformations (Barndorff-Nielsen and Blaesild [1981] [6]) we have that the conditional distribution of rtgiven Ft−1under the P-dynamics is

GH  λ, α σt , β σt , δσt, mt+ µσt  .

Consequently, the conditional MGF of rtgiven Ft−1under the P-measure is given by

GPrt|Ft−1(u) = exp(u(mt+ µσt))  _α2_{− β}2 α − (β + uσt)2 λ2 _K λ(δpα − (β + uσt)2) Kλ(δ p α2_{− β}2₎ where |β + σtu| < α. Finally, the following proposition provides the Q-dynamics of rt

Proposition 5.1 (Chorro et al. (2012) [19]). The conditional distribution of rtgiven Ft−1 is given by GH  λ, α σt , β σt + θt, δσt, mt+ µσt  . (25)

Equivalently the Q-dynamics of rtgiven Ft−1are given by

rt= mt+ σtzt, zt∼ GH(λ, α, β + σtθt, δ, µ).

It is interesting to note that the variable ξt does not appear in the specification of the Q-dynamics here, nor in the Radon-Nikodym derivative (24) from the previous subsection. Consequently, the only equation from (23) that needs to be solved is the first one. Finally, the domain of Φtin (23) is given by ΘGt =

The pricing algorithm that unfolds is as follows. Consider the sample of prices and log returns of a cryptocurrency, where the returns are modelled using an ARMA-GARCH spec-ification with GHD or NIGD innovations, given by respectively {Pt}t∈T and {rt}t∈T \{0}, where T is as before. The first date t = 0 shall differ per cryptocurrency, the final date T is always the 31st of March 2018. Let PT be the last cryptocurrency price of the sample, let τ be the time to maturity of the call option contract expressed in days, let M be the moneyness of the contract such that the strike price is given by K = M · PT. Then the price of the European call option at time T , given by CT = exp(−ρτ )EQT(max{PT +τ− K, 0}), is estimated as follows.

Algorithm 1

1. The ARMA-GARCH model is estimated on the log returns over the sample T . The model yields estimates for {σt, εt}_{t∈ ˜T}, given by {ˆσt, ˆεt}_{t∈ ˜T}, where ˜T = {k, . . . , T } ⊂ T , k ∈ N depends on the order of the ARMA-GARCH model.

2. Use {ˆσt, ˆεt}_{t∈ ˜T} to obtain ˆσT +1using the estimated GARCH model. (a) Start with the estimated value of ˆσT +1.

(b) Solve the first equation of 23 to obtain θT +1.

(c) Sample the next innovation zT +1from a GH(λ, α, β + σT +1θT +1, δ, µ) distribution. (d) Compute the log return ˆrT +1using the estimated ARMA equation and compute

ˆ

σT +2 using the estimated GARCH equation.

(e) Return to (a) and replace T with T +1 and repeat untill the sample { ˆrt}t∈{T +1,...,τ +T } is obtained.

3. The result of the simulation in the previous step will give a sample path of the cryp-tocurrency returns (under the Q measure) from which the crypcryp-tocurrency price at time T + τ is obtained given by ˆ PT +τ = PT T +τ Y t=T +1 exp(ˆrt).

4. Steps 1. and 2. are repeated until the sample { ˆPi

T +τ}i∈{1,...,N }, where N = 1000, is obtained. Finally, the estimated call price of a European call option with strike price K = M · Pt and time to maturity τ is estimated by

ˆ CT = exp(−ρτ ) 1 N N X i=1 max{ ˆPT +τi − K, 0}.

5.1.2

Obtaining the Q-dynamics: the sGED case

The methodology of the previous subsection requires that the moment generating function (21) exists. However, if the distribution of the innovations ztin (17) is the sGED, the moment generating function does not exist. Hence, in this case we rely on the methodology proposed by Elliott and Madan [1998] [27], which is based on a static Girsanov-type transformation that shifts the conditional mean of the returns. The advantage here is that the algorithm used is quicker than the one used in the previous subsection. The caveat is that this method only shifts the conditional mean whereas θt in the method by Chorro et al. [2012] [19] also shifts the contitional variance, kurtosis and skewness. This has been cited as a reason that this method has a poor pricing performance for longer maturities (Badescu and Kulperger [2008] [3]).

The methodology starts with a multiplicative Doob decomposition of the discount cryp-tocurrency price { ˜Pt}t∈T into a predictable finite variation component and a martingale component

˜

Pt= ˜P0AtNt

deduce that ˜ Pt ˜ Pt−1 = AtNt At−1Nt−1 ⇐⇒ E P t−1 ˜ Pt ˜ Pt−1 ! = At At−1 ⇐⇒ At= EPt−1 ˜ Pt ˜ Pt−1 ! At−1 ⇒ At= t Y i=1 EPi−1 ˜ Pi ˜ Pi−1 ! Consequently, ˜ Pt ˜ Pt−1 = Qt i=1EPi−1  _˜ Pi ˜ Pi−1  Qt−1 i=1EPi−1  _˜ Pi ˜ Pi−1  Nt Nt−1 = EP t−1 ˜ Pt ˜ Pt−1 ! Nt Nt−1 = exp(−ρ)EP t−1(exp(rt)) ˜Pt−1 Nt Nt−1 = exp(vt) ˜Pt−1 Nt Nt−1

where exp(vt) = exp(−ρ)EPt−1(exp(rt)). Now let ftP be the density of Nt/Nt−1 given Ft−1 under P. Now let the process Z˜t,tbe defined by

Zt,t˜ = t Y k=˜t fP k  _˜ Pk ˜ Pk−1  exp(vt) fP k  _˜ Pk ˜ Pk−1 exp(−vt)  . (26)

Then, a contingent claim maturing at T + τ can be priced at time T by using the Radon-Nikodym derivative given by

dQ

dP = ZT +1,T +τ. (27)

In the ARMA-GARCH setting the processes ˜Pt/ ˜Pt−1 and Nt/Nt−1are simply given by re-spectively exp(rt− ρ) and exp(rt− ρ − vt).

Note that in order to obtain the Radon-nikodym derivative (27) the density fP

t is required, as well as the conditional density of exp(rt) given Ft−1under P, say gtP. Firstly recall that exp(rt) = exp(mt+ σtzt) where mtand rtare Ft−1predictable, which entails the only ran-dom part in the expression of exp(rt) is zt∼ sGED(κ, ξ). The density gPt can be obtained by a simple transformation technique. Define yt: R → R++ by yt(zt) = exp(mt+ σtzt). Since ytis invertible, zt in terms of yt is given by zt(yt) = (log yt− mt)/σt. Now let fsGED(·|κ, ξ) be the pdf of the sGED, so the density gP

t is given by gP

t(y) = fsGED(y|κ, ξ)     dzt(y) dy     = fsGED(y|κ, ξ) 1 ytσt where the absolute bars can be omitted since the support of gP

t is strictly positive. The density fP

t can be obtained analogously by using the function yt= exp(mt+ σtzt− ρ − vt). The algorithm that unfolds here is as follows. Similarly to the previous subsection con-sider the sample of cryptocurrency prices and returns given by {Pt}t∈T and {rt}t∈T \{0}, respectively. Again, the aim is to price a European call option at time T with time to ma-turity τ , again expressed in days. Like in the previous subsection, the moneyness given by M induces the strike price K = M · PT. The price of the European call option is as follows.

Algorithm 2

1. The ARMA-GARCH model is estimated on the log returns over the sample T . The model yields estimates for {σt, εt}t∈ ˜T, given by {ˆσt, ˆεt}t∈ ˜T, where ˜T = {k, . . . , T } ⊂ T , k ∈ N depends on the order of the ARMA-GARCH model.

3. {exp(ˆvt,i)}t∈{T +1,...,T +τ },i∈{1,...,N } are computed by calculating exp(ˆvi,t) = exp(−ρ)

Z ∞

0 ygP

t(y)dy.

4. The final estimated stock prices (using the P-dynamics) { ˆPi

T +τ}i∈{1,...,N }are given by ˆ P_{T +τ}i = PT ,i T +τ Y t=T +1 exp(ˆrt,i).

The estimated Radon-Nikodym derivatives { ˆZ_{T +1,T +τ}i } (cf. (26) and (27)) are given by ˆ Z_{T +1,T +τ}i = T +τ Y t=T +1 fP k(exp[ ˆrt− ρ]) exp[ˆvi] fP k(exp[ ˆrt− ρ] exp[−ˆvi]) .

5. Finally the estimated European call price is given by ˆ CT = exp(−ρτ ) 1 N N X i=1 max{ ˆP_{T +τ}i − K, 0} ˆZT +1,T +τ

5.2

Implementation and results

In this section the results from the algorithms proposed in subsection 5.1 are presented. Algo-rithm 1 is applicable to the ARMA-GARCH specifications with GHD or NIGD innovations, i.e. Bitcoin and Bitcoin Cash. Algorithm 2 is applicable to the ARMA-GARCH specifica-tions with sGED innovaspecifica-tions, i.e. CRIX, Ethereum, Ripple and EOS. The simulated sample paths used to estimate the European call prices start from the first of April 2018, which lies out-of-sample. As a consequence there is no accurate data (or estimate) for the daily continuously compounded risk free rate ρ. However, since there is no data on cryptocurrency option prices and consequently no data to verify the validity of the resulting call prices, there is also no need for an accurate risk free rate. Hence, the daily continuously compounded risk free rate will be arbitrarily set to ρ = 0.001. Admittedly, ρ can be set to any arbitrary value, as long as it is not too high as this would yield unrealistically low option prices. The sets of moneyness and maturities for which the algorithms are applied are respectively given by ¯M = {0.8, 0.85, 0.9, 0.95, 1.0, 1.05, 1.1, 1.15, 1.2} and ¯τ = {30, 90, 180, 360, 720}, i.e. the maturities range from a month to two years. Finally, the P-dynamics are estimated using all the available data. The reason for this is that it is common practice to use larger datasets. Chorro et al. [2012] [19] and Badescu and Kulperger [2008] [3] use respectively 4000 observations and 4290 observations, whereas the largest cryptocurrency timeseries here - Bitcoin - has 2811 observations.

That the expected call prices of Bitcoin Cash are not ‘well-behaved’ like the expected prices of Bitcoin is not surprising. In the dataset there are 2812 observations available to calibrate the P-dynamics of Bitcoin, whereas there are only 244 for Bitcoin Cash. It is therefore very likely that the P-dynamics are highly unreliable and naturally for larger maturities the sim-ulated price paths are longer than the time series from which the P-dynamics are calibrated. Moreover, note that Chorro et al. [2012] [19] prove in their paper that for a GHD(λ, α, β, δ, γ) with α > 1/2, (23) has a solution. However, the distribution in proposition 5.1 given by (25) need not exist, which was a frequently encountered problem when applying algorithm 1. Since the algorithm samples from a GHD(λ, α, β + θtσ, δ, µ), several iterations had to be discarded when α < |β + θtσ|, which falls outside of the support of the GHD. This could also partially explain the incoherent results for the Bitcoin Cash call prices. The Black-Scholes-Merton call prices for Bitcoin Cash are priced rather optimistically on the other hand, a 2-year at-the-money call option costs almost the same amount as the cryptocurrency itself. In conclusion, it is safe to say that the algorithmically expected call prices for Bitcoin Cash can be safely disregarded.

The second method used is the method by proposed by Elliot and Madan [1998] [27] us-ing a discrete Girsanov-type equivalent Martus-ingale measure which is found usus-ing algorithm 2. The prices generated by this algorithm for respectively CRIX, Ethereum and ripple are displayed in tables 11, 13 and 14. Their Black-Scholes-Merton counterparts are displayed in respectively table 16, 18 and 19. Recall from subsection 5.1.2 that this method tends to pro-duce poor pricing results for longer maturities. Indeed, when examining the prices generated by algorithm 2, the most striking observation is how quickly the prices diverge for higher maturities. All of the entries displaying “∼ ∞” represent numbers exceeding a million, often by a factor of thousands or more. The prices generated for the CRIX call option (table 11) appear to be within reason for maturities of 30 days and 90 days, compared the current price Pt= 19662.32 USD. The standard deviations of the call prices are rather high though, and start diverging extremely fast for larger maturities, which is not in favour of the reliability of the generated call prices. When comparing the expected call options prices with their Black-Scholes-Merton counterparts, where we restrict the comparison to maturities {30, 90}, it can be seen that the Black-Scholes-Merton prices are higher for 1-month call options with a low level of moneyness and vice-versa for high moneyness options. This is not the case anymore for 90-day options, as the Black-Scholes-Merton price is only higher for a money-ness of 80 %. Consequently, comparing the algorithmically generated CRIX option pricing to their Black-Scholes-Merton counterparts does not lead to consistent conclusions. Finally, the generated option prices for Ethereum Ripple, respectively given in tables 14 and 13, are unsatisfactory in a similar way.

6

Concluding remarks

In summary, in this study an econometrically tractable direct modelling approach is pro-posed in a first attempt to find a useful option pricing strategy for pricing cryptocurrency derivatives. Specifically, a direct modelling approach is proposed where the risk-neutral Q-dynamics are a by-product of the historical dynamics P-dynamics and the exponentially affine stochastic discount factor. The historical dynamics follow an ARMA-GARCH type process, where the specification is selected in two steps. In the first step the ARMA-GARCH orders are determined using the AIC and in the second step the GARCH type and the inno-vation distribution is selected using the log-likelihood. From the estimated P-dynamics, two methods are used to obtain the Q-dynamics. The first method is based on the conditional Esscher transform and the second method is based on a static Girsanov-type EMM. The former method only produced sensible results for Bitcoin whereas the latter methods failed to produce useful prices at all.

7

References

[1] Angel, J.J. and McCabe, D., 2015. The ethics of payments: Paper, plastic, or Bitcoin?. Journal of Business Ethics, 132(3), pp.603-611.

[2] Akaike, H., 1974. A new look at the statistical model identification. IEEE transactions on automatic control, 19(6), pp.716-723.

[3] Badescu, A.M. and Kulperger, R.J., 2008. GARCH option pricing: A semiparametric approach. Insurance: Mathematics and Economics, 43(1), pp.69-84.

[4] Barndorff-Nielsen, O. and Halgreen, C., 1977. Infinite divisibility of the hyperbolic and generalized inverse Gaussian distributions. Zeitschrift f¨ur Wahrscheinlichkeitstheorie und verwandte Gebiete, 38(4), pp.309-311.

[5] Barndorff-Nielsen, O., 1978. Hyperbolic distributions and distributions on hyperbolae. Scandinavian Journal of statistics, pp.151-157.

[6] Barndorff-Nielsen, O. and Blaesild, P., 1981. Hyperbolic distributions and ramifications: Contributions to theory and application. In Statistical distributions in scientific work (pp. 19-44). Springer, Dordrecht.

[7] Bertholon, H., Monfort, A. and Pegoraro, F., 2008. Econometric asset pricing modelling. Journal of Financial Econometrics, 6(4), pp.407-458.

[8] Blæsild, P., 1999. Generalized hyperbolic and generalized inverse Gaussian distributions. Working Paper, University of ˚Arhus.

[9] Bollerslev, T., 1986. Generalized autoregressive conditional heteroskedasticity. Journal of econometrics, 31(3), pp.307-327.

[10] Bollerslev, T., 1987. A conditionally heteroskedastic time series model for speculative prices and rates of return. The review of economics and statistics, pp.542-547.

[11] Bollerslev, T., 2009. Glossary to ARCH (GARCH. In Volatility and Time Series Econo-metrics: Essays in Honour of Robert F. Engle.

[12] Bozdogan, H., 1987. Model selection and Akaike’s information criterion (AIC): The general theory and its analytical extensions. Psychometrika, 52(3), pp.345-370.

[13] Burnham, K.P. and Anderson, D.R., 2004. Multimodel inference: understanding AIC and BIC in model selection. Sociological methods & research, 33(2), pp.261-304. [14] Catania, L., Grassi, S. and Ravazzolo, F., 2018. Predicting the Volatility of

Cryptocur-rency Time Series.

[15] Cheah, E.T. and Fry, J., 2015. Speculative bubbles in Bitcoin markets? An empirical investigation into the fundamental value of Bitcoin. Economics Letters, 130, pp.32-36. [16] Chaum, D., 1983. Blind signatures for untraceable payments. In Advances in cryptology

(pp. 199-203). Springer, Boston, MA.

[17] Chen, S., Chen, C.Y.H., H¨ardle, W.K., Lee, T.M. and Ong, B., 2017. Econometric Analysis of a Cryptocurrency Index for Portfolio Investment. In Handbook of Blockchain, Digital Finance, and Inclusion, Volume 1 (pp. 175-206).

[18] Chen, C.Y.H., H¨ardle, W.K., Hou, A.J. and Wang, W., 2018. Pricing Cryptocurrency options: the case of CRIX and Bitcoin.

[19] Chorro, C., Gu´egan, D. and Ielpo, F., 2012. Option pricing for GARCH-type models with generalized hyperbolic innovations. Quantitative Finance, 12(7), pp.1079-1094. [20] Chu, J., Chan, S., Nadarajah, S. and Osterrieder, J., 2017. GARCH modelling of

cryp-tocurrencies. Journal of Risk and Financial Management, 10(4), p.17.

[21] Ciaian, P., Rajcaniova, M. and Kancs, D.A., 2016. The economics of BitCoin price formation. Applied Economics, 48(19), pp.1799-1815.

[22] Dalang, R.C., Morton, A. and Willinger, W., 1990. Equivalent martingale measures and no-arbitrage in stochastic securities market models. Stochastics: An International Journal of Probability and Stochastic Processes, 29(2), pp.185-201.

[23] Ding Z, Granger CW, Engle RF. A long memory property of stock market returns and a new model. Journal of empirical finance. 1993 Jun 1;1(1):83-106.

[24] Dion, D.A., 2013. I’ll gladly trade you two bits on Tuesday for a byte today: Bitcoin, regulating fraud in the e-conomy of Hacker-cash. U. Ill. JL Tech. & Pol’y, p.165. [25] Dwyer, G.P., 2015. The economics of Bitcoin and similar private digital currencies.

[26] Dyhrberg, Anne Haubo. ”Bitcoin, gold and the dollar?A GARCH volatility analysis.” Finance Research Letters 16 (2016): 85-92.

[27] Elliott, R.J. and Madan, D.B., 1998. A discrete time equivalent martingale measure. Mathematical finance, 8(2), pp.127-152.

[28] Engle, R.F., 1982. Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica: Journal of the Econometric Society, pp.987-1007.

[29] Engle, R.F. and Bollerslev, T., 1986. Modelling the persistence of conditional variances. Econometric reviews, 5(1), pp.1-50.

[30] Engle, R.F. and Lee, G., 1999. A long-run and short-run component model of stock return volatility. Cointegration, Causality, and Forecasting: A Festschrift in Honour of Clive WJ Granger, pp.475-497.

[31] Fern´andez, C. and Steel, M.F., 1998. On Bayesian modeling of fat tails and skewness. Journal of the American Statistical Association, 93(441), pp.359-371.

[32] Glosten, L.R., Jagannathan, R. and Runkle, D.E., 1993. On the relation between the expected value and the volatility of the nominal excess return on stocks. The journal of finance, 48(5), pp.1779-1801.

[33] Hannan, E.J. and Quinn, B.G., 1979. The determination of the order of an autoregres-sion. Journal of the Royal Statistical Society. Series B (Methodological), pp.190-195. [34] Hansen, L.P. and Richard, S.F., 1987. The role of conditioning information in deducing

testable restrictions implied by dynamic asset pricing models. Econometrica: Journal of the Econometric Society, pp.587-613.

[35] Hurvich, C.M. and Tsai, C.L., 1989. Regression and time series model selection in small samples. Biometrika, 76(2), pp.297-307.

[36] Ly, M.K.M., 2013. Coining Bitcoin’s Legal-Bits: Examining the Regulatory Framework for Bitcoin and Virtual Currencies. Harv. JL & Tech., 27, p.587.

[37] Madan, D.B., Carr, P.P. and Chang, E.C., 1998. The variance gamma process and option pricing. Review of Finance, 2(1), pp.79-105.

[38] Mandelbrot, B.B., 1963. The variation of certain speculative prices. Journal of Business 36,pp. 394-419.

[39] McNeil, A.J., Frey, R. and Embrechts, P., 2005. Quantitative risk management: Con-cepts, techniques and tools (Vol. 3). Princeton: Princeton university press.

[40] Nadarajah, S., 2005. A generalized normal distribution. Journal of Applied Statistics, 32(7), pp.685-694.

[41] Nakamoto, S., 2008. Bitcoin: A peer-to-peer electronic cash system.

[42] Nelson, D.B., 1991. Conditional heteroskedasticity in asset returns: A new approach. Econometrica: Journal of the Econometric Society, pp.347-370.

[43] Nelson, D.B. and Cao, C.Q., 1992. Inequality constraints in the univariate GARCH model. Journal of Business & Economic Statistics, 10(2), pp.229-235.

[44] Schwarz, G., 1978. Estimating the dimension of a model. The annals of statistics, 6(2), pp.461-464.

[45] Taylor, S., 1986, Modelling financial time series, New York, John Wiley & Sons. [46] Trimborn, S. and H¨ardle, W.K., 2016. CRIX an Index for blockchain based Currencies. [47] Zakoian, J.M., 1994. Threshold heteroskedastic models. Journal of Economic Dynamics

A

Appendix: Tables and figures

mean

median

sd

min

max

ex. kurtosis

skewness

observations

crix

5762.84

1010.64

11120.33

342.07

62895.26

10.16

2.73

1340.00

btc

1058.67

254.68

2609.81

0.05

19345.49

19.16

3.90

2812.00

eth

167.45

12.39

277.62

0.42

1385.02

6.57

2.02

969.00

xrp

0.15

0.01

0.37

0.00

2.78

20.12

3.82

1166.00

bch

1102.74

1009.49

720.00

220.30

3715.91

3.37

0.96

244.00

eos

4.42

1.92

4.03

0.49

16.00

2.57

0.85

277.00

Table 1: Summary statistics of the cryptocurrency prices (USD).

mean

median

sd

min

max

ex. kurtosis

skewness

observations

crix

0.00

0.00

0.04

-0.25

0.20

10.04

-0.76

1339.00

btc

0.00

0.00

0.07

-0.85

1.47

89.75

2.90

2811.00

eth

0.00

0.00

0.08

-0.92

0.38

20.84

-1.21

968.00

xrp

0.00

-0.00

0.12

-1.00

1.03

22.36

0.75

1165.00

bch

0.00

-0.01

0.11

-0.47

0.44

6.55

0.43

243.00

eos

0.01

-0.01

0.12

-0.35

1.05

21.53

2.39

276.00

Table 2: Summary statistics of the cryptocurrency log returns.

mean

median

sd

min

max

kurtosis

skewness

observations

crix

-0.00

-0.00

0.04

-0.45

0.50

18.66

-0.05

5000.00

btc

0.01

0.00

0.09

-1.07

3.55

512.33

13.68

5000.00

eth

0.00

-0.00

0.24

-6.86

3.09

187.48

-3.59

5000.00

xrp

-0.00

-0.00

0.08

-1.15

0.69

25.75

-0.52

5000.00

eos

-0.00

-0.01

0.10

-0.64

1.24

10.07

0.82

5000.00

Parameter restrictions

Special case

Reference

λ = ν/2, α = β = 0, δ =

√

ν

Student t-distribution

Barndorff-Nielsen [1978] [5]

α, δ → ∞, α/β < ∞

Normal distribution

Barndorff-Nielsen [1977] [4]

λ = 1

Hyperbolic distribution

Barndorff-Nielsen [1978] [5]

λ > 0, δ → 0

Variance Gamma

Madan et al. [1998] [37]

α → ∞, δ → 0, αδ

2

_{< ∞ λ = 1, δ → 0}

_{Skewed/shifted Laplace distribution}

_{Blaesild [1990] [8]}

λ = −1/2, α = |β| → 0

Scaled/shifted Cauchy distribution

Blaesild [1999] [8]

λ = −δ

2

/2, α = |β| → 0

Shifted Student t-distribution

Blaesild [1999] [8]

2015 2016 2017 2018 0 20000 40000 60000 crix time pr ice 2012 2014 2016 2018 0 5000 10000 15000 20000 btc time pr ice 2016 2017 2018 0 200 600 1000 1400 eth time pr ice 2015 2016 2017 2018 0.0 0.5 1.0 1.5 2.0 2.5 xrp time pr ice

Sep Nov Jan Mar

500 1500 2500 3500 bch time pr ice

Jul Sep Nov Jan Mar

0 5 10 15 eos time pr ice

0 10000 20000 30000 40000 50000 60000 0e+00 2e−04 4e−04 crix pr ice density 0 5000 10000 15000 20000 0.0000 0.0005 0.0010 0.0015 btc pr ice density 0 500 1000 1500 0.000 0.002 0.004 eth pr ice density 0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 xrp pr ice density 0 1000 2000 3000 4000 0e+00 2e−04 4e−04 6e−04 bch pr ice density 0 5 10 15 20 0.00 0.05 0.10 0.15 eos pr ice density

2015 2016 2017 2018 −0.2 −0.1 0.0 0.1 0.2 crix time Log retur n 2012 2014 2016 2018 −0.5 0.0 0.5 1.0 1.5 btc time Log retur n 2016 2017 2018 −0.8 −0.4 0.0 0.2 0.4 eth time Log retur n 2015 2016 2017 2018 −1.0 −0.5 0.0 0.5 1.0 xrp time Log retur n

Sep Nov Jan Mar

−0.4 −0.2 0.0 0.2 0.4 bch time Log retur n

Jul Sep Nov Jan Mar

−0.4 0.0 0.4 0.8 eos time Log retur n