Interest rate derivatives for the fractional Cox-Ingersoll-Ross model

Abstract

We obtain the bond price formula for the fractional Cox-Ingersoll-Ross model. Then we obtain option price formula for the bond. Finally we apply it to derive option price formula in fractional Heston model.

Keywords

Wick-Itô stochastic differential equation affine models fractional Cox-Ingersoll-Ross model interest rate bond price fractional Heston model stochastic volatility Monte Carlo method

1 Introduction and preliminaries

There are many connections between models in biology and finance. Feller (1951) reached at the square-root process as the weak limit of Galton-Watson branching process with immigration while studying a problem in genetics. Using the Feller’s square-root process, Cox et al. (1985) studied the theory of term structure of interest rates and the model is now known as the Cox-Ingersoll-Ross model. Overbeck and Ryden (1997) studied asymptotics of conditional least squares esimtimators of Cox-Ingersoll-Ross process from discrete observations using an auto-regressive type representation of the model with non-Gaussian error. Dehtiar et al. (2021) studied strong consistency for the maximum likelihood method and an alternative method of estimation of the drift parameters of the Cox-Ingersoll-Ross process based on continuous observations.

On the other hand, anomalous diffusions with long range dependence of the extended Vasicek type have been used in short-rate modeling (Shea, 1991; Backus and Zin, 1993; Booth and Tse, 1995; Tsay, 2000; Cajueiro and Tabak, 2007, 2008; Hao et al., 2014; Zhang et al., 2017). Backus and Zin (1993) studied long memory inflation uncertainty as an evidence from the term structure of interest rates. Biagini et al. (2013) studied a fractional credit model with long range dependent default rate. Ohashi (2009) studied fractional term structure models and obtained no arbitrage and consistency properties. Recently, fractional Cox-Ingersoll-Ross model has received some attention in finance. Fractional Cox-Ingersoll-Ross model catches the realistic feature of the positivity of the interest rates in addition to preserving the long memory. Mishura and Yurchenko-Tytarenko (2018) studied hitting probability of fractional Cox-Ingersoll-Ross model which involves long memory. Fallah et al. (2019) studied fractional CIR model and obtained pricing formula for double barrier option under transaction cost with Hurst index H ∈ (2/3, 1). They also verified the effect of the parameters of the model on the value of the option. Gao et al. (2022) studied Euler-Maruyama approximation method for fractional CIR model and obtained H-order rate of the root mean square error by Lamperti transformation using Malliavin calculus. Mpanda et al. (2022) studied a generalization of fractional CIR model extending Mishura and Yurchenko-Tytarenko (2011). Araneda (2020) studied a fractional and mixed-fractional Constant Elasticity of Volatility (CEV) model. Recently, Bishwal (2022a) introduced many non-affine fractional models including fractional Black-Karasinski model of term structure of interest rates.

Let $(Ω, F, {F_{t}}_{t \geq 0}, P)$ be a stochastic basis. Let ${W_{t}^{H}}_{t \geq 0}$ be a persistent (Hurst parameter H > 0.5) fractional Brownian motion with the filtration ${F_{t}}_{t \geq 0}$ and α > 0, β > 0, σ > 0 are the unknown parameters. Recall that a fractional Brownian motion (fBM) has the covariance ${\tilde{C}}_{H} (s, t) = \frac{1}{2} [s^{2 H} + t^{2 H} - | s - t |^{2 H}], s, t > 0 .$

For H > 0.5 the process has long range dependence or long memory. For H ≠ 0.5, the process is neither a Markov process nor a semimartingale. For H = 0.5, the process reduces to standard Brownian motion.

Recently Ichiba et al. (2021, 2022) studied generalized fractional Brownian motion (GFBM). A generalized fractional Brownian motion is a Gaussian self-similar process whose increments are not necessarily stationary. It appears in the scaling limit of a shot-noise process with a power law shape function and non-stationary noises with a power law variance function. They studied semimartingale properties of the mixed process made up of an independent Brownian motion and a GFBM for the persistent Hurst parameter.

On the stochastic basis the fractional Ornstein-Uhlenbeck (fOU) process X_t is defined satisfying the Wick-Itô stochastic differential equation ${dX}_{t} = θ X_{t} dt + {dW}_{t}^{H}, t \geq 0, X_{0} = ξ$ (1.1) where ${W_{t}^{H}}$ is a fractional Brownian motion with H > 1/2 with the filtration ${F_{t}}_{t \geq 0}$ and θ < 0 is the unknown parameter.

Davydov (1970) gave an AR(1) approximation of the fOU model: $\begin{matrix} y_{j} = ρ y_{j - 1} + v_{j}, (1 - L)^{H - 1 / 2} v_{j} = ϵ_{j}, \\ ρ = 1, y_{0} = 0, j = 1, 2 . \dots, T \end{matrix}$ where L is the lag operator, ϵ_j∼ i.i.d. (0, σ²) with $E (ϵ_{j}^{4}) < \infty$ , whereas {v_j} is a stationary long-memory process generated by $v_{j} = (1 - L)^{H - 1 / 2} ϵ_{j} = \sum_{k = 0}^{\infty} \frac{Γ (k + H - 1 / 2)}{Γ (H - 1 / 2) Γ (k + 1)} ϵ_{j - k} .$

Davydov (1970) proved that $\frac{λ_{H}^{1 / 2}}{σ T^{H}} y_{[Tt]} \to^{D} W_{t}^{H} as T \to \infty .$ Next we present the martingale approximation to the nonsemimartingale fBm.

Define $\begin{matrix} κ_{H} : = 2 H Γ (3 / 2 - H) Γ (H + 1 / 2), k_{H} (t, s) : = κ_{H}^{- 1} (s (t - s))^{\frac{1}{2} - H} \\ λ_{H} : = \frac{2 H Γ (3 - 2 H) Γ (H + \frac{1}{2})}{Γ (3 / 2 - H)}, v_{t} \equiv v_{t}^{H} : = λ_{H}^{- 1} t^{2 - 2 H}, M_{t}^{H} : = \int_{0}^{t} k_{H} (t, s) {dW}_{s}^{H} . \end{matrix}$ From Norros et al. (1999) it is well known that $M_{t}^{H}$ is a Gaussian martingale, called the fundamental martingale whose variance function 〈M^H〉_t is $v_{t}^{H}$ . The natural filtration of the martingale M^H coincides with the natural filtration of the fBm W^H since $W_{t}^{H} : = \int_{0}^{t} K_{H} (t, s) {dM}_{s}^{H}$ holds for H ∈ (1/2, 1) where $K_{H} (t, s) : = H (2 H - 1) \int_{s}^{t} r^{h - \frac{1}{2}} (r - s)^{H - \frac{3}{2}}, 0 \leq s \leq t$ and for H = 1/2, the convention K_1/2 ≡ 1 is used. Observe that the increments of $M_{t}^{H}$ are independent of $W_{t}^{H}$ and $Cov (M_{t}^{H}, M_{s}^{H}) = λ_{H}^{- 1 / 2} (s \land t)^{2 - 2 H}, Cov (W_{t}^{H}, M_{t}^{H}) = λ_{H}^{1 / 2},$ $M_{t}^{H} = b_{H} \int_{0}^{t} u^{1 / 2 - H} {dW}_{u} =^{D} M_{t}^{H}, b_{H} = \sqrt{2 - 2 H} λ_{H} .$ Define $Q_{t} : = \frac{d}{{dv}_{t}} \int_{0}^{t} k_{H} (t, s) X_{s} ds .$ It is easy to see that $Q_{t} = \frac{λ_{H}}{2 (2 - 2 H)} {t^{2 H - 1} Z_{t} + \int_{0}^{t} r^{2 H - 1} {dZ}_{s}} .$ (1.2) The process X admits the representation $X_{t} = \int_{0}^{t} K_{H} (t, s) {dZ}_{s} .$ The natural filtration generated by the fundamental semimartingale process $Z_{t} = θ \int_{0}^{t} Q_{s} {dv}_{s} + M_{t}^{H}$ and the process X coincide (Kleptsyna and Le Breton, 2002). The available information for X and Z are strictly equivalent.

