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

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 ≥ 0 , P ) be a stochastic basis. Let { W t H } t ≥ 0 be a persistent (Hurst parameter H > 0.5) fractional Brownian motion with the filtration { F t } t ≥ 0 and α > 0, β > 0, σ > 0 are the unknown parameters. Recall that a fractional Brownian motion (fBM) has the covariance C ˜ H ( s , t ) = 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 ≥ 0 , X 0 = ξ(1.1) where { W t H } is a fractional Brownian motion with H > 1/2 with the filtration { F t } t ≥ 0 and θ < 0 is the unknown parameter.

Davydov (1970) gave an AR(1) approximation of the fOU model: y j = ρ y j - 1 + v j , ( 1 - L ) H - 1 / 2 v j = ϵ j , ρ = 1 , y 0 = 0 , j = 1 , 2 . … , T where L is the lag operator, ϵ _j ∼ i.i.d. (0, σ²) with E ( ϵ j 4 ) < ∞ , whereas {v _j } is a stationary long-memory process generated by v j = ( 1 - L ) H - 1 / 2 ϵ j = ∑ k = 0 ∞ Γ ( k + H - 1 / 2 ) Γ ( H - 1 / 2 ) Γ ( k + 1 ) ϵ j - k .

Davydov (1970) proved that λ H 1 / 2 σ T H y [ Tt ] → D W t H as T → ∞ . Next we present the martingale approximation to the nonsemimartingale fBm.

Define κ H : = 2 H Γ ( 3 / 2 - H ) Γ ( H + 1 / 2 ) , k H ( t , s ) : = κ H - 1 ( s ( t - s ) ) 1 2 - H λ H : = 2 H Γ ( 3 - 2 H ) Γ ( H + 1 2 ) Γ ( 3 / 2 - H ) , v t ≡ v t H : = λ H - 1 t 2 - 2 H , M t H : = ∫ 0 t k H ( t , s ) dW s H . 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 : = ∫ 0 t K H ( t , s ) dM s H holds for H ∈ (1/2, 1) where K H ( t , s ) : = H ( 2 H - 1 ) ∫ s t r h - 1 2 ( r - s ) H - 3 2 , 0 ≤ s ≤ 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 ∧ t ) 2 - 2 H , Cov ( W t H , M t H ) = λ H 1 / 2 , M t H = b H ∫ 0 t u 1 / 2 - H dW u = D M t H , b H = 2 - 2 H λ H . Define Q t : = d dv t ∫ 0 t k H ( t , s ) X s ds . It is easy to see that Q t = λ H 2 ( 2 - 2 H ) { t 2 H - 1 Z t + ∫ 0 t r 2 H - 1 dZ s } .(1.2) The process X admits the representation X t = ∫ 0 t K H ( t , s ) dZ s . The natural filtration generated by the fundamental semimartingale process Z t = θ ∫ 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 ( θ ) : = dP θ T dP 0 T ( Z 0 T ) = exp { θ ∫ 0 T Q t dZ t - θ 2 2 ∫ 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 ( θ ) : = ∫ 0 T Q t dZ t - θ ∫ 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 Λ T : = log dP θ 1 T dP θ 0 T ( Z 0 T ) = { ( θ 1 - θ 0 ) ∫ 0 T Q t dM t - ( θ 1 - θ 0 ) 2 2 ∫ 0 T Q t 2 dv t } . ( 1.5 ) 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 [ dP θ 1 T dP θ 0 T ( Z 0 T ) ] . We have ∫ 0 T Q t dM t = 1 2 ( U T + 2 θ V T - T )(1.6) where U T : = λ H 2 - 2 H Z T ∫ 0 T t 2 H - 1 dZ s(1.7) and V T : = ∫ 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 ) = { 4 sin π H θ 2 + 2 u e - θ T π T D T H ( θ ; θ 2 + 2 u ) } 1 / 2(1.9) where D T H ( θ ; α ) : = [ α cosh ( α 2 T ) - θ sinh ( α 2 T ) ] 2 J - H ( α 2 T ) J H - 1 ( α 2 T ) - [ α sinh ( α 2 T ) - θ cosh ( α 2 T ) ] 2 J 1 - H ( α 2 T ) J H ( α 2 T ) ( 1.10 ) for α > 0 and J _ν is the modified Bessel function of first kind of order ν.

Remark When H = 1/2, J - 1 / 2 ( z ) = 2 / ( π z ) cos z , J 1 / 2 ( z ) = 2 / ( π z ) sin z .

