4.1. Finite mixture model of normal distributions

Assume that, for each $i = 1, 2, \ldots, N$ , the ith-component CDF $F_{i} (x)$ is the CDF of a normal distribution with mean $\mu_i$ and variance $\sigma_i^2$ (i.e. $N (\mu_i, \sigma_i^2)$ ). Then

(4.1) \begin{equation} M_i (\eta) = \exp \bigg ( \mu_i \eta + \frac{1}{2} \eta^2 \sigma_i^2 \bigg).\end{equation}

From (3.14) and (4.1),

(4.2) \begin{equation} h^{\dagger} (i) = \frac{r - \mu_i - \frac{1}{2} \sigma^2_i}{\sigma^2_i}, \quad i = 1, 2, \ldots, N.\end{equation}

Then (4.2) ensures the existence of an equivalent martingale measure $\mathbb{P}^{h^{\dagger}}$ .

From (3.16), (4.1), and (4.2),

(4.3) \begin{equation} M^{h^{\dagger}}_{i} (\eta) = \exp \bigg ( \bigg (r - \frac{1}{2} \sigma^2_i \bigg ) \eta + \frac{1}{2} \eta^2 \sigma^2_i \bigg ).\end{equation}

This implies that under $\mathbb{P}^{h^{\dagger}}$ , the conditional distribution of $R_t$ given $\Theta = i$ is a normal distribution with mean $r - \frac{1}{2} \sigma^2_i$ and variance $\sigma^2_i$ , (i.e. $R_t \mid \{ \Theta = i \} \sim N(r - \frac{1}{2} \sigma^2_i, \sigma^2_i)$ ), for each $i = 1, 2, \ldots, N$ and each $t = 1, 2, \ldots, T$ . Define

(4.4) \begin{equation} R_{1, T} \,:\!=\, \sum^{T}_{t = 1} R_t.\end{equation}

Note that by Lemma 3.2, $R_1, R_2, \ldots, R_T \mid \Theta$ are (conditionally) independent under $\mathbb{P}^{h^{\dagger}}$ . Furthermore, for each $t = 1, 2, \ldots, T$ and each $i = 1, 2, \ldots, N$ , $R_t \mid \{ \Theta = i \} \sim N\big(r - \frac{1}{2} \sigma^2_i, \sigma^2_i\big)$ under $\mathbb{P}^{h^{\dagger}}$ . Consequently, $R_{1, T} \mid \{ \Theta = i \} \sim N \big(\big(r - \frac{1}{2} \sigma^2_i\big) T, \sigma^2_i T\big)$ under $\mathbb{P}^{h^{\dagger}}$ . It may be seen that, for each $i = 1, 2, \ldots, N$ ,

(4.5) \begin{equation} C_0 (i) = \mathbb{E}^{h^{\dagger}}[{\mathrm{e}}^{-rT}(S_T - K)^+ \mid \Theta = i] = \mathrm{BSM}(S_0,K,T,r,\sigma_i).\end{equation}

Note that $\mathrm{BSM} (S_0, K, T, r, \sigma_i)$ is the Black–Scholes–Merton call price with volatility parameter $\sigma_i$ :

(4.6) \begin{equation} \mathrm{BSM}(S_0,K,T,r,\sigma_i) = S_0\Phi(d_1(i)) - K{\mathrm{e}}^{-rT}\Phi(d_2(i)),\end{equation}

where

\begin{equation*} d_1(i) = \frac{\ln({S_0}/{K}) + \big(r + \frac{1}{2}\sigma^2_i\big)T}{\sigma_i\sqrt{T}}, \qquad d_2(i) = d_1(i) - \sigma_i\sqrt{T}.\end{equation*}

Then, from (3.21) and (4.5), an analytical expression for the call price is

(4.7) \begin{equation} C_0 = \sum^{N}_{i = 1} \mathrm{BSM} (S_0, K, T, r, \sigma_i) p_i.\end{equation}

4.2. Finite mixture model of negative-shifted gamma distributions

From the empirical studies in Section 5, the skewness of the logarithmic returns of Bitcoin in USD is negative. Consequently, a shifted Gamma distribution may not be suitable to model the logarithmic returns of Bitcoin in USD since the skewness of a shifted Gamma distribution is positive [Reference Gerber and Shiu26, p. 117]. Consequently, instead of considering a finite mixture model of shifted Gamma distributions, we suppose that the logarithmic returns follow a finite mixture model of negative-shifted Gamma distributions. Specifically, we assume that under the real-world probability measure $\mathbb{P}$ , for each $i = 1, 2, \ldots, N$ and each $t = 1, 2, \ldots, T$ ,