Let the realization {X_t, 0 ≤ t ≤ T} or equivalently {Z_t, 0 ≤ t ≤ T} be denoted by $Z_{0}^{T}$ . Let $P_{θ}^{T}$ be the measure generated on the space (C_T, B_T) of continuous functions on [0, T] with the associated Borel σ-algebra B_T generated under the supremum norm by the process $Z_{0}^{T}$ and $P_{0}^{T}$ be the standard Wiener measure. Applying fractional Girsanov formula, when θ is the true value of the parameter, $P_{θ}^{T}$ is absolutely continuous with respect to $P_{0}^{T}$ and the Radon-Nikodym derivative (likelihood) of $P_{θ}^{T}$ with respect to $P_{0}^{T}$ based on $Z_{0}^{T}$ is given by $L_{T} (θ) : = \frac{{dP}_{θ}^{T}}{{dP}_{0}^{T}} (Z_{0}^{T}) = exp {θ \int_{0}^{T} Q_{t} {dZ}_{t} - \frac{θ^{2}}{2} \int_{0}^{T} Q_{t}^{2} {dv}_{t}} .$ (1.3) Consider the score function, the derivative of the log-likelihood function, which is given by $Y_{T} (θ) : = \int_{0}^{T} Q_{t} {dZ}_{t} - θ \int_{0}^{T} Q_{t}^{2} {dv}_{t} .$ (1.4) Consider the problem of testing hypotheses $H_{0} : θ = θ_{0}$ against the alternative $H_{1} : θ = θ_{1}$ Under the hypothesis $H_{0}$ , the log-likelihood ratio process admits the representation $\begin{matrix} Λ_{T} : = log \frac{{dP}_{θ_{1}}^{T}}{{dP}_{θ_{0}}^{T}} (Z_{0}^{T}) \\ = {(θ_{1} - θ_{0}) \int_{0}^{T} Q_{t} {dM}_{t} - \frac{(θ_{1} - θ_{0})^{2}}{2} \int_{0}^{T} Q_{t}^{2} {dv}_{t}} . (1.5) \end{matrix}$ The Hellinger integral of order ɛ ∈ (- ∞ , ∞) is defined as $h_{T} (ɛ) = E_{θ_{0}} [exp (ɛ Λ_{T})] .$ Note that $h_{T} (ɛ) : = h_{T} (ɛ; P_{θ_{1}}^{T}, P_{θ_{0}}^{T}) = h_{T} (1 - ɛ; P_{θ_{0}}^{T}, P_{θ_{1}}^{T}) .$ With ɛ = 1, $h_{T} (1) = E_{θ_{0}} [exp (Λ_{T})] = E_{θ_{0}} [\frac{{dP}_{θ_{1}}^{T}}{{dP}_{θ_{0}}^{T}} (Z_{0}^{T})] .$ We have $\int_{0}^{T} Q_{t} {dM}_{t} = \frac{1}{2} (U_{T} + 2 θ V_{T} - T)$ (1.6) where $U_{T} : = \frac{λ_{H}}{2 - 2 H} Z_{T} \int_{0}^{T} t^{2 H - 1} {dZ}_{s}$ (1.7) and $V_{T} : = \int_{0}^{T} Q_{s}^{2} {dv}_{s} .$ (1.8)

To test the hypotheses $H_{0}$ and $H_{1}$ , the test statistic is given by $δ_{t}^{c, ɛ} = I (X_{t} > c) + ϵ I (X_{t} = c)$ where c ɛ (0, ∞) and ɛ ∈ [0, 1] are parameters of the test. From the Neyman-Pearson lemma, for any α ∈ [0, 1], there exists a test $δ_{t}^{c (α), ɛ (α)}$ of level α, where (c (α) , ɛ (α)) is some solution to the equation $δ_{t}^{c (α), ɛ (α)} = α$ with respect to (c, α). For any α ∈ [0, 1], the test $δ_{t}^{c (α), ɛ (α)}$ with ɛ (0) =1 is the most powerful test. A test $δ_{t}^{c (α), ɛ (α)}$ with ɛ (0) =1 is called the Neyman-Pearson test.

Let $φ_{T}^{H} (u) : = E exp (- {uV}_{T}), u > 0$ . Kleptsyna and Le Breton (2002) showed that $φ_{T}^{H} (u)$ is given by $φ_{T}^{H} (u) = E exp (- {uV}_{T}) = {\frac{4 sin π H \sqrt{θ^{2} + 2 u} e^{- θ T}}{π T D_{T}^{H} (θ; \sqrt{θ^{2} + 2 u})}}^{1 / 2}$ (1.9) where $\begin{matrix} D_{T}^{H} (θ; α) & : = & [α cosh (\frac{α}{2} T) - θ sinh (\frac{α}{2} T)]^{2} J_{- H} (\frac{α}{2} T) J_{H - 1} (\frac{α}{2} T) \\ - [α sinh (\frac{α}{2} T) - θ cosh (\frac{α}{2} T)]^{2} J_{1 - H} (\frac{α}{2} T) J_{H} (\frac{α}{2} T) (1.10) \end{matrix}$ for α > 0 and J_ν is the modified Bessel function of first kind of order ν.

Remark When H = 1/2, $J_{- 1 / 2} (z) = \sqrt{2 / (π z)} cos z$ , $J_{1 / 2} (z) = \sqrt{2 / (π z)} sin z$ .

In this case we have the well-known Cameron-Martin formula: $E exp (- \int_{0}^{T} W_{t}^{2} dt) = (cosh \sqrt{2} T)^{- 1 / 2} .$ (1.11) Tanaka (2014) studied distributions of quadratic functionals of the fractional Brownian motion based on a martingale approximation.

Bishwal (2011a) generalized the above Cameron-Martin type formula and obtained the joint moment generating function of U_T and V_T.

Theorem 1.1 Let $Ψ_{T}^{H} (z_{1}, z_{2}) : = E exp (z_{1} V_{T} + z_{2} U_{T}), z_{1}, z_{2} \in ℂ$ . Then $Ψ_{T}^{H} (z_{1}, z_{2})$ exists for |z_i| ≤ δ, 1 = 1,2 for some δ > 0 and is given by

$Ψ_{T}^{H} (z_{1}, z_{2}) = exp (\frac{- θ T}{2}) {[\frac{4 (sin π H) γ}{π T D_{T}^{H} (γ; z_{2})}]}^{1 / 2}$ (1.12)

where γ : = (θ² + 2z₁) ^1/2 and we choose the principal branch of the square root and

$\begin{matrix} D_{T}^{H} (γ; z_{2}) & : = & [γ cosh (\frac{γ}{2} T) - (θ - 2 z_{2}) sinh (\frac{γ}{2} T)]^{2} J_{- H} (\frac{γ}{2} T) J_{H - 1} (\frac{γ}{2} T) \\ - [γ sinh (\frac{γ}{2} T) - (θ - 2 z_{2}) cosh (\frac{γ}{2} T)]^{2} J_{1 - H} (\frac{γ}{2} T) J_{H} (\frac{γ}{2} T) (1.13) \end{matrix}$