In this case we have the well-known Cameron-Martin formula: E exp ( - ∫ 0 T W t 2 dt ) = ( cosh 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 ∈ ℂ . 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 ( - θ T 2 ) [ 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

D T H ( γ ; z 2 ) : = [ γ cosh ( γ 2 T ) - ( θ - 2 z 2 ) sinh ( γ 2 T ) ] 2 J - H ( γ 2 T ) J H - 1 ( γ 2 T ) - [ γ sinh ( γ 2 T ) - ( θ - 2 z 2 ) cosh ( γ 2 T ) ] 2 J 1 - H ( γ 2 T ) J H ( γ 2 T ) ( 1.13 )

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 ( T n ) 2 - 2 H [ i 2 - 2 H - ( i - 1 ) 2 - 2 H ] .(1.14)

We have

Q t = d dv t ∫ 0 t k H ( t , s ) X s ds = κ H - 1 d dv t ∫ 0 t s 1 / 2 - H ( t - s ) 1 / 2 - H X s ds = κ H - 1 λ H t 2 H - 1 d dt ∫ 0 t s 1 / 2 - H ( t - s ) 1 / 2 - H X s ds = κ H - 1 λ H t 2 H - 1 ∫ 0 t d dt s 1 / 2 - H ( t - s ) 1 / 2 - H X s ds = κ H - 1 λ H t 2 H - 1 ∫ 0 t s 1 / 2 - H ( t - s ) - 1 / 2 - H X s ds .(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

Q ˜ n = κ H - 1 λ H n 2 H - 1 ∑ j = 0 n - 1 j 1 / 2 - H ( n - j ) - 1 / 2 - H X j .(1.16)

It is easy to show that Q ˜ n → 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

Q ˜ t i = κ H - 1 λ H t i 2 H - 1 ∑ 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 Q ˜ t i → 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.

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 ≥ 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 X t dW t , t ≥ 0 , X 0 = 1(2.1) where {W _t } is a standard Wiener process with the filtration { F t } t ≥ 0 , α > 0 and β < 0 are the unknown parameters.

The true transition density is given by q ( t , x , y , β ) : = − 2 β ( y x ) α − 1 2 e ( 1 2 − α ) β t } 1 − e β t exp [ 2 β ( x + y ) e − β t − 1 ] I ν ( − 2 β x y sinh ( − 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 ( β , α ) : = dP β , α T dP 0 T ( X 0 T ) = exp { ∫ 0 T α + β X t 4 X t dX t - ∫ 0 T ( α + β 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 ( β , α ) : = { ∫ 0 T α + β X t 4 X t dX t - ∫ 0 T ( α + β X t ) 2 8 X t dt } .(2.4)

Let I T : = ∫ 0 T X t dt , J T : = ∫ 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 ∈ ℝ . Then E exp ( iu I T ) = exp ( 2 iu β + γ coth γ T 2 ) [ cosh γ T 2 + β γ sinh γ T 2 ] - 1(2.6) where γ : = (β² - 2iu) ^1/2. Observe that E exp ( iu I T ) = exp ( 2 iu ( e γ T 2 - e - γ T 2 ) e γ T 2 ( γ + β ) + e - γ T 2 ( γ - β ) ) 2 [ e γ T 2 ( γ + β ) + e - γ T 2 ( γ - β ) ] - 1 .(2.7) Further lim T → ∞ E ( e - uI T ) = lim T → ∞ exp ( - ab 2 σ T ( 1 + 4 u σ b 2 - 1 ) ) .(2.8) E exp ( - u X T - v I T ) = ( 2 ρ e ( b - ρ ) T / 2 2 σ u ( 1 - e - ρ T ) + ( ρ - b ) e - ρ T + ( ρ + b ) ) a / σ exp ( u ( ( ρ + b ) e - ρ T + ( ρ - b ) + 2 v ( 1 - e - ρ T ) 2 σ u ( 1 - e - ρ T ) + ( ρ - b ) e - ρ T + ( ρ + b ) )(2.9) where ρ = : b 2 + 4 σ v . Finally E exp ( - u J T ) = Γ ( k + ν 2 + 1 2 ) Γ ( ν + 1 ) ( x α ) - k β ν 2 + 1 2 exp ( b 2 σ [ at - 2 x e bt - 1 ] ) 1 F 1 ( k + ν 2 + 12 , ν + 1 , β )(2.10) where k = : a 2 σ , α = : be bt σ ( e bt - 1 ) , ν = : 1 σ ( α - σ ) 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 ≥ 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 ) : = 1 - e - a ( T - t ) a , A ( t , T ) : = exp ( ( B ( t , T ) - T + t ) ( a 2 b - σ 2 / 2 ) a 2 - σ 2 B ( t , T ) 2 4 a ) .(2.14)

For the CIR model dX t = a ( b - X t ) dt + σ 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 P ˜ ( exp ( - ∫ t T X t dt ) ) = A ( t , T ) e - B ( t , T ) X t where B ( t , T ) : = 2 ( e γ ( T - t ) - 1 ) ( γ + a ) ( e γ ( T - t ) - 1 ) + 2 γ , A ( t , T ) : = ( 2 γ e ( a + γ ) ( T - t ) / 2 ( γ + a ) ( e γ ( T - t ) - 1 ) + 2 γ ) 2 ab / σ 2(2.16) where γ = a 2 + 2 σ 2 and P ˜ is the risk-neutral measure.

Note that B (t, T) can be written as B ( t , T ) = 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 + σ 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 P ˜ ( exp ( - ∫ t T X t dt ) ) = A H ( t , T ) e - B H ( t , T ) X t(2.19) where B H ( t , T ) = 2 ( e γ ( T - t ) - 1 ) ( γ + a ) ( e γ ( T - t ) - 1 ) + 2 γ , A H ( t , T ) = ( 2 γ e ( a + γ ) ( T - t ) / 2 ( γ + a ) ( e γ ( T - t ) - 1 ) + 2 γ ) 2 ab / σ 2(2.20) where γ = a 2 + 2 σ 2 and P ˜ is the risk-neutral measure.

Now consider the fCIR model dX t = a ( b - X t ) dt + σ 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 ∈ [ 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 + σ X t dW t H ,(2.23) let a = -2θ,  b  - 1/(2θ) ,  σ = 2, R T : = ∫ 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 ) = { 4 sin π H θ 2 + 2 u e - θ T π T D T H ( θ ; θ 2 + 2 u ) } 1 / 2(2.24) where D T H ( θ ; α ) : = [ α cosh ( α 2 T ) - θ sinh ( α 2 T ) ] 2 J - H ( α 2 T ) J H - 1 ( α 2 T ) - [ α sinh ( α 2 T ) - θ cosh ( α 2 T ) ] 2 J 1 - H ( α 2 T ) J H ( α 2 T ) ( 2.25 ) 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 ) .

First we consider the modeling of stock price. Let ( Ω , F , { F t } t ≥ 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 ≥ 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 { ( ( μ - 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 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 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 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 ∂ C t ∂ t + rS t ∂ C t ∂ S t + 1 2 σ 2 S t 2 ∂ 2 C t ∂ 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 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 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 : = log ( S K ) + ( r + σ 2 2 ) ( T - t ) σ T - t , d 2 : = log ( S K ) + ( r - σ 2 2 ) ( T - t ) σ T - t = d 1 - σ 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 ≥ 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 - 1 2 σ 2 ) t + σ W t } , t ≥ 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 : = 1 σ P ln LP ( 0 , s ) P ( 0 , T ) K + σ P 2 , σ P : = σ a ( 1 - e - a ( s - T ) ) 1 - e - 2 aT 2 a .(3.12) When a = 0, σ P = σ ( s - T ) 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 ≥ 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 ≥ 0 ,(3.13) where {W _t } is a standard Wiener process with the filtration { F t } t ≥ 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 = T n = Δ , i = 1 , 2 ⋯ , 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 ≤ t ≤ 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 ( θ ) : = dP θ T dP 0 T ( D 0 T ) = exp { θ ∫ 0 T d σ t 2 - θ 2 2 ∫ 0 T σ t 2 dt } .(3.14) Consider the score function, the derivative of the log-likelihood function, which is given by γ T ( θ ) : = ∫ 0 T d σ t 2 - θ ∫ 0 T σ t 2 dt .(3.15)

A solution of the estimating equation γ _T  (θ) =0 provides the maximum likelihood estimate (MLE) θ ˆ T : = σ T 2 - σ 0 2 - T / 2 ∫ 0 T σ t 2 dt .(3.16) The minimum contrast estimate (MCE) is given by θ ˜ T : = - T 2 ∫ 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 + X t S t dW t ,(3.17) d X t = ( α + β X t ) dt + σ 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 : = ∫ 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 : = ∑ i = 1 n Δ s t i - 1 2 .(3.21) It is well known that P - lim n → ∞ 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 θ ˜ T .

Hence from the definition of MCE, the approximate minimum contrast estimate (AMCE) of θ would be θ ˜ n , T : = - 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 log ( A ( t , s ) / K ) B ( T , s ) [ φ + ψ + B ( T , s ) ] ; 4 ab σ 2 , 2 φ 2 X t e γ ( T - t ) φ + ψ + B ( T , s ) ) - KP ( t , T ) χ 2 ( 2 log ( A ( t , s ) / K ) B ( T , s ) [ φ + ψ ] ; 4 ab σ 2 , 2 φ 2 X t e γ ( T - t ) φ + ψ )(3.24) where φ = : 2 γ σ 2 ( e γ ( T - t ) - 1 ) , ψ = : 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 + X t S t dW t(3.26) d X t = ( 1 + 2 θ X t ) dt + 2 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 = T n , i = 1 , 2 ⋯ , n . Let s _t  = log S _t be the log-price process. Then ds t = ( μ - 1 2 X t ) dt + 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 + S t dW t 1 ,(3.29) dV t = κ ( θ - V t ) dt + σ 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 P ˜ , we have dS t = r S t dt + V t S t d W ˜ t 1(3.32) dV t = κ * ( θ * - V t ) dt + σ V t d W ˜ t 2(3.33) d W ˜ t 1 d W ˜ t 2 = ρ dt , κ * = κ + λ , θ * = κ θ κ + λ(3.34) d W ˜ t 1 = dW t 1 + θ t dt , d W ˜ t 2 = dW t 2 + k V t dt(3.35) and the Radon-Nikodym derivative is given by d P ˜ dP = exp { - 1 2 ∫ 0 t ( θ s 2 + k 2 V s ) ds - ∫ 0 t θ s dW s 1 - ∫ 0 t k V s dW s 2 } ,(3.36) where θ t : = ( μ - r ) / V t is the risk premium.

The characteristic function of s _T under P ˜ is given by ψ T ( u ) = E P ˜ ( exp ( i s T u ) ) = A ( u , T ) e B ( u , T )(3.37) where A ( u , T ) : = exp ( iu ( s 0 + rT ) ) , B ( u , T ) : = - ( u 2 + iu ) ( 1 - e γ T ) V 0 2 γ - ( γ - ( κ - ρ σ ui ) ( 1 - e γ T ) - κ θ σ 2 [ 2 log ( 2 γ - ( γ - κ - ρ σ ui ) ( 1 - e γ T ) 2 γ ) + ( γ - κ - ρ σ ui ) T ] ( 3.38 ) with γ = : ( κ - ρ σ ui ) 2 + ( u 2 + iu ) and we choose the principal branch of the square root.

Option price is defined as E 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 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 ) : = 12 + 1 π ∫ 0 ∞ Re ( 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 ) + a σ 2 [ ( b j - ρ σ ui + d ) ( T - t ) - 2 ln ( 1 - ge d ( T - t ) 1 - g ) ] ,(3.43) D ( T - t , u ) : = b j - ρ σ ui + d σ 2 ( 1 - e d ( T - t ) 1 - ge d ( T - t ) ) , g = b j - ρ σ ui + d b j - ρ σ ui - d ,(3.44) d = ( ρ σ ui - b j ) 2 - σ 2 ( 2 ν j ui - u 2 ) , ν 1 : = 12 , ν 2 : = - 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 + σ X t dW t H(3.46) the European call option is given by C t H = P H ( t , s ) χ 2 ( 2 log ( A H ( t , s ) / K ) B H ( T , s ) [ φ + ψ + B H ( T , s ) ] ; 4 ab σ 2 , 2 φ 2 X t e γ ( T - t ) φ + ψ + B H ( T , s ) )(3.47) - K P H ( t , T ) χ 2 ( 2 log ( A H ( t , s ) / K ) B H ( T , s ) [ φ + ψ ] ; 4 ab σ 2 , 2 φ 2 X t e γ ( T - t ) φ + ψ )(3.48) where φ = : 2 γ σ 2 ( e γ ( T - t ) - 1 ) , ψ = : 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 + V t S t dW t 1 ,(3.50) dV t = a ( b - V t ) dt + σ 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 ) = 12 + 1 π ∫ 0 ∞ Re ( 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 ) + a σ 2 [ ( b j - ρ σ ui + d ) ( T - t ) - 2 ln ( 1 - ge d ( T - t ) 1 - g ) ] ,(3.56) D H ( T - t , u ) : = b j - ρ σ ui + d σ 2 ( 1 - e d ( T - t ) 1 - ge d ( T - t ) ) , g : = b j - ρ σ ui + d b j - ρ σ ui - d ,(3.57) d : = ( ρ σ ui - b j ) 2 - σ 2 ( 2 ν j ui - u 2 ) , ν 1 = 12 , ν 2 = - 12 , a : = κ θ , b 1 : = κ + λ - ρ σ , b 2 : = κ + λ .(3.58)