(4.8) \begin{equation} R_t \mid \{ \Theta = i \} \overset{{\mathrm{d}}}{=} k_i - Y_i,\end{equation}

where $Y_i$ follows a Gamma distribution $\mathrm{Ga} (\alpha_i, \beta_i)$ with shape parameter $\alpha_i$ and rate parameter $\beta_i$ ; $\overset{{\mathrm{d}}}{=}$ means equality in distribution. In this case, for each $i = 1, 2, \ldots, N$ and each $t = 1, 2, \ldots, T$ , the conditional distribution of $R_t$ given $\{ \Theta = i \}$ under $\mathbb{P}$ is a negative-shifted Gamma distribution with shape parameter $\alpha_i$ , rate parameter $\beta_i$ , and shifted parameter $k_i$ (i.e. $\mathrm{NegSGa} (\alpha_i, \beta_i, k_i)$ ), where $k_i$ , $\alpha_i$ , and $\beta_i$ are all positive.

Let $G (x \mid \alpha, \beta)$ denote the CDF of a Gamma distribution with shape parameter $\alpha$ and rate parameter $\beta$ evaluated at $x > 0$ . Then, for each $i = 1, 2, \ldots, N$ , the ith-component CDF $F_{i} (x)$ is given by

(4.9) \begin{equation} F_i (x) = 1 - G (k_i - x \mid \alpha_i, \beta_i), \quad x < k_i.\end{equation}

Similarly to [Reference Gerber and Shiu26, (4.1.3)], for each $i = 1, 2, \ldots, N$ the MGF of $F_i (x)$ evaluated at $\eta$ is given by

(4.10) \begin{equation} M_i (\eta) = \bigg ( \frac{\beta_i}{\beta_i + \eta} \bigg )^{\alpha_i} {\mathrm{e}}^{k_i \eta}, \quad \eta > - \beta_i.\end{equation}

From (3.14) and (4.10),

(4.11) \begin{equation} {\mathrm{e}}^r = \bigg ( \frac{\beta_i + h^{\dagger} (i)}{\beta_i + h^{\dagger} (i) + 1} \bigg )^{\alpha_i}{\mathrm{e}}^{k_i}, \quad i = 1, 2, \ldots, N.\end{equation}

Take, for each $i = 1, 2, \ldots, N$ ,

(4.12) \begin{equation} \beta^{\dagger}_i = \beta_i + h^{\dagger} (i).\end{equation}

Then, from (4.11) and (4.12),

(4.13) \begin{equation} \beta^{\dagger}_i = \frac{1}{{\mathrm{e}}^{({k_i - r})/{\alpha_i}} - 1}, \quad k_i > r.\end{equation}

Indeed, (4.12) and (4.13) ensure the existence of an equivalent martingale measure $\mathbb{P}^{h^{\dagger}}$ provided that $k_i > r$ . The condition that $k_i > r$ is required to ensure that $\beta^{\dagger}_i > 0$ .

From (3.16), (4.10), and (4.12), for each $i = 1, 2, \ldots, N$ the conditional MGF of $R_t$ given $\Theta = i$ under $\mathbb{P}^{h^{\dagger}}$ evaluated at $\eta$ is given by

(4.14) \begin{equation} M^{h^{\dagger}}_{i} (\eta) = \bigg ( \frac{\beta^{\dagger}_i}{\beta^{\dagger}_i + \eta} \bigg )^{\alpha_i}{\mathrm{e}}^{k_i \eta}, \quad \eta > - \beta^{\dagger}_i.\end{equation}

Consequently, under $\mathbb{P}^{h^{\dagger}}$ , the conditional distribution of $R_t$ given $\Theta = i$ is a negative-shifted Gamma distribution with shape parameter $\alpha_i$ , rate parameter $\beta^{\dagger}_i$ , and shifted parameter $k_i$ . Since $R_1, R_2, \ldots, R_T \mid \Theta$ are (conditionally) independent under $\mathbb{P}^{h^{\dagger}}$ , $R_{1, T} \mid \{ \Theta = i \} \sim \mathrm{NegSGa} (\alpha_i T, \beta^{\dagger}_i, k_i T)$ under $\mathbb{P}^{h^{\dagger}}$ for each $i = 1, 2, \ldots, N$ .