and J_ν is the modified Bessel function of first kind of order ν. Note that for equally spaced data

$v_{t_{i}} - v_{t_{i - 1}} = λ_{H}^{- 1} {(\frac{T}{n})}^{2 - 2 H} [i^{2 - 2 H} - (i - 1)^{2 - 2 H}] .$ (1.14)

We have

$\begin{matrix} Q_{t} & = & \frac{d}{{dv}_{t}} \int_{0}^{t} k_{H} (t, s) X_{s} ds = κ_{H}^{- 1} \frac{d}{{dv}_{t}} \int_{0}^{t} s^{1 / 2 - H} (t - s)^{1 / 2 - H} X_{s} ds \\ = & κ_{H}^{- 1} λ_{H} t^{2 H - 1} \frac{d}{dt} \int_{0}^{t} s^{1 / 2 - H} (t - s)^{1 / 2 - H} X_{s} ds \\ = & κ_{H}^{- 1} λ_{H} t^{2 H - 1} \int_{0}^{t} \frac{d}{dt} s^{1 / 2 - H} (t - s)^{1 / 2 - H} X_{s} ds \\ = & κ_{H}^{- 1} λ_{H} t^{2 H - 1} \int_{0}^{t} s^{1 / 2 - H} (t - s)^{- 1 / 2 - H} X_{s} ds . \end{matrix}$ (1.15)

The process Q depends continuously on X and therefore, the discrete observations of X does not allow one to obtain the discrete observations of Q. The process Q can be approximated by

${\tilde{Q}}_{n} = κ_{H}^{- 1} λ_{H} n^{2 H - 1} \sum_{j = 0}^{n - 1} j^{1 / 2 - H} (n - j)^{- 1 / 2 - H} X_{j} .$ (1.16)

It is easy to show that ${\tilde{Q}}_{n} \to Q_{t}$ almost surely as n→ ∞ (Tudor and Viens, 2007).

Define a new partition 0 ≤ u₁ < u₂ < u₃ < ⋯ < um_i = t_i, i = 1, 2, ⋯ , n . Define

${\tilde{Q}}_{t_{i}} = κ_{H}^{- 1} λ_{H} t_{i}^{2 H - 1} \sum_{j = 1}^{m_{i}} u_{j}^{1 / 2 - H} (u_{m_{i}} - u_{j})^{- 1 / 2 - H} X_{u_{j}} (u_{j} - u_{j - 1}),$ (1.17)

i = 1, 2, ⋯ , n . It is easy to show that ${\tilde{Q}}_{t_{i}} \to Q_{t}$ almost surely as m_i→ ∞ for each i = 1, 2, ⋯ , n . In practice, one can use this approximate observation in the calculation of the bond price and option price.

2 Bond price

The Cox-Ingersoll-Ross (CIR) model, also called the square-root process, is used as a short rate mean reverting model in term structure of interest rates and a stochastic volatility process in the Heston model (Bishwal, 2022b).

Let $(Ω, F, {F_{t}}_{t \geq 0}, P)$ be a stochastic basis on which is defined the Cox-Ingersoll-Ross process {X_t} satisfying the Itô stochastic differential equation ${dX}_{t} = (α + β X_{t}) dt + 2 {\sqrt{X}}_{t} {dW}_{t}, t \geq 0, X_{0} = 1$ (2.1) where {W_t} is a standard Wiener process with the filtration ${F_{t}}_{t \geq 0}$ , α > 0 and β < 0 are the unknown parameters.

The true transition density is given by $q (t, x, y, β) : = - 2 β {(\frac{y}{x})}^{α - \frac{1}{2}} \frac{e^{(\frac{1}{2} - α) β t}}}{1 - e^{β t}} \exp [\frac{2 β (x + y)}{e^{- β t} - 1}] I_{ν} (\frac{- 2 β \sqrt{x y}}{\sinh (- \frac{1}{2} β t)})$ (2.2) where I_ν is the modified Bessel function with index ν which is noncentral chi-square density. The invariant density as t→ ∞ is gamma. Sankaran (1963) obtained approximation to the non-central chi-square distribution.

Let the continuous realization {X_t, 0 ≤ t ≤ T} be denoted by $X_{0}^{T}$ . Let $P_{β, α}^{T}$ be the measure generated on the space (C_T, B_T) of continuous functions on [0, T] with the associated Borel σ-algebra B_T generated under the supremum norm by the process $X_{0}^{T}$ and let $P_{0}^{T}$ be the standard Wiener measure. It is well known that when β, α are the true values of the parameters, $P_{β, α}^{T}$ is absolutely continuous with respect to $P_{0}^{T}$ and the Radon-Nikodym derivative (likelihood) of $P_{β, α}^{T}$ with respect to $P_{0}^{T}$ based on $X_{0}^{T}$ is given by $L_{T} (β, α) : = \frac{{dP}_{β, α}^{T}}{{dP}_{0}^{T}} (X_{0}^{T}) = exp {\int_{0}^{T} \frac{α + β X_{t}}{4 X_{t}} {dX}_{t} - \int_{0}^{T} \frac{(α + β X_{t})^{2}}{8 X_{t}} dt} .$ (2.3) Consider the score function, the derivative of the log-likelihood function, which is given by $γ_{T} (β, α) : = {\int_{0}^{T} \frac{α + β X_{t}}{4 X_{t}} {dX}_{t} - \int_{0}^{T} \frac{(α + β X_{t})^{2}}{8 X_{t}} dt} .$ (2.4)

Let $I_{T} : = \int_{0}^{T} X_{t} dt, J_{T} : = \int_{0}^{T} X_{t}^{- 1} dt .$ (2.5)

The following characteristic function of I_T is closely associated with Levy’s stochastic area formula and is well known from Brownian motion literature and also from the work of Cox et al. (1985).

Let $φ_{T} (u) : = E exp (iu I_{T}), u \in ℝ$ . Then $E exp (iu I_{T}) = exp (\frac{2 iu}{β + γ coth \frac{γ T}{2}}) {[cosh \frac{γ T}{2} + \frac{β}{γ} sinh \frac{γ T}{2}]}^{- 1}$ (2.6) where γ : = (β² - 2iu) ^1/2. Observe that $E exp (iu I_{T}) = exp (\frac{2 iu (e^{\frac{γ T}{2}} - e^{\frac{- γ T}{2}})}{e^{\frac{γ T}{2}} (γ + β) + e^{- \frac{γ T}{2}} (γ - β)}) 2 {[e^{\frac{γ T}{2}} (γ + β) + e^{- \frac{γ T}{2}} (γ - β)]}^{- 1} .$ (2.7) Further $lim_{T \to \infty} E (e^{- {uI}_{T}}) = lim_{T \to \infty} exp (- \frac{ab}{2 σ} T (\sqrt{1 + \frac{4 u σ}{b^{2}}} - 1)) .$ (2.8) $\begin{matrix} E exp (- u X_{T} - v I_{T}) \\ = & {(\frac{2 ρ e^{(b - ρ) T / 2}}{2 σ u (1 - e^{- ρ T}) + (ρ - b) e^{- ρ T} + (ρ + b)})}^{a / σ} \\ exp (\frac{u ((ρ + b) e^{- ρ T} + (ρ - b) + 2 v (1 - e^{- ρ T})}{2 σ u (1 - e^{- ρ T}) + (ρ - b) e^{- ρ T} + (ρ + b)}) \end{matrix}$ (2.9) where $ρ = : \sqrt{b^{2} + 4 σ v}$ . Finally $\begin{matrix} E exp (- u J_{T}) = \frac{Γ (k + \frac{ν}{2} + \frac{1}{2})}{Γ (ν + 1)} (\frac{x}{α})^{- k} β^{\frac{ν}{2} + \frac{1}{2}} \\ exp (\frac{b}{2 σ} [at - \frac{2 x}{e^{bt} - 1}])_{1} F_{1} (k + \frac{ν}{2} + \frac{12}{,} ν + 1, β) \end{matrix}$ (2.10) where $k = : \frac{a}{2 σ}, α = : \frac{{be}^{bt}}{σ (e^{bt} - 1)}, ν = : \frac{1}{σ} \sqrt{(α - σ)^{2} + 4 u σ}$ (2.11) and ₁F₁ is Kummer’s confluent hypergeometric function.

The Vasicek model for short rate is given by ${dr}_{t} = a (b - r_{t}) dt + σ {dW}_{t}, t \geq 0$ (2.12) The price of a zero coupon bond at time t maturing at time T is given by $P (t, T) = A (t, T) e^{- B (t, T) r (t)}$ (2.13) where $B (t, T) : = \frac{1 - e^{- a (T - t)}}{a}, A (t, T) : = exp (\frac{(B (t, T) - T + t) (a^{2} b - σ^{2} / 2)}{a^{2}} - \frac{σ^{2} B (t, T)^{2}}{4 a}) .$ (2.14)

For the CIR model ${dX}_{t} = a (b - X_{t}) dt + σ \sqrt{X_{t}} {dW}_{t}$ (2.15) the price at time t of a zero-coupon bond that pays $1 at time T is given by $P (t, T) = E_{\tilde{P}} (exp (- \int_{t}^{T} X_{t} dt)) = A (t, T) e^{- B (t, T) X_{t}}$ where $B (t, T) : = \frac{2 (e^{γ (T - t) - 1})}{(γ + a) (e^{γ (T - t) - 1}) + 2 γ}, A (t, T) : = {(\frac{2 γ e^{(a + γ) (T - t) / 2}}{(γ + a) (e^{γ (T - t) - 1}) + 2 γ})}^{2 ab / σ^{2}}$ (2.16) where $γ = \sqrt{a^{2} + 2 σ^{2}}$ and $\tilde{P}$ is the risk-neutral measure.

Note that B (t, T) can be written as $B (t, T) = \frac{2 (e^{γ (T - t) - 1})}{(γ + a) (e^{γ (T - t) - 1}) + 2 γ} .$ (2.17)

Now we focus on the fractional CIR (fCIR) model with H = 0.5

Theorem 2.1 For the fCIR model with H = 0.5 ${dX}_{t} = a (b - X_{t}) dt + σ \sqrt{X_{t}} {dW}_{t}^{H}$ (2.18) the price at time t of a zero-coupon bond that pays $1 at time T is given by $P^{H} (t, T) = E_{\tilde{P}} (exp (- \int_{t}^{T} X_{t} dt)) = A^{H} (t, T) e^{- B^{H} (t, T) X_{t}}$ (2.19) where $B^{H} (t, T) = \frac{2 (e^{γ (T - t) - 1})}{(γ + a) (e^{γ (T - t) - 1}) + 2 γ}, A^{H} (t, T) = {(\frac{2 γ e^{(a + γ) (T - t) / 2}}{(γ + a) (e^{γ (T - t) - 1}) + 2 γ})}^{2 ab / σ^{2}}$ (2.20) where $γ = \sqrt{a^{2} + 2 σ^{2}}$ and $\tilde{P}$ is the risk-neutral measure.

Now consider the fCIR model ${dX}_{t} = a (b - X_{t}) dt + σ \sqrt{X_{t}} {dW}_{t}^{H},$ (2.21) Then by Proposition 5.7 of Buchmann and Kluppelberg (2006), we have X_t = f (Y_t) where ${dY}_{t} = a (b - Y_{t}) dt + {dW}_{t}^{H}, Y_{0} = f^{- 1} (X_{0}), t \in [0, T]$ (2.22) and f (x) = sign (x) σ²x²/4.

Fink et al. (2013) studied the bond price formula expression of which involves fractional Riemann-Liouville type integral.

Theorem 2.2 For the fCIR model with H ≥ 0.5 ${dX}_{t} = a (b - X_{t}) dt + σ \sqrt{X_{t}} {dW}_{t}^{H},$ (2.23) let a = -2θ, b - 1/(2θ) , σ = 2, $R_{T} : = \int_{0}^{T} Q_{s} {dv}_{s}$ and $φ_{T}^{H} (u) : = E exp (- {uR}_{T}), u > 0$ . Then $φ_{T}^{H} (u)$ is given by $φ_{T}^{H} (u) = E exp (- {uR}_{T}) = {\frac{4 sin π H \sqrt{θ^{2} + 2 u} e^{- θ T}}{π T D_{T}^{H} (θ; \sqrt{θ^{2} + 2 u})}}^{1 / 2}$ (2.24) where $\begin{matrix} D_{T}^{H} (θ; α) & : = & [α cosh (\frac{α}{2} T) - θ sinh (\frac{α}{2} T)]^{2} J_{- H} (\frac{α}{2} T) J_{H - 1} (\frac{α}{2} T) \\ - [α sinh (\frac{α}{2} T) - θ cosh (\frac{α}{2} T)]^{2} J_{1 - H} (\frac{α}{2} T) J_{H} (\frac{α}{2} T) (2.25) \end{matrix}$ for α > 0 and J_ν is the modified Bessel function of first kind of order ν.

With H ≥ 0.5, for the general case, the bond price is given by (2.19) where the denominators in B^H and A^H are replaced by $π {TD}_{T}^{H} (a; γ) / (4 sin π H)$ .

3 Interest rate derivative

First we consider the modeling of stock price. Let $(Ω, F, {F_{t}}_{t \geq 0}, P)$ be a stochastic basis on which is defined the following process {S_t, t ≥ 0} where {W_t} is a standard Wiener process with the filtration ${F_{t}}_{t \geq 0}$ .

For the Black-Scholes model for stock price ${dS}_{t} = μ S_{t} dt + σ S_{t} {dW}_{t},$ (3.1) using Itô formula to log S_t, one obtains the solution $S_{t} = S_{0} exp {((μ - {\frac{12}{σ}}^{2}) t + σ W_{t}}$ (3.2) which is known as geometric Brownian motion. The parameter μ is known as the mean rate of return and σ as the volatility.

Call option (buyer’s option) at time t is the expected discounted (at the risk free interest rate r) pay-off $C_{t} = E_{\tilde{P}} [e^{- r (T - t)} max (S_{T} - K, 0) | F_{t}]$ (3.3) where K is the strike price of the option and T is the time of maturity of the option and $\tilde{P}$ is the risk-neutral measure.

Put option (seller’s option) at time t is the expected discounted (at the risk free interest rate r) pay-off $P_{t} = E_{\tilde{P}} [e^{- r (T - t)} max (K - S_{T}, 0) | F_{t}]$ (3.4) where K is the strike price of the option and T is the time of maturity of the option.

Using no arbitrage principle, Black-Scholes derived the PDE which is given by $\frac{\partial C_{t}}{\partial t} + {rS}_{t} \frac{\partial C_{t}}{\partial S_{t}} + \frac{1}{2} σ^{2} S_{t}^{2} \frac{\partial^{2} C_{t}}{\partial^{2} S_{t}} - {rS}_{t} = 0 .$ (3.5) Black and Scholes (1973) calculated the above expectation by solving the PDE for C_t and is known as the famous Black-Scholes option pricing formula. The Black-Scholes option price formulae for European call and put options are given respectively by: $C_{t} = E_{\tilde{P}} [e^{- r (T - t)} (S_{T} - K)^{+} | F_{t}] = S_{t} Φ (d_{1}) - K e^{- r (T - t)} Φ (d_{2})$ (3.6) $P_{t} = E_{\tilde{P}} [e^{- r (T - t)} (K - S_{T})^{+} | F_{t}] = K e^{- r (T - t)} Φ (- d_{2}) - S_{t} Φ (- d_{1})$ (3.7) where $d_{1} : = \frac{log (\frac{S}{K}) + (r + \frac{σ^{2}}{2}) (T - t)}{σ \sqrt{T - t}}, d_{2} : = \frac{log (\frac{S}{K}) + (r - \frac{σ^{2}}{2}) (T - t)}{σ \sqrt{T - t}} = d_{1} - σ \sqrt{T - t}$ (3.8) and Φ (·) is the cummulative distribution function of standard normal distribution. In a risk neutral world, where all expectations are calculated under the risk-neutral measure or the martingale measure, the stock price S_t at time t follows the following linear Itô stochastic differential equation, known as the Black-Scholes model ${dS}_{t} = r S_{t} dt + σ S_{t} {dW}_{t}, t \geq 0$ (3.9) where {W_t} is a standard Brownian motion, r is the risk-free interest rate and σ is the volatility. A simple application of Itô’s formula to log S_t provides the exact solution of the equation given by $S_{t} = S_{0} exp {(r - \frac{1}{2} σ^{2}) t + σ W_{t}}, t \geq 0$ (3.10) where S₀ is the initial price of the stock. S_t is called Geometric Brownian motion. One can generate the path of the stock price by the exact method and the Euler method and calculate the price of the European call option at time 0 based on both the methods.

In the Black-Scholes model, the interest rate r and the volatility σ are constant. In practice, both interest rates and the volatility are stochastic processes. First we consider stochastic interest rate, which is known as short rate Vasicek model: dr_t = a (b - r_t) dt + σdW_t, a, b, σ > 0 .

For the Vasicek model of short rate, the interest rate derivative, European call option is given by $C = LP (0, s) Φ (h) - KP (0, T) Φ (h - σ_{p})$ (3.11) where L is the bond principal, s is the bond maturity, T is the option maturity, K is the strike price, P (0, s) is the bond price of zero-coupon bond with maturity s, $h : = \frac{1}{σ_{P}} ln \frac{LP (0, s)}{P (0, T) K} + \frac{σ_{P}}{2}, σ_{P} : = \frac{σ}{a} (1 - e^{- a (s - T)}) \sqrt{\frac{1 - e^{- 2 aT}}{2 a}} .$ (3.12) When a = 0, $σ_{P} = σ (s - T) \sqrt{T}$ . Parameter estimation in Vasicek model is extensively studied in Bishwal (2008).

The drawback of the Vasicek model is that interest rate can be negative since the transition density of the process is normal. Next we consider a positive interest rate model which also serves as stochastic volatility model. Let $(Ω, F, {F_{t}}_{t \geq 0}, P)$ be a stochastic basis on which is defined the Cox-Ingersoll-Ross process ${σ_{t}^{2}}$ satisfying the Itô stochastic differential equation $d σ_{t}^{2} = (1 + 2 θ σ_{t}^{2}) dt + 2 σ_{t} {dW}_{t}, t \geq 0,$ (3.13) where {W_t} is a standard Wiener process with the filtration ${F_{t}}_{t \geq 0}$ and consider the classical direct estimation problem where θ < 0 is the unknown parameter to be estimated on the basis of discrete observations of the process ${σ_{t}^{2}}$ at times 0 = t₀ < t₁ < ⋯ t_n = T with $t_{i} - t_{i - 1} = \frac{T}{n} = Δ, i = 1, 2 \dots, n$ . There are two approaches to the asymptotic framework: 1) Δ → 0, n → ∞ , 2) Δ fixed, n→ ∞.