Define

(4.15) \begin{equation} \kappa \,:\!=\, \ln \bigg ( \frac{K}{S_0} \bigg).\end{equation}

Then, similarly to [Reference Gerber and Shiu26, (4.1.7)],

(4.16) \begin{align} C_0(i) & = \mathbb{E}^{h^{\dagger}}[{\mathrm{e}}^{-rT}(S_T - K)^+ \mid \Theta = i] \nonumber \\ & = S_0G(k_iT - \kappa \mid \alpha_iT, \beta^{\dagger}_i + 1) - K{\mathrm{e}}^{-rT}G(k_iT - \kappa \mid \alpha_iT, \beta^{\dagger}_i) \nonumber \\ & = \mathrm{NSGa}(S_0,K,T,r,\alpha_i,\beta^{\dagger}_i,k_i), \quad \mbox{say}. \end{align}

Note that for (4.16) to hold, it is required that $k_i T > \kappa$ for each $i = 1, 2, \ldots, N$ .

From (3.21) and (4.16), an analytical expression for the call price is

(4.17) \begin{equation} C_0 = \sum^{N}_{i = 1} \mathrm{NSGa} (S_0, K, T, r, \alpha_i, \beta^{\dagger}_i, k_i) p_i.\end{equation}

From (4.13), $\beta^{\dagger}_i$ does not depend on $\beta_i$ for each $i = 1, 2, \ldots, N$ . Consequently, the call price $C_0$ in (4.17) does not depend on $\beta_i$ .

4.3. Finite mixture model of negative-shifted inverse Gaussian distributions

Let $J (x \mid a, b)$ denote the CDF of an inverse Gaussian (IG) distribution with parameters a and b evaluated at $x > 0$ . Then, from [Reference Gerber and Shiu26, (4.2.1)],

(4.18) \begin{equation} J (x \mid a, b) = \Phi\bigg({-}\frac{a}{\sqrt{2x}\,} + \sqrt{2bx}\bigg) + {\mathrm{e}}^{2a\sqrt{b}}\Phi\bigg({-}\frac{a}{\sqrt{2x}\,} - \sqrt{2bx}\bigg),\end{equation}

where $\Phi(\cdot)$ is the CDF of a standard normal distribution. That is, if $\mu^{\scriptstyle{sig}}$ and $\lambda$ are the mean and the shape parameter of the IG distribution, respectively, then $\mu^{\scriptstyle{sig}} = {a}/{2\sqrt{b}}$ and $\lambda = {a^2}/{2}$ .

We suppose that, under $\mathbb{P}$ , for each $i = 1, 2, \ldots, N$ and each $t = 1, 2, \ldots, T$ ,

(4.19) \begin{equation} R_t \mid \{ \Theta = i \} \overset{{\mathrm{d}}}{=} k_i - Y_i,\end{equation}

where $Y_i$ follows an IG distribution with parameters $a_i$ and $b_i$ . Note that $a_i$ , $b_i$ , and $k_i$ are all positive. In this case, the conditional distribution of $R_t$ given $\{ \Theta = i \}$ under $\mathbb{P}$ is a negative-shifted inverse Gaussian distribution with IG parameters $a_i$ and $b_i$ , as well as the shifted parameter $k_i$ . That is, under $\mathbb{P}$ , $R_t \mid \{ \Theta = i \} \sim \mathrm{NegSIG} (a_i, b_i, k_i)$ .

Then, for each $i = 1, 2, \ldots, N$ , the ith-component CDF $F_i (x)$ is given by

(4.20) \begin{equation} F_i (x) = 1 - J (k_i - x \mid a_i, b_i), \quad x < k_i.\end{equation}

Similarly to [Reference Gerber and Shiu26, (4.2.3)], for each $i = 1, 2, \ldots, N$ the MGF of $F_i (x)$ evaluated at $\eta$ is given by

(4.21) \begin{align} M_i(\eta) = \exp\left(a_i(\sqrt{b_i} - \sqrt{b_i + \eta}) + k_i\eta\right), \quad \eta > - b_i.\end{align}

From (3.14) and (4.21),