For the moment assume that a continuous realization ${σ_{t}^{2}, 0 \leq t \leq T}$ be denoted by $D_{0}^{T}$ . Let $P_{θ}^{T}$ be the measure generated on the space $(C_{T}, B_{T})$ of continuous functions on [0, T] with the associated Borel σ-algebra $B_{T}$ generated under the supremum norm by the process $D_{0}^{T}$ and let $P_{0}^{T}$ be the standard Wiener measure. It is well known that when θ is the true value of the parameter $P_{θ}^{T}$ is absolutely continuous with respect to $P_{0}^{T}$ and the Radon-Nikodym derivative (likelihood) of $P_{θ}^{T}$ with respect to $P_{0}^{T}$ based on $D_{0}^{T}$ is given by $L_{T} (θ) : = \frac{{dP}_{θ}^{T}}{{dP}_{0}^{T}} (D_{0}^{T}) = exp {θ \int_{0}^{T} d σ_{t}^{2} - \frac{θ^{2}}{2} \int_{0}^{T} σ_{t}^{2} dt} .$ (3.14) Consider the score function, the derivative of the log-likelihood function, which is given by $γ_{T} (θ) : = \int_{0}^{T} d σ_{t}^{2} - θ \int_{0}^{T} σ_{t}^{2} dt .$ (3.15)

A solution of the estimating equation γ_T (θ) =0 provides the maximum likelihood estimate (MLE) ${\hat{θ}}_{T} : = \frac{σ_{T}^{2} - σ_{0}^{2} - T / 2}{\int_{0}^{T} σ_{t}^{2} dt} .$ (3.16) The minimum contrast estimate (MCE) is given by ${\tilde{θ}}_{T} : = - \frac{T}{2 \int_{0}^{T} σ_{t}^{2} dt} .$ Minimum contrast estimation based on direct observation was studied in Bishwal (2006, 2011a). Berry-Esseen inequalities for fractional Black-Karasinski model was studied in Bishwal (2022a). Tanaka (2013) showed the advantage of this estimator over least squares estimator in terms of efficiency. Note that the volatility which is given by the CIR process is not observed.

Consider the Heston stochastic volatility model ${dS}_{t} = μ S_{t} dt + \sqrt{X_{t}} S_{t} {dW}_{t},$ (3.17) $d X_{t} = (α + β X_{t}) dt + σ \sqrt{X_{t}} {dZ}_{t}$ (3.18) where {W_t} is a standard Brownian motion independent of another standard Brownian motion {Z_t}, and the parameters α, σ are positive and β is negative. The integrated volatility is given by $I_{T} : = \int_{0}^{T} X_{t} dt .$ (3.19)

Denote $s_{t} : = ln S_{t}, Δ s_{t_{i - 1}}^{2} : = (s_{t_{i}} - s_{t_{i - 1}})^{2} .$ (3.20) The realized volatility is defined as $R_{n, T} : = \sum_{i = 1}^{n} Δ s_{t_{i - 1}}^{2} .$ (3.21) It is well known that $\underset{n \to \infty}{P - \lim} R_{n, T} = I_{T} .$ (3.22)

Thus the realized volatility estimates the integrated volatility. Bishwal (2011b) obtained several higher order new estimators of integrated volatility using kernel method. Maximum quasi-likelihood estimation in fractional Levy stochastic volatility model was studied in Bishwal (2011c) based on observations of the asset return data.

Note that the volatility which is given by the CIR process is not observed. In the following we obtain nonparametric estimators of the minimum contrast estimator of the mean reversion parameter in the Heston model using approximations to ${\tilde{θ}}_{T}$ .

Hence from the definition of MCE, the approximate minimum contrast estimate (AMCE) of θ would be ${\tilde{θ}}_{n, T} : = - \frac{T}{2 R_{n, T}} .$ (3.23)

The characteristic function of I_T is closely associated with Levy’s stochastic area formula and is well known from Brownian motion literature and also from the work of Cox et al. (1985).

The European call option for the CIR model is given by $C_{t} = P (t, s) χ^{2} (2 \frac{log (A (t, s) / K)}{B (T, s)} [φ + ψ + B (T, s)]; \frac{4 ab}{σ^{2}}, \frac{2 φ^{2} X_{t} e^{γ (T - t)}}{φ + ψ + B (T, s)})$ $- KP (t, T) χ^{2} (2 \frac{log (A (t, s) / K)}{B (T, s)} [φ + ψ]; \frac{4 ab}{σ^{2}}, \frac{2 φ^{2} X_{t} e^{γ (T - t)}}{φ + ψ})$ (3.24) where $φ = : \frac{2 γ}{σ^{2} (e^{γ (T - t) - 1})}, ψ = : \frac{a + γ}{σ^{2}}$ (3.25) and χ² (x ; d, λ) is the noncentral chi-square distribution with d degrees of freedom and noncentrality parameter λ.

Next consider the special Heston stochastic volatility model (Heston, 1993): ${dS}_{t} = μ S_{t} dt + \sqrt{X_{t}} S_{t} {dW}_{t}$ (3.26) $d X_{t} = (1 + 2 θ X_{t}) dt + 2 \sqrt{X_{t}} {dZ}_{t}$ (3.27) where {W_t} is a standard Brownian motion independent of another standard Brownian motion {Z_t} and θ < 0. The advantage of this representation with α = 1, β = 2θ and σ = 1 is that it represents the squared Vasicek process. The integrated volatility is given by (3.19).

Integrated volatility has to be estimated on the basis of discrete observations of the price process {S_t} at times 0 = t₀ < t₁ < ⋯ t_n = T with $t_{i} - t_{i - 1} = \frac{T}{n}, i = 1, 2 \dots, n$ . Let s_t = log S_t be the log-price process. Then ${ds}_{t} = (μ - \frac{1}{2} X_{t}) dt + \sqrt{X_{t}} {dW}_{t} .$ (3.28) Thus the drift term depends on the volatility.

Consider the Heston model with correlated noises: Under the real world measure P, we have ${dS}_{t} = μ S_{t} dt + \sqrt{S_{t}} {dW}_{t}^{1},$ (3.29) ${dV}_{t} = κ (θ - V_{t}) dt + σ \sqrt{V_{t}} {dW}_{t}^{2},$ (3.30) ${dW}_{t}^{1} {dW}_{t}^{2} = ρ dt .$ (3.31)

The correlation ρ < 0 is called the leverage effect.

Under the risk neutral measure $\tilde{P}$ , we have ${dS}_{t} = r S_{t} dt + \sqrt{V_{t}} S_{t} d {\tilde{W}}_{t}^{1}$ (3.32) ${dV}_{t} = κ^{*} (θ^{*} - V_{t}) dt + σ \sqrt{V_{t}} d {\tilde{W}}_{t}^{2}$ (3.33) $d {\tilde{W}}_{t}^{1} d {\tilde{W}}_{t}^{2} = ρ dt, κ^{*} = κ + λ, θ^{*} = \frac{κ θ}{κ + λ}$ (3.34) $d {\tilde{W}}_{t}^{1} = {dW}_{t}^{1} + θ_{t} dt, d {\tilde{W}}_{t}^{2} = {dW}_{t}^{2} + k \sqrt{V_{t}} dt$ (3.35) and the Radon-Nikodym derivative is given by $\frac{d \tilde{P}}{dP} = exp {- \frac{1}{2} \int_{0}^{t} (θ_{s}^{2} + k^{2} V_{s}) ds - \int_{0}^{t} θ_{s} {dW}_{s}^{1} - \int_{0}^{t} k \sqrt{V_{s}} {dW}_{s}^{2}},$ (3.36) where $θ_{t} : = (μ - r) / \sqrt{V_{t}}$ is the risk premium.