(4.22) \begin{equation} r = a_i \left(\sqrt{b_i + h^{\dagger} (i)} - \sqrt{b_i + h^{\dagger} (i) + 1}\right) + k_i.\end{equation}

Take, for each $i = 1, 2, \ldots, N$ ,

(4.23) \begin{equation} b^{\dagger}_i = b_i + h^{\dagger} (i).\end{equation}

Then, from (4.22) and (4.23),

(4.24) \begin{equation} \sqrt{b^{\dagger}_i} - \sqrt{b^{\dagger}_i + 1} = \frac{r - k_i}{a_i}.\end{equation}

Consequently, from (4.23), the existence of an equivalent martingale measure $\mathbb{P}^{h^{\dagger}}$ becomes the existence of a solution to the non-linear equation in (4.24) provided that $b_i > - h^{\dagger} (i)$ . The latter condition is required to ensure that $b^{\dagger}_i > 0$ .

From (3.16), (4.21), and (4.23), for each $i = 1, 2, \ldots, N$ the conditional MGF of $R_t$ given $\Theta = i$ under $\mathbb{P}^{h^{\dagger}}$ evaluated at $\eta$ is given by

(4.25) \begin{equation} M^{h^{\dagger}}_{i}(\eta) = \exp\Big(a_i\Big(\sqrt{b^{\dagger}_i} - \sqrt{b^{\dagger}_i + \eta}\Big) + k_i\eta\Big), \quad \eta > - b^{\dagger}_i.\end{equation}

Consequently, under $\mathbb{P}^{h^{\dagger}}$ , the conditional distribution of $R_t$ given $\Theta = i$ is a negative-shifted IG distribution with parameters $a_i$ and $b^{\dagger}_i$ , as well as the shifted parameter $k_i$ . Since $R_1, R_2, \ldots, R_T \mid \Theta$ are (conditionally) independent under $\mathbb{P}^{h^{\dagger}}$ , $R_{1, T} \mid \{ \Theta = i \} \sim \mathrm{NegSIG} (a_i T, b^{\dagger}_i, k_i T)$ under $\mathbb{P}^{h^{\dagger}}$ for each $i = 1, 2, \ldots, N$ . Then, similarly to [Reference Gerber and Shiu26, (4.2.7)],

(4.26) \begin{align} C_0(i) & = \mathbb{E}^{h^{\dagger}}[{\mathrm{e}}^{-rT}(S_T - K)^+ \mid \Theta = i] \nonumber \\ & = S_0J(k_iT - \kappa \mid a_iT,b^{\dagger}_i + 1) - K{\mathrm{e}}^{-rT}J(k_iT - \kappa \mid a_iT,b^{\dagger}_i), \nonumber \\ & = \mathrm{NSIG}(S_0, K, T, r, a_i, b^{\dagger}_i, k_i). \end{align}

Note that for (4.26) to hold, it is required that $k_i T > \kappa$ for each $i = 1, 2, \ldots, N$ .

From (3.21) and (4.26), an analytical expression for the call price is

(4.27) \begin{equation} C_0 = \sum^{N}_{i = 1} \mathrm{NSIG} (S_0, K, T, r, a_i, b^{\dagger}_i, k_i) p_i.\end{equation}

This analytical pricing formula is up to solving the non-linear equation for $b^{\dagger}_i$ in (4.24). Furthermore, from (4.24), $b^{\dagger}_i$ does not depend on $b_i$ for each $i = 1, 2, \ldots, N$ . Consequently, the call price $C_0$ in (4.27) does not depend on $b_i$ .

4.4. Hybrid finite mixture model

In this subsection, a hybrid finite mixture model of normal, negative-shifted Gamma and negative inverse Gaussian distributions is considered. The hybrid model allows for the flexibility that the component distributions may come from different parametric classes of distributions. It provides a feasible way to describe both model uncertainty and parameter uncertainty. The former is interpreted as uncertainty about the parametric form of distributions. Specifically, to describe model uncertainty we consider more than one parametric form for the component distributions and use the mixing mechanism in a finite mixture model to take a ‘weighted average’ of the component distributions across different parametric distribution classes. The idea may perhaps be related to the use of Bayesian averaging to capture model uncertainty in the context of smooth ambiguity pioneered in [Reference Klibanoff, Marinacci and Mukerji43]. To incorporate parameter uncertainty, the mixing mechanism in a finite mixture model is used to take a ‘weighted average’ of the component distributions within the same parametric class of distributions. Informally speaking, model uncertainty is incorporated through considering the variation or heterogeneity between different parametric classes of distributions. Parameter uncertainty is captured by considering the variation or heterogeneity within the same parametric class of distributions.

Suppose, for each $i = 1, 2, \ldots, N_1$ , the ith-component CDF $F_{i} (x)$ is the CDF of a normal distribution with mean $\mu_i$ and variance $\sigma^2_i$ . Also, for each $i = N_1 + 1, N_1 + 2, \ldots, N_1 + N_2$ , the ith-component CDF $F_{i} (x)$ is the CDF of a negative-shifted Gamma distribution with shape parameter $\alpha_i$ , rate parameter $\beta_i$ , and shifted parameter $k_i$ . Assume, for each $i = N_1 + N_2 + 1, N_1 + N_2 + 2, \ldots, N$ , the ith-component CDF $F_{i} (x)$ is the CDF of a negative-shifted IG distribution with IG parameters $a_i$ and $b_i$ , and the shifted parameter $k_i$ . Then the hybrid finite mixture model can be written as

(4.28) \begin{align} F(x) & = \sum^{N_1}_{i=1}p_iN(x \mid \mu_i,\sigma^2_i) + \sum^{N_1+N_2}_{i=N_1+1}p_i(1 - G(k_i - x \mid \alpha_i,\beta_i)) \nonumber \\ & \quad + \sum^{N}_{i = N_1 + N_2 + 1}p_i(1 - J(k_i - x \mid a_i, b_i))\end{align}

for $x < \min \{ k_{N_1+1}, k_{N_1 + 2}, \ldots, k_{N} \}$ .

Using the results in Sections 4.1–4.3, an analytical expression for the call price is

(4.29) \begin{align} C_0 & = \sum^{N_1}_{i=1}\mathrm{BSM}(S_0,K,T,r,\sigma_i)p_i + \sum^{N_1+N_2}_{i = N_1+1}\mathrm{NSGa}(S_0,K,T,r,\alpha_i,\beta^{\dagger}_i,k_i)p_i \nonumber \\ & \quad + \sum^{N}_{i=N_1+N_2+1}\mathrm{NSIG}(S_0,K,T,r,a_i,b^{\dagger}_i,k_i)p_i.\end{align}

Again, the analytical pricing formula in (4.29) is up to solving the non-linear equation for $b^{\dagger}_i$ in (4.24). The existence of an equivalent martingale measure $\mathbb{P}^{h^{\dagger}}$ follows directly from the respective discussions in Sections 4.1–4.3.

5.1. Estimation procedures

The ECF estimation method [Reference Feuerverger and McDunnough21, Reference Feuerverger and Mureika22, Reference Knight and Yu44, Reference Madan and Seneta48, Reference Madan and Seneta49, Reference Singleton61, Reference Tran79, Reference Yu83] has been applied to estimate unknown parameters of finite mixture models of normal distributions [Reference Bryant and Paulson10, Reference Tran79, Reference Yu83]. Its key idea is to estimate the unknown parameters by minimizing the ‘distance’ between the empirical characteristic function and its theoretical counterpart. It was noted in [Reference Yu83] that the ECF estimation method can be considered when the maximum likelihood estimation method is difficult and a closed-form solution of the (theoretical) characteristic function is available. This is why the ECF estimation method is adopted here. Specifically, for the cases of a finite mixture model of negative-shifted Gamma distributions and a finite mixture model of negative-shifted IG distributions, it is possible that the observed logarithmic returns may not be bounded above by the shifted parameters. This renders the maximum likelihood estimation method intractable. Furthermore, we note that closed-form expressions for the (theoretical) characteristic functions of a negative-shifted Gamma distribution and a negative-shifted IG distribution are available and that the (theoretical) characteristic function of a finite mixture model is a mixture of the (theoretical) characteristic functions of individual (theoretical) characteristic functions. Consequently, closed-form expressions for the (theoretical) characteristic functions of the finite mixture model of negative-shifted Gamma distributions and the finite mixture model of negative-shifted IG distributions are obtained. The ECF estimation method is intimately linked with the generalized method of moment (GMM) estimation [Reference Hansen35]. See [Reference Yu83] for a detailed discussion on the link between the GMM method and the ECF method. [Reference Gerber and Shiu26] proposed the use of the method of moments to estimate unknown parameters in a shifted Gamma distribution and a shifted IG distribution, and derived the respective closed-form estimators. In the following, the ECF estimation method is briefly discussed in the context of a general finite mixture model.