The characteristic function of s_T under $\tilde{P}$ is given by $ψ_{T} (u) = E_{\tilde{P}} (exp (i s_{T} u)) = A (u, T) e^{B (u, T)}$ (3.37) where $\begin{matrix} A (u, T) & : = & exp (iu (s_{0} + rT)), \\ B (u, T) & : = & \frac{- (u^{2} + iu) (1 - e^{γ T}) V_{0}}{2 γ - (γ - (κ - ρ σ ui) (1 - e^{γ T})} \\ - & \frac{κ θ}{σ^{2}} [2 log (\frac{2 γ - (γ - κ - ρ σ ui) (1 - e^{γ T})}{2 γ}) + (γ - κ - ρ σ ui) T] (3.38) \end{matrix}$ with $γ = : \sqrt{(κ - ρ σ ui)^{2} + (u^{2} + iu)}$ and we choose the principal branch of the square root.

Option price is defined as $E_{\tilde{P}, t}^{T} [e^{r (T - t)} ϖ (T)]$ (3.39) where ϖ (T) is the payoff at time T where the expectation is under the risk neutral measure $\tilde{P}$ and r is the risk-free interest rate. The call option with strike price K with expiry T is given by $C (S_{t}, V_{t}, t, T) = S_{t} P_{1} - K e^{- r (T - t)} P_{2}$ (3.40) where for j = 1, 2 $P_{j} (x, V_{t}, T, K) : = \frac{12}{+} \frac{1}{π} \int_{0}^{\infty} Re (\frac{e^{- iu ln K} f_{j} (x, V_{t}, T, u)}{iu}) du,$ (3.41) $x : = ln S_{t}, f_{j} (x, V_{t}, T, u) : = exp {G (T - t, u) + D (T - t, u) V_{t} + iux},$ (3.42) $G (T - t, u) : = rui (T - t) + \frac{a}{σ^{2}} [(b_{j} - ρ σ ui + d) (T - t) - 2 ln (\frac{1 - {ge}^{d (T - t)}}{1 - g})],$ (3.43) $D (T - t, u) : = \frac{b_{j} - ρ σ ui + d}{σ^{2}} (\frac{1 - e^{d (T - t)}}{1 - {ge}^{d (T - t)}}), g = \frac{b_{j} - ρ σ ui + d}{b_{j} - ρ σ ui - d},$ (3.44) $d = \sqrt{(ρ σ ui - b_{j})^{2} - σ^{2} (2 ν_{j} ui - u^{2})}, ν_{1} : = \frac{12}{,} ν_{2} : = - \frac{12}{,} a : = κ θ, b_{1} : = κ + λ - ρ σ, b_{2} : = κ + λ .$ (3.45) Carr and Madan (1999) used fast Fourier transform (FFT) to evaluate the integral in this option price formula. On can also use numerical quadrature, like adaptive Simpson’s rule to evaluate the integral. Then Matlab or other programs can be used to calculate the option price. Alternatively, one can calculate the option price by using the Monte Carlo method after using second order discretization of the Heston model (Glasserman, 2004, pp. 356-357.)

Theorem 3.1 For the fCIR model ${dX}_{t} = a (b - X_{t}) dt + σ \sqrt{X_{t}} {dW}_{t}^{H}$ (3.46) the European call option is given by $C_{t}^{H} = P^{H} (t, s) χ^{2} (2 \frac{log (A^{H} (t, s) / K)}{B^{H} (T, s)} [φ + ψ + B^{H} (T, s)]; \frac{4 ab}{σ^{2}}, \frac{2 φ^{2} X_{t} e^{γ (T - t)}}{φ + ψ + B^{H} (T, s)})$ (3.47) $- K P^{H} (t, T) χ^{2} (2 \frac{log (A^{H} (t, s) / K)}{B^{H} (T, s)} [φ + ψ]; \frac{4 ab}{σ^{2}}, \frac{2 φ^{2} X_{t} e^{γ (T - t)}}{φ + ψ})$ (3.48) where $φ = : \frac{2 γ}{σ^{2} (e^{γ (T - t) - 1})}, ψ = : \frac{a + γ}{σ^{2}}$ (3.49) and χ² (x ; d, λ) is the noncentral chi-square distribution with d degrees of freedom and noncentrality parameter λ, and $P^{H} (t, T) = φ_{t}^{H} (1)$ given in (2.24).

Next consider the fractional Heston model with correlated noises: Under the real world measure P, we have the asset price and volatility processes respectively given by ${dS}_{t} = μ S_{t} dt + \sqrt{V_{t}} S_{t} {dW}_{t}^{1},$ (3.50) ${dV}_{t} = a (b - V_{t}) dt + σ \sqrt{V_{t}} {dW}_{t}^{H, 2},$ (3.51) ${dW}_{t}^{1} {dW}_{t}^{H, 2} = ρ_{H} dt .$ (3.52) Note that we have long-memory (H > 1/2) only in the volatility process.

Theorem 3.2 The call option price for the fractional Heston model with strike price K is given by $C^{H} (S_{t}, V_{t}, t, T) = S_{t} P_{1}^{H} - K e^{- r (T - t)} P_{2}^{H}$ (3.53) where for j = 1, 2, $P_{j}^{H} (x, V_{t}, T, K) = \frac{12}{+} \frac{1}{π} \int_{0}^{\infty} Re (\frac{e^{- iu ln K} f_{j} (x, V_{t}, T, u)}{iu}) du,$ (3.54) $x = ln S_{t}, f_{j} (x, V_{t}, t, u) = exp {G^{H} (T - t, u) + D^{H} (T - t, u) V_{t} + iux},$ (3.55) $G^{H} (T - t, u) : = rui (T - t) + \frac{a}{σ^{2}} [(b_{j} - ρ σ ui + d) (T - t) - 2 ln (\frac{1 - {ge}^{d (T - t)}}{1 - g})],$ (3.56) $D^{H} (T - t, u) : = \frac{b_{j} - ρ σ ui + d}{σ^{2}} (\frac{1 - e^{d (T - t)}}{1 - {ge}^{d (T - t)}}), g : = \frac{b_{j} - ρ σ ui + d}{b_{j} - ρ σ ui - d},$ (3.57) $d : = \sqrt{(ρ σ ui - b_{j})^{2} - σ^{2} (2 ν_{j} ui - u^{2})}, ν_{1} = \frac{12}{,} ν_{2} = - \frac{12}{,} a : = κ θ, b_{1} : = κ + λ - ρ σ, b_{2} : = κ + λ .$ (3.58)

The characteristic function here is calculated by fractional Itô’s lemma and fractional Fokker-Plank forward equation.

4 Concluding remarks

1) Hawkes processes are an efficient generalization of the Poisson processes to model a sequence of arrivals over time of some types of events, that present self-exciting feature, in the sense that each arrival increases the rate of future arrivals for some period of time (Hawkes, 1971). This class of counting processes allows one to capture self-exciting phenomena in a more accurate way compared to inhomogeneous Poisson processes or Cox processes. In finance, they are accurate to model for example credit risk contagion, order book or microstructure noises’s feature of financial markets.

A Hawkes process is a counting process A_t with stochastic intensity λ_t given by $λ_{t} = μ + \int_{0}^{t} Φ (t - s) {dA}_{s}$ where μ > 0 and $Φ : ℝ \to ℝ_{+}$ are two parameters. The parameter μ > 0 is called the background intensity and the function Φ is called the excitation function. When Φ = 0, this a homogeneous Poisson process.

A fractional Hawkes process {A_H (t) , t > 0} with Hurst parameter H ∈ (1/2, 1) is defined as $A_{H} (t) = \frac{1}{Γ (H - \frac{1}{2})} \int_{0}^{t} u^{\frac{1}{2} - H} (\int_{u}^{t} τ^{H - \frac{1}{2}} (τ - u)^{H - \frac{3}{2}} d τ) dR (u)$ where $R (u) = \frac{A (u)}{{\sqrt{λ}}_{t}} - {\sqrt{λ}}_{t} u$ and A (u) is a Hawkes process with stochastic intensity λ_t.

A CIR process with Hawkes jumps has been proposed in Zhu (2014): $d X_{t} = (α + β X_{t}) dt + σ \sqrt{X_{t}} {dW}_{t} + γ {dN}_{t}$ (3.59) where (N_t) is a Hawkes process (Bacry et al., 2013a, 2013b). This model helps in understanding default clustering effect, i.e, one default tends to trigger more defaults.

2) A real asset price model should be of the following hybrid type with 14 parameters. We consider the hybrid stochastic volatility, stochastic interest rate, stochastic leverage and stochastic elasticity model under the risk neutral measure which is given by ${dS}_{t} = X_{t} dt + \sqrt{V_{t -}} S_{t} {dW}_{t} + ρ_{λ t} {dL}_{τ_{λ t}},$ (3.60) ${dV}_{t} = - λ V_{t} dt + υ_{λ t -} {dL}_{τ_{λ t}},$ (3.61) ${dX}_{t} = α (β - X_{t}) dt + σ X_{t}^{γ_{t}} {dW}_{t}^{H},$ (3.62) $d ρ_{t} = ((2 ζ - η) - η ρ_{t}) dt + θ \sqrt{(1 + ρ_{t}) (1 - ρ_{t})} {dZ}_{t},$ (3.63) $d ξ_{t} = κ (μ - ξ_{t}) dt + ς \sqrt{ξ_{t}} {dB}_{t},$ (3.64) $d γ_{t} = ϖ (ψ - δ)) dt + \sqrt{χ} {dM}_{t},$ (3.65) $d τ_{t} = ξ_{t -} dt .$ (3.66) where L_t is a Levy process, W^H is a subfractional Brownian motion, B_t, Z_t and M_t are standard Brownian motions. Here S_t is the asset price which a geometric jump-diffusion, V_t is the stochastic volatility which is a Levy O-U process, X_t is the stochastic interest rate which is a sub-fractional Chan-Karloyi-Longstaff-Sanders (CKLS) process, ρ_t is the stochastic leverage Jacobi (Beta) process, ξ_t is a volatility modulation (stochastic time change) of the driving Levy subordinator which is a Cox-Ingersoll-Ross (CIR) process, γ_t is the stochastic elasticity models which is another CIR process, and all the 14 parameters λ, α, β, σ, ξ, η, θ, κ, μ, ς, ϖ, ψ, δ, χ are positive (Bishwal, 2021).

References

Araneda,

A.A.

2020 A fractional and mixed-fractional CEV model, Journal of Computational and Applied Mathematics 363, 106–120.

Bacry, E., Delattre, S., Hoffmann, M., Muzy, J.F. 2013a. Modelling microstructure noise with mutually exciting point processes, Quantitative Finance 13, 65–77.