Let $\{ R_1, R_2, \ldots, R_n \}$ denote a set of n observations about the logarithmic returns of Bitcoin. Then the ECF of $\{ R_1, R_2, \ldots, R_n \}$ evaluated at $u \in {\mathbb{R}}$ is defined as

(5.1) \begin{equation} C_n (u) \,:\!=\, \frac{1}{n} \sum^{n}_{j = 1} \exp (\mathrm{i} u R_j),\end{equation}

where $\mathrm{i} \,:\!=\, \sqrt{-1}$ .

Let $C (u;\, \varphi)$ denote the (theoretical) characteristic function of the general finite mixture distribution F(x) in (2.1) evaluated at $u \in {\mathbb{R}}$ under $\mathbb{P}$ , where $\varphi$ is a vector of unknown parameters in the finite mixture distribution $F (x) \,:\!=\, F (x;\, \varphi)$ and $\varphi \in \Phi \subset {\mathbb{R}}^k$ for some Euclidean space ${\mathbb{R}}^k$ . For each $j = 1, 2, \ldots, N$ , let $C_j (u; \varphi)$ denote the (theoretical) characteristic function of the jth-component CDF $F_j (x)$ of the general finite mixture distribution F(x) evaluated at $u \in {\mathbb{R}}$ under $\mathbb{P}$ . That is, $C_j (u;\, \varphi)$ is the conditional (theoretical) characteristic function of $R_t$ given $\{ \Theta = j \}$ evaluated at $u \in {\mathbb{R}}$ under $\mathbb{P}$ . Note that for each $j = 1, 2, \ldots, N$ , the vector of unknown parameters in $F_j (x)$ is obtained by putting some components of $\varphi$ equal to zero. To simplify the notation, we still adopt $\varphi$ to denote the vector of unknown parameters in $F_j (x)$ . By putting $\eta = \mathrm{i} u$ in (2.7), we have

(5.2) \begin{align} C (u, \varphi) & = \sum^{N}_{j = 1} p_{j} C_{j} (u, \varphi), \end{align}

(5.3) \begin{align} C (u, \varphi) & = \mathbb{E} [{\mathrm{e}}^{\mathrm{i} u R_t} ] = \int_{\mathcal{S}}{\mathrm{e}}^{\mathrm{i} u x}\,{\mathrm{d}} F(x), \end{align}

(5.4) \begin{align} C_j(u,\varphi) & = \mathbb{E}[{\mathrm{e}}^{\mathrm{i}uR_t} \mid \Theta = j] = \int_{\mathcal{S}_j}{\mathrm{e}}^{\mathrm{i}ux}\,{\mathrm{d}} F_j (x). \end{align}

Note that $\mathcal{S}$ is the support of the probability density function (PDF) of the distribution F(x); $\mathcal{S}_j$ is the support of the PDF of the distribution $F_j (x)$ for each $j = 1, 2, \ldots, N$ . In the case of a finite mixture of negative-shifted Gamma distributions, a closed-form expression for $C_j (u, \varphi)$ is obtained by putting $\eta = \mathrm{i} u$ in (4.10). In the case of a finite mixture of negative-shifted IG distributions, a closed-form expression for $C_j (u, \varphi)$ is obtained by putting $\eta = \mathrm{i} u$ in (4.21).

Then the ECF estimation method is to estimate $\varphi$ by solving the following minimization problem:

(5.5) \begin{eqnarray} \min_{\varphi \in \Phi} \int^{\infty}_{-\infty} | C_n (u, \varphi) - C (u, \varphi) |^2 g (u)\, {\mathrm{d}} u,\end{eqnarray}

where g(u) is some continuous weighting function. Here, as in [Reference Yu83], we choose $g (u) = \exp (\!- u ^ 2)$ .

In practice, it is often the case that there is no closed-form expression for the integral in (5.5). Consequently, the integral is evaluated via some numerical integration procedures. Here we shall use the function integrate in R to compute the integral numerically, which adopts adaptive quadrature of functions and involves a set of discrete points for the variable u of the integration. To solve the minimization problem in (5.5), the differential evolution algorithm is implemented using the DEoptim R package [Reference Ardia, Mullen, Peterson, Ulrich and Boudt2]. Differential evolution is an evolutionary global optimization algorithm that adopts biology-inspired operations on a population to minimize an objective function over successive generations [Reference Ardia, Mullen, Peterson, Ulrich and Boudt2, Reference Mitchell53, Reference Price, Storn and Lampinen56, Reference Storn and Price72]. More detail about the use of differential evolution to estimate the finite mixture models considered in the current paper is provided in Section 5.2. It will be seen that human or subjective judgement may need to be exercised even when this advanced optimization algorithm is adopted. Finally, it may be noted that the ECF estimator enjoys some desirable asymptotic properties such as strong consistency and asymptotic normality. See [Reference Yu83] for a more detailed discussion.

5.2. Estimation results

A dataset on daily Bitcoin-USD (BTC-USD) adjusted close prices is used in empirical studies on the eight models in Table D.1. Each of Models I–V has two mixing components. The time period for the BTC-USD prices data is from 1 January 2020 to 29 May 2025 (1976 observations). The dataset was extracted from Yahoo Finance using R. It covers some periods of the Covid pandemic. In [Reference Siu64, Reference Siu66], Bitcoin price data covering some periods of the Covid pandemic were considered. However, [Reference Siu64] focused on risk evaluation based on Value at Risk and Expected Shortfall. [Reference Siu66] studied the impacts of a long memory in the conditional volatility and conditional non-normality on pricing Bitcoin options. With the Bitcoin logarithmic returns data (1975 observations), where the empirical studies were performed using R. The estimation results are discussed in this subsection.

The estimation results for the five two-component finite mixture models (Models I–V) are provided in Tables D.2–D.6 in the supplementary material, respectively. Specifically, for each of the five two-component finite mixture models, the estimates of the unknown parameters, the optimal value(s) of the objective function corresponding to the parameter estimates, and the lower and upper boundaries of the parameters are presented. For the finite mixture model of normal distributions (Model I), one set of parameter estimates is provided based on one set of lower and upper boundaries of the parameters. For the other four finite mixture models with non-normal mixing components (Models II–V), two sets of parameter estimates are provided based on two sets of lower and upper boundaries of the parameters. In each of Tables D.2–D.6, ‘fn’ represents the optimal value of the objective function in differential evolution evaluated at the parameter estimates. All the parameter estimates in Tables D.2–D.6 are in daily units. It will be seen from the results that the parameter estimates obtained from differential evolution may be different when different sets of lower and upper boundaries of the parameters are adopted.

In Table D.2, the lower and upper boundaries of $\mu_1$ and $\mu_2$ are given by $-100$ and 100 times the sample mean of the logarithmic returns, respectively. The upper boundaries of $\sigma_1$ and $\sigma_2$ are 100 times the sample standard deviation of the logarithmic returns. From the estimation results, around $34.87\%$ weight is allocated to the first component. The first component corresponds to a ‘bad’ economic regime with a negative mean return and high volatility. The second component corresponds to a ‘good’ economic regime with a positive mean return and low volatility.

There are two sets (Sets 1 and 2) of estimation results in each of Tables D.3–D.6. In each of Tables D.3 and D.4, the upper boundaries of parameters in Set 1 are 5 times their respective method of moments estimates from the sample logrithmic returns data. The upper boundaries of parameters in Set 2 are 10 times their respective method of moments estimates from the sample data. The estimates of $a_i$ and $b_i$ , $i = 1, 2$ , in Table D.4 are obtained from the respective estimates of $\mu^{\scriptstyle{sig}}_i$ and $\lambda_i$ . In each of Tables D.5 and D.6, for each of Sets 1 and 2, the lower and upper boundaries of $\mu$ are $-100$ and 100 times the sample mean of the logarithmic returns, respectively; the upper boundaries of $\sigma$ in Sets 1 and 2 are 100 times the sample standard deviation of the logarithmic returns. Again, the upper boundaries of negative SGa and SIG parameters in Set 1 are five times their respective method of moment estimates from the sample logrithmic returns data. The upper boundaries of negative SGa and SIG parameters in Set 2 are 10 times their respective method of moments estimates from the sample data. The estimates of a and b in Table D.6 are obtained from the respective estimates of $\mu^{\scriptstyle{sig}}$ and $\lambda$ . There are two reasons why the upper boundaries of negative SGa and SIG parameters are set as either 5 or 10 times their respective method of moment estimates. Firstly, some of the parameters may become large if their upper boundaries are large. Specifically, the parameters $\beta_i$ (or $\beta$ ) in the negative SGa distributions may become large if their upper boundaries are large. Secondly, the optimal value ‘fn’ may not improve significantly even if the upper boundaries are large. However, the computational time increases if the upper boundaries are large. In each of Tables D.3–D.6, the estimation results from Sets 1 and 2 could be quite different. This illustrates that the parameter estimates obtained from differential evolution depend on the chosen boundaries of the parameters. We may need to exercise subjective judgements even though the advanced optimization algorithm is adopted.