Bacry, E., Delattre, S., Hoffmann, M., Muzy, J.F. 2013b. Scaling limits of Hawkes processes and application to financial statistics, Stochastic Processes and their Applications 123, 2475–2499.

Backus,

D.K

, Zin,

S.E.

1993 Long memory inflation uncertainty: evidence from the term structure of interest rates, Journal of Money, Credit and Banking 25, 681–700.

Biagini,

, Fink,

, Kluppelberg,

2013 A fractional credit model with long range dependent default rate, Stochastic Processes and thir Applications 123, 1319–1347.

Bishwal,

J.P.N.

2006 Rate of weak convergence of the approximate minimum contrast estimators for discretely observed Ornstein-Uhlenbeck process, Statistics and Probability Letters 76(13), 1397–1409.

Bishwal,

J.P.N.

2008 Parameter Estimation in Stochastic Differential Equations, Berlin: Springer-Verlag.

Bishwal, J.P.N.

2011a. Minimum contrast estimation in fractional Ornstein-Uhlenbeck process: continuous and discrete sampling, Fractional Calculus and Applied Analysis 14(3), 375–410.

Bishwal, J.P.N.

2011b. Some new estimators of integrated volatility, Open Journal of Statistics 1(2), 74–80.

10.

Bishwal, J.P.N.

2011c. Maximum quasi-likelihood estimation in fractional Levy stochastic volatility model, Journal of Mathematical Finance 1(3), 12–15.

11.

Bishwal,

J.P.N.

2021 A new algorithm for approximate maximum likelihood estimation in sub-fractional Chan-Karolyi-Longstaff-Sanders model, Asian Journal of Probability and Statistics 13(3), 62–88.

12.

Bishwal, J.P.N.

2022a. Berry-Esseen inequalities for the fractional Black-Karasinski model of term structure of interest rates, Monte Carlo Methods and Applications 28(2), 111–124.

13.

Bishwal, J.P.N.

2022b. Parameter Estimation in Stochastic Volatility Models. Cham: Springer Nature, 2022.

14.

Black,

, Scholes,

1973 The pricing of options and corporate liabilities, Journal of Political Economy 81, 637–659.

15.

Booth,

G.G

, Tse,

Y.M.

1995 Long memory in interest rate futures market: A fractional cointegration analysis, Journal of Futures Markets 15(5), 573–584.

16.

Buchmann,

, Kluppelberg,

2006 Fractional integral equations and state space transforms, Bernoulli 12, 431–456.

17.

Cajueiro,

D.O

, Tabak,

B.M.

2007 Time-varying long-range dependence in US interest rates, Chaos, Solitions and Frls 34, 360–367 acta.

18.

Cajueiro,

D.O

, Tabak,

B.M.

2008 Long range dependence in interest rates and monetary policy, Physical Letters A 372, 181–184.

19.

Carr,

, Madan,

D.B.

1999 Option valuation using fast Fourier transform, Journal of Computational Finance 2(4), 61–73.

20.

Cox,

J.C

, Ingersoll,

J.E

, Ross,

S.A.

1985 A theory of term structure of interest rates, Econometrica 53, 385–407.

21.

Davydov,

1970 The invariance principle for stationary processes, Theory of Probability and its Applications 15, 487–498.

22.

Dehtiar,

, Mishura,

, Ralchenko,

2021 Two methods of estimation of the drift parameters of the Cox-Ingersoll-Ross process; Continuous observations, Communications in Statistics-Theory and Methods 50, 1–16.

23.

Fallah,

, Najafi,

A.R

, Mehrdoust,

2019 A fractional version of the Cox-Ingersoll-Ross interest rate model and pricing double barrier option with Hurst index H ∈ (2/3, 1). Communications in Statistics-Theory and Methods 48(9), 2254–2266.

24.

Feller,

1951 Two singular diffusion problems, Annals of Mathematics 54, 173–182.

25.

Fink,

, Kluppelberg,

, Zahle,

2013 Conditional distributions of processes related to fractional Brownian motion, Journal of Applied Probability 50, 166–183.

26.

Gao,

, Wang,

, Yang,

2022 The truncated Euler-Maruyama method for CIR model driven by fractional Brownian motion, Statistics and Probability Letters 189, 109573.

27.

Glasserman,

2004. Monte Carlo Methods in Financial Engineering. New York: Springer.

28.

Hao,

, Liu,

, Wang,

2014 Pricing credit default swap under fractional Vasicek interest rate model, Journal of Mathematical Finance 4, 10–20.

29.

Hawkes,

A.G.

1971 Spectra of some self-exciting and mutually exciting point processes, Biometrika 58, 83–90.

30.

Kleptsyna,

M.L

, Le Breton,

2002 Statistical analysis of the fractional Ornstein-Uhlenbeck type process, Statistical Inference for Stochastic Processes 5, 229–248.

31.

Mishura,

, Yurchenko-Tytarenko,

2018 Fractional Cox-Ingersoll-Ross model with non-zero mean, Modern Stochastics: Theory and Methods 5(1), 99–111.

32.

Mpanda,

M.M

, Mukeru,

, Mulaudzi,

2022 Genaralisation of of fractional Cox-Ingersoll-Ross process, Results in Applied Mathematics 15, 100322.

33.

Norros,

, Valkeila,

, Virtamo,

1999 An elementary approach to a Girsanov formula and other analytical results on fractional Brownian motion, Bernoulli 5, 571–587.

34.

Ohashi,

2009 Fractional term structure models: no arbitrage and consistency, Annals of Applied Probability 19(4), 1533–1580.

35.

Overbeck,

, Ryden,

1997 Estimation in the Cox-Ingersoll-Ross model, Econometric Theory 13, 430–461.

36.

Schroder,

1989 Computing the constant elasticity of variance option price formula, Journal of Finance 44(1), 211–219.

37.

Sankaran,

1963 Approximation to the non-central chi-square distribution, Biometrika 50, 199–204.

38.

Ichiba,

, Pang,

, Taqqu,

M.S.

2021 Semimartingale properties of generalized fractional Brownian motion and its mixtures with applications in finance, ArXiv:2012.00975v2.v.

39.

Ichiba,

, Pang,

, Taqqu,

M.S.

2022 Path properties of a genaralized fractional Brownian motion, Journal of Theoretical Probability 35, 550–574.

40.

Tanaka,

2013 Distributions of the maximum likelihood and minimum contrast estimators associated with the fractional Ornstein-Uhlenbeck process, Statistical Inference for Stochastic Processes 36, 173–192.

41.

Tanaka,

2014 Distributions of quadratic functionals of the fractional Brownian motion based on a martingale approximation, Econometric Theory 30, 1078–1109.

42.

Tsay,

W.J.

2000 Long memory story of the real interest rate, Economics Letters 67(3), 325–330.

43.

Tudor,

C.A

, Viens,

2007 Statistical aspects of the fractional stochastic calculus, Annals of Statistics 35(3), 1183–1212.

44.

Shea,

G.S.

1991 Uncertainty and implied variance bounds in long memory models of interest rate term structure, Empirical Economics 3(16), 287–312.

45.

Zhang,

, Xiuchun,

, Li,

Rong

, Zhang, S. 2017 Pricing credit derivatives under fractional stochastic interest rate models with jumps, Journal of Systems Science and Complexity 30, 645–659.

46.

Zhu,

2014 Limit theorems for Cox-Ingersoll-Ross process with Hawkes jumps, Journal of Applied Probability 51(3), 699–712.