In Table D.3, the optimal value ‘fn’ (i.e. the objective function in (5.5)) increases from $1.699\,928\times10^{-10}$ to $2.519\,128\times10^{-9}$ when moving from Set 1 to Set 2. This indicates that the optimal value ‘fn’ may not improve when parameter bounds are relaxed. This happens because the solver was not able to find the global optimum and led to suboptimal convergence. Although the estimates of $\beta_1$ in Sets 1 and 2 are large, they are still interior optimums (i.e. strictly less than the respective upper boundaries). Although the parameter estimates in Set 1 provide a better fit to the data than those in Set 2 according to the optimal value ‘fn’, the parameter estimates in Set 2 will be used to calculate the option prices in Section 6 to illustrate the mixture effect. In Table D.4, the optimal value ‘fn’ increases from $1.512\,852 \times 10^{-9}$ to $3.604\,965\times10^{-9}$ when moving from Set 1 and Set 2. This also indicates that the optimal value ‘fn’ may not improve when the upper boundaries become larger. Since the parameter estimate of b in Set 2 looks more reasonable than that in Set 1, the former will be used to calculate the option prices in Section 6. In Table D.5, the parameter estimates in Set 1 provide a better fit to the data than those in Set 2. Furthermore, the parameter estimate of $\sigma$ in Set 2 is very small (0.002 944 167), which may result in unreasonably small option prices. Consequently, the parameter estimates in Set 1 will be used to calculate the option prices in Section 6. In Table D.6, the estimates of p in Sets 1 and 2 are $0.012\,423\,47$ and $0.578\,201\,806$ , respectively. To illustrate the mixture effect, the parameter estimates in Set 2 may be better than those in Set 1. However, the estimate of $\sigma$ in Set 2 is very small ( $0.002\,358\,598$ ), which may result in unreasonably small option prices. Furthermore, the parameter estimates in Set 1 provide a better fit to the data than those in Set 2. Consequently, the parameter estimates in Set 1 will be used to calculate the option prices in Section 6.

Tables D.7–D.9 present the estimation results for the three non-mixture models (Models VI–VIII), respectively.

In Table D.7, both the maximum likelihood estimates (MLE) and the ECF estimates are provided. For the ECF estimates, the lower and upper boundaries of $\mu$ are $-100$ and 100 times the sample mean of the logarithmic returns, respectively. The upper boundary of $\sigma$ is 100 times the sample standard deviation of the logarithmic returns. The optimal value ‘fn’ (i.e. the objective function in (5.5)) evaluated at the ECF estimates is smaller than the optimal value ‘fn’ evaluated at the MLE, though the optimal values evaluated at both the MLE and the ECF estimates are reasonably small. This indicates that the ECF estimation method may be better than the MLE method when estimating the mean and standard deviation of a normal distribution. However, it is possible that the solver was not able to find the global optimum and led to suboptimal convergence. The ECF estimates will be used to calculate the option prices in Section 6. In each of Tables D.8 and D.9, the upper boundaries for the parameters in Sets 1 and 2 are 5 and 10 times their respective method of moment estimates, respectively. The estimates of a and b in Table D.9 are obtained from those of $\mu^{\scriptstyle{sig}}$ and $\lambda$ . In Table D.8, the parameter estimates in Set 2 provide a better fit to the data and will be used to calculate the option prices in Section 6. In Table D.9, the parameter estimates in Set 1 provide a better fit to the data and will be used to calculate the option prices in Section 6. Using ‘fn’ as a criterion, the fitting performances of Models I–VIII look similar, though Model VIII is the best among them.