Eric Djeutcha
email: djeutchaeric@yahoo.fr

Jules Sadefo Kamdem
email: jules.sadefo-kamdem@umontpellier.com

Pricing for a vulnerable bull spread options using a mixed modified fractional Hull-White-Vasicek model

Keywords: Bull spread option, Hull-White-Vasicek model, Double Mellin transform. MSC 2010: 60G22, 60G18

In this paper, in order to serve credit risk management, we introduce a pricing model for a vulnerable Bull Spread options in a Mixed Modified Fractional Hull-White-Vasicek stochastic volatility and stochastic interest rate model. We use Milstein scheme to find the sample paths of asset price and its volatility, and the sample paths of interest rates of asset price movement. We use the double Mellin transform to obtain an analytical vulnerable bull spread call option formula and an analytical vulnerable bull spread put option formula under fractional stochastic volatility and fractional stochastic interest rates.

Introduction

Credit risk management is crucial for financial institutions (banks, insurance companies, mutual funds and pension funds, etc.), but also for national or international regulatory institutions. Indeed, within the framework of market surveillance and the risk of bankruptcy of financial institutions, the bank for international settlements requires compliance (within the framework of the application of the Basel III or Solvency II systems), that certain prudential rules are respected. Thus, hedging financial risks through the use of financial derivatives is increasingly common. It is enough to see the extension works of the option pricing models proposed by Black and Scholes [START_REF] Black | The pricing of options and corporate liabilities[END_REF], to understand the growing interest in the valuation of derivative products. With regard to the risks of credit default on the over-the-counter markets, some authors have, under certain assumptions, proposed pricing models of vulnerable options. We must mention the recent occurrence of economic crises (for example the subprime crisis in 2008), with, as a corollary, the increase in the probability of default by certain financial institutions, companies or municipalities. Among the authors who have proposed pricing models for vulnerable options, the holder of which is subject to the risk of default, in addition to the work of Johnson and Stuly [START_REF] Johnson | The pricing of options with default risk[END_REF] where the first to substituted default risk in option pricing and put forward a new definition called the vulnerable option, Klein [START_REF] Klein | Pricing Black-Scholes options with correlated credit risk[END_REF] has derived an analytical pricing formula which not only allows a correlation between the underlying asset of the option and the default risk of the counterparty, but also for the writer of the option to have other liabilities in the capital structure. Ammann [START_REF] Ammann | Credit Risk Valuation: Methods, Models, and Applications[END_REF] developed Klein's credit risk model on the basis of a structural approach. He finally obtained the explicit expression for the vulnerable option under the assumption of interest rate and default intensity obeying a stochastic differential equation.

Other scholars like Chang and Hung [START_REF] Chang | Valuation of vulnerable American options with correlated credit risk[END_REF] have also discussed this problem, while all the above mentioned discussions are in the environment of geometric Brownian motion. Due to the inadequacies of geometric Brownian motion in describing self-similarity and long-term dependence on stock prices, fractional Brownian motion and mixed fractional Brownian motion are widely used in the valuation of assets. Hu and Oksendal [START_REF] Hu | Fractional white noise calculus and applications to finance, Infinite Dimensional Analysis[END_REF], Djeutcha et al. [START_REF] Djeutcha | Pricing for Options in a Hull-White-Vasicek volatility and Interest rate Model[END_REF] developed the structural approach, respectively, provided that the stock prices followed a mixed modified fractional Brownian motion and mixed fractional Brownian motion and they proved that the correspondence with the fractional market and the mixed modified fractional Black-Scholes market had no arbitrage effect for any Hurst parameter greater than 1/2.

For more literature on fractional Brownian motions and mixed modified fractional Brownian motions, among others , we can refer to Oksendal [START_REF] Ãksendal | Fractional Brownian Motion in Finance[END_REF], Djeutcha and Sadefo Kamdem [START_REF] Djeutcha | Local and implied volatilities with the mixed-modifiedfractional-Dupire model[END_REF] and Djeutcha et al. [START_REF] Djeutcha | Pricing for Options in a Hull-White-Vasicek volatility and Interest rate Model[END_REF]. But there is another problem that the Black-Scholes fractional market does not have an equivalent martingale measure according to Sottinen and Valkeila [START_REF] Sottinen | Fractional Brownian motion as a model in finance[END_REF]. Necula [START_REF] Necula | Option Pricing in a Fractional Brownian Motion Environment[END_REF] applied the quasi-martingale method to risk-neutral measurement. Huang et al. [START_REF] Huang | Pricing formulae for European options under the fractional Vasicek interest rate model[END_REF] obtained the explicit expression for the price of the European option under the hypothesis of a Black-Scholes fractional market. Su and Wang [START_REF] Su | Pricing options with credit risk in a reduced form model[END_REF] and Li and Ma [START_REF] Li | Pricing options with credit risk in Markovian regime-switching markets[END_REF] derived the closed form formula of the price of the vulnerable European option by the method of changing measures.

The constant volatility assumption of the Black-Scholes [START_REF] Black | The pricing of options and corporate liabilities[END_REF] model is unreliable, because in the market volatility is in the form of a smile. Hence the need to introduce in the market stochastic volatility models which give more realistic results in relation with the historical asset returns and interest rate variation, compared to the classical Black-Scholes model. Pioneers in the field are: Hull-White [START_REF] Hull | The pricing of options with stochastic volatilities[END_REF], Stein-Stein [START_REF] Stein | Stock price distributions with stochastic volatility: An analytic approach[END_REF], Heston [START_REF] Heston | A closed-form solution for options wth stochastic volatility with applications to bond and currency options[END_REF], Njamen et al. [START_REF] Njomen | Existence and Uniqueness Solution under Non-lipschiz condition of the mixed fractional Heston's model[END_REF], Bakshi, Cao-Chen [START_REF] Bakshi | Pricing and hedging long-term options[END_REF] and Schobel-Zhu [START_REF] Schöbel | Stochastic volatility with an Ornstein-Uhlenbeck process: an extension[END_REF]. Each of these authors model the volatility of the underlying and propose models that present analytical solutions and closed formulas for option pricing.

We therefore appeal to the class of stochastic volatility models and the pioneers in this field are among these stochastic volatility models, the best known and studied is that of Hull-White. But, if in the Hull-White model we consider the interest rate to be stochastic and follow the Vasicek model, we can define a new kind of three-factor model. In general, these classes of three-factor models are not studied in the literature, at least there is a lack of method to determine the analytical explicit expression for the option hedging problem.

These competing models modeled the interest rate of the Black and Scholes model and authors such as Vasicek [START_REF] Vasicek | An Equilibrium Characterisation of the Term Structure[END_REF], Dothan [START_REF] Dothan | On the term structure of interest rates[END_REF], Cox-Ingersoll-Ross [START_REF] Cox | A theory of the term structure of interest rates[END_REF], Djeutcha et al. [START_REF] Djeutcha | Pricing for Options in a Hull-White-Vasicek volatility and Interest rate Model[END_REF], Yousouf et al. [START_REF] Yousuf | Partial differential integral equation model for pricing American option under multi state regime switching with jumps[END_REF], Mohamadinejad et al. [START_REF] Mohamadinejad | Spread option pricing using two jump-diffusion interest rates[END_REF] and Black-Derman-Toy [START_REF] Black | A one-factor model of interest rates and its applications[END_REF] have worked on these models.

In this article, we will be using the double Mellin transform method as described in Ji-Hun and Jeong-Hoon [START_REF] Ji-Hun | The pricing of vulnerable options with double Mellin transforms[END_REF], so that we can derive the general pricing formula for the vulnerable bull spread option under the stock price assumption, with the stock's volatility obeying the Hull and White diffusion model and the underlying interest rate obeying the Vasicek model which are driven by the mixed modified fractional Brownian motion.

The competition over bull spread options pricing is more intense comparison with vulnerable types. This is due to the uncertainty of the execution time in bull spread option and in result. The vulnerable option is a kind of option with credit risk which refers to a risk that a borrower will default on any type of debt by failing to make the required payments.

The bull spread option is an optimistic option strategy designed to take advantage of a moderate rise in the price of a security or asset. If we replace in the European vulnerable option, the payoff function of the European option with the payoff function of the bull spread option, we get a new option called the vulnerable bull spread option.

Our paper looks at these two families of models (stochastic interest rate) and stochastic volatility). But our motivation comes from the fact that we model the two families of models by assuming that their random parts are described by mixed modified fractional Brownian motions in such a way that our market is without arbitrage. The combination of the stochastic interest rate model (Vasicek[33]) and the stochastic volatility model (Hull-White [START_REF] Hull | The pricing of options with stochastic volatilities[END_REF]) allows us to define a three-factors model. This new model has a double advantage in that, it can correctly price any bond unlike the one-factor models mentioned above and can allow any options to be valued. In this paper, we will therefore be talking about simulating of our model using Milstein's discretization approach and proposing a closed formula for the pricing of the particular case of vulnerable options and the bull spread option using the double Mellin's transform method.

Nammed MMFHWV model for "Mixed-Modified-Fractional-Hull-White-Vasiseck ", this new model is a combination of the Hull and White [START_REF] Hull | The pricing of options with stochastic volatilities[END_REF] and the Vasicek [START_REF] Vasicek | An Equilibrium Characterisation of the Term Structure[END_REF] model. In this model, the volatility process and asset model are not correlated, while the interest rate process and asset model, the volatility process and interest rates process are correlated, with each other and they are controlled by a distinct diffusion process. In MMFHWV model, the existence of the mean reversion process causes the adjustment of the volatility process and the interest rates behavior in the financial markets and it is a benefit of the MMFHWV model. The rest of this paper is organized as follows. Section 2 gives the preliminaries and basic assumptions. In Section 3, we introduce the MMFHWV model. We simulate the MMFHWV Model based on Milstein discretization. In section 5, we study the pricing of a vulnerable bull spread option under the MMFHWV model and conclusions are presented in the last section.

Preliminaries and Basic Assumptions

Preliminaries

In this subsection, we shall briefly review the definition and some main properties of the mixed modified fractional Brownian motion, the double Mellin transform, we also recall the Hull-White model and the Vasicek model. Fractional Brownian motion or Mellin transform have also been used for option pricing in Djeutcha et al. [START_REF] Djeutcha | Solving Arbitrage Problem on the Financial Market Under the Mixed Fractional Brownian Motion With Hurst Parameter H ∈]1/2, 3/4[END_REF], Thao [START_REF] Thao | An approximate approach to fractional analysis for finance[END_REF], Sadefo Kamdem [START_REF] Sadefo Kamdem | Option pricing with Levy process using Mellin Transform[END_REF],Ji-Hun [START_REF] Ji-Hun | The pricing of vulnerable options with double Mellin transforms[END_REF],Hull-White [START_REF] Hull | The pricing of options with stochastic volatilities[END_REF] and Vasicek [START_REF] Vasicek | An Equilibrium Characterisation of the Term Structure[END_REF].

Definition 2.1. The double Mellin transform is defined by

M xy (f (x, y), w 1 , w 2 ) = f (w 1 , w 2 ) = ∞ 0 ∞ 0 f (x, y)x w 1 -1 y w 2 -1 dxdy (1)
where f (x, y) is a locally Lebesgue integrable function and w 1 and w 2 are complex numbers. Also, if a < Re(w 1 ), Re(w 2 ) < b and if c 1 and c 2 are such that a < c 1 < b and a < c 2 < b, then the inverse of the double Mellin transform is given by

f (x, y) = M -1 xy ( f (w 1 , w 2 )) = 1 (2π) 2 a 1 +i∞ a 1 -i∞ a 2 +i∞ a 2 -i∞ f (w 1 , w 2 )x -w 1 y -w 2 dw 1 dw 2 (2)
We have the following properties for derivative

M xy (x 2 ∂f (x, y) ∂x 2 , w 1 , w 2 ) = w 1 (w 1 + 1) f (w 1 , w 2 ) (3)
M xy (y 2 ∂f (x, y) ∂y 2 , w 1 , w 2 ) = w 2 (w 2 + 1) f (w 1 , w 2 ) (4)
Lemma 2.1. Ji-Hun and Jeong-Hoon [START_REF] Ji-Hun | The pricing of vulnerable options with double Mellin transforms[END_REF]Given complex numbers α and β with Re(α) ⩾ 0, let

f (x) = 1 2πi a-i∞ a-i∞ f (w)x -w dw, where f (w)x -w = e α(w+β) 2 . Then f (x) = 1 2 1 √ πα x β e 1 4α (ln x) 2 holds.
Lemma 2.2. Ji-Hun and Jeong-Hoon [START_REF] Ji-Hun | The pricing of vulnerable options with double Mellin transforms[END_REF] Let f and g be functions from R n + into C. If f (w 1 , w 2 ) and g(w 1 , w 2 ) are the double Mellin transform of f (x, y) and g(x, y), respectively, given by

■(w 1 , w 2 ) = ∞ 0 ∞ 0 ■(x, y)x w 1 -1 y w 2 -1 dxdy. ( 5
)
Then the double Mellin convolution of f and g is given by inverse double Mellin transform of f (w 1 , w 2 ) g(w 1 , w 2 ) as follows

f (x, y) ⋆ g(x, y) = M -1 w 1 w 2 f (w 1 , w 2 ) g(w 1 , w 2 ); x, y = 1 (2π) 2 ∞ 0 ∞ 0 f x u , y w g(u, w)u -1 w -1 dudw (6)
Definition 2.2. Hull and White [START_REF] Hull | The pricing of options with stochastic volatilities[END_REF]The Hull-White model assumes that the underlying stock price, S t and variance V t follow a Black-Scholes-type stochastic process. Hence, the Hull-White model is represented by the bivariate system of stochastic differential equations (SDEs):

dS t = µS t dt + S t σ s dB t,s , dV t = µV t dt + σ v V t dB t,v . (7)
Where E P [dB

dr t = κ(θ -r t )dt + r t σ r dB t,r , (8)
where B .,r is the wiener process.

The parameters of the model are described in Appendix A.1.

Definition 2.4. Djeutcha et al. [START_REF] Djeutcha | Solving Arbitrage Problem on the Financial Market Under the Mixed Fractional Brownian Motion With Hurst Parameter H ∈]1/2, 3/4[END_REF][START_REF] Djeutcha | Local and implied volatilities with the mixed-modifiedfractional-Dupire model[END_REF]A Mixed Modified Fractional Brownian Motion process whose parameters a, b, ε and H is a linear combination of Brownian Motion B t and independent semimartingale process B H,ε t , defined on a probability space (Ω, F, P) by:

M H,a,b,ε t = M H,ε t = aB t + bB H,ε t , ∀(a, b) ∈ R * + × R + , t ∈ [0, T ]. (9)
Where H ∈] 1 2 , 1[ is the Hurst parameter.

In this paper, our aim is to establish an extension of Hull-White model Eq.( 7) and Vasicek model Eq. [START_REF] Dothan | On the term structure of interest rates[END_REF], by exchanging in Eq. [START_REF] Cox | A theory of the term structure of interest rates[END_REF] and Eq.( 8) the Brownian motion B .,s , B .,v and B .,r by the Mixed Modified Fractional Brownian process defined by the following system:

     M H 1 ,ε t,s = aB t,s (t) + bB H 1 ,ε t,s M H 2 ,ε t,v = aB t,v (t) + bB H 2 ,ε t,v M H 3 ,ε t,r = aB t,r (t) + bB H 3 ,ε t,r . (10)
We suppose that: B .,s &B H 1 ,B .,v &B H 2 and B .,r &B H 3 are independents.

From Thao [START_REF] Thao | An approximate approach to fractional analysis for finance[END_REF], we also deduce that

B H 1 ,ε s , B H 2 ,ε v
and B H 3 ,ε r are three local martingale defined as follows:

     B H 1 ,ε t,s = t 0 (t + ε -s * ) H 1 -1 2 dB s * ,s B H 2 ,ε t,v = t 0 (t + ε -s * ) H 2 -1 2 dB s * ,v B H 3 ,ε t,r = t 0 (t + ε -s * ) H 3 -1 2 dB s * ,r . (11)
From Eq.( 10), we define

     dM H 1 ,ε t,s = adB t,s (t) + bdB H 1 ,ε t,s dM H 2 ,ε t,v = adB t,v (t) + bdB H 2 ,ε t,v dM H 3 ,ε t,r = adB t,r (t) + bdB H 3 ,ε t,r . (12)
We know that

         dB H 1 ,ε t,s = ϕ ε t,s dt + a + bε H 1 -1 2 dB t,s dB H 2 ,ε t,v = ϕ ε t,v dt + a + bε H 2 -1 2 dB t,v dB H 3 ,ε t,r = ϕ ε t,r dt + a + bε H 3 -1 2 dB t,r . (13)
where

     ϕ ε t,s = t 0 (H 1 -1 2 )(t + ε -s * ) H 1 -3 2 dB s * ,s , ϕ ε t,v = t 0 (H 2 -1 2 )(t + ε -s * ) H 2 -3 2 dB s * ,v , ϕ ε t,r = t 0 (H 3 -1 2 )(t + ε -s * ) H 3 -3 2 dB s * ,r . (14)
By Letting χ

H i ,ε = a + bε H i -1 2 2
, i = 1, 2, 3. Eq.( 12) helps to define the following expressions:

     dM H 1 ,ε t,s = ϕ ε t,s dt + χ H 1 ,ε dB s * ,s , dM H 2 ,ε t,v = ϕ ε t,v dt + χ H 2 ,ε dB s * ,v , dM H 3 ,ε t,r = ϕ ε t,r dt + χ H 3 ,ε dB s * ,r . (15)

Basic assumptions

In this subsection, we introduce the basic assumptions. Some of them are those of Djeutcha et al. [START_REF] Djeutcha | Solving Arbitrage Problem on the Financial Market Under the Mixed Fractional Brownian Motion With Hurst Parameter H ∈]1/2, 3/4[END_REF].

(a) Consider the probability space (Ω, F, P) where F = (F t ) t⩾0 is a complete and right continuous filtration generated by

F t = σ M H 1 ,ε t,s ; M H 2 ,ε t,s ; M H 3 ,ε t,s
and P is a probability measure on Ω.

(b) The dynamic of stock price S t and its volatility V t are given by a Mixed Modified Fractional Hull-WHite no-arbitrage model, namely :

dS t = µS t dt + σ s S t dM H 1 ,ε t,s , dV t = µV t dt + σ v V t dM H 2 ,ε t,v . (16)
where M H is given by the following matrix defined from Eq.(15) as follows:

E P [dM H i ,ε t,s dM H j ,ε t,s ] =          χ H i ,ε dt, if i = j ρ sv χ H i ,ε χ H j ,ε dt, if i ̸ = j, i, j = 1, 2. ρ sr χ H i ,ε χ H j ,ε dt, if i ̸ = j, i, j = 1, 3. ρ vr χ H i ,ε χ H j ,ε dt, if i ̸ = j, i, j = 2, 3. (18)
Where

χ H i ,ε χ H j ,ε = a 2 + abε H i -1 2 + abε H j -1 2 + b 2 ε H i +H j -1 , ∀i ̸ = j. (19)
(f) The correlation coefficients ρ sv , ρ sr and ρ vr are the parameters satisfying

ρ 2 sv < 1; ρ 2 sr < 1; ρ 2 vr < 1; ρ 2 sv + ρ 2 sr + ρ 2 vr -2ρ sv ρ sr ρ vr < 1. ( 20
)
3 The MMFHWV Model Framework

In the Mixed Modified Fractional Hull-White model Eq.( 16), if we select the parameter µ as a stochastic process(see Mixed Modified Fractional Vasicek model( 16)), we obtain the Mixed Modified Fractional Hull-White-Vasicek(MMFHWV) model defined system of trivariate fractional stochastic differential equations (SDEs):

     dS t = r t S t dt + σ s S t dM H 1 ,ε t,s , dV t = r t V t dt + σ v V t dM H 2 ,ε t,v , dr t = κ(θ -r t )dt + r t σ r dM H 3 ,ε t,r . (21)
Where M H 1 ,ε t,s , M H 2 ,ε t,v , and M H 3 ,ε t,r are three correlated standard Mixed Modified Fractional Brownian motions.

The stock price and variance follow the processes in Eq.( 16) and the interest rate follow the processes in Eq.( 17) under the historical measure P also called the physical measure. For pricing purposes, however, we need the processes for (S t , V t , r t ) under the risk-neutral measure Q. In the MMFHWV model Eq.( 21), this is done by modifying each SDE in Eq.( 16) and Eq.( 17) separately by an application of Girsanov's theorem. The risk-neutral process for the MMFHWV model Eq.( 21) is defined by

     dS t = r t S t dt + σ s S t d M H 1 ,ε t,s , dV t = r t V t dt + σ v V t d M H 2 ,ε t,v , dr t = κ * (θ * -r t )dt + r t σ r d M H 3 ,ε t,r . (22)
With

       M H 1 ,ε t,s = M H 1 ,ε t,s + µ-r σs t M H 2 ,ε t,v = M H 2 ,ε t,v + µ-r σv t M H 3 ,ε t,r = M H 3 ,ε t,r + µ-r σr √ χ H 3 ,ε t. (23)
Where M H 1 ,ε t,s , M H 2 ,ε t,v , and M H 3 ,ε t,r are three correlated standard Mixed Modified Fractional Brownian motions. Its correlation matrix is defined from Eq.( 15) as follows:

E Q [d M H i ,ε t,s d M H j ,ε t,s ] =          χ H i ,ε dt, if i = j ρ sv χ H i ,ε χ H j ,ε dt, if i ̸ = j, i, j = 1, 2. ρ sr χ H i ,ε χ H j ,ε dt, if i ̸ = j, i, j = 1, 3. ρ vr χ H i ,ε χ H j ,ε dt, if i ̸ = j, i, j = 2, 3. (24)
Where the expression of χ H i ,ε χ H j ,ε is given by Eq.( 19), and Q is the risk-neutral measure. Using Eq.( 15), the risk-neutral process for the stock price is

dS t = S t r + bσ s ϕ ε t,s dt + σ s χ H 1 ,ε d B t,s . (25)
where

B t,s = B t,s + µ -r σ s χ H 1 ,ε t (26)
The risk-neutral process for the variance of stock price is

dV t = V t r + bσ v ϕ ε t,v dt + σ v χ H 2 ,ε d B t,v . (27)
where

B t,v = B t,v + µ -r σ v χ H 2 ,ε t (28)
The risk-neutral process for the interest rate is

dr t = κ * (θ * -r t ) + σ r ϕ ε t,r dt + σ r χ H 3 ,ε d B t,r . (29)
where

B t,r = B t,r + λr t σ r χ H 3 ,ε t (30)
with

κ * = κ + λ σr √ χ H 3 ,ε and θ * = κθ κ+ λ σr √ χ H 3 ,ε
, are the risk-neutral parameter and λ the interest rate risk parameter which is given by Breeden [START_REF] Breeden | An Intertemporal Asset Pricing Model With Stochastic Consumption and Investment Opportunities[END_REF].

Using Eq.( 25), Eq.( 27) and Eq.( 29), we propose the new model appointed Mixed-Modified-Fractional-Hull-White-Vasicek(MMFHWV) model defined by the trivariate system of stochastic diferential equations(SDEs) given in the following definition: 25) and Eq.( 27), if we select the parameter µ as a stochastic process( see Mixed Modified Fractional Vasicek model Eq.( 29)), then we will obtain the new model appointed Mixed-Modified-Fractional-Hull-White-Vasicek (MMFHWV) model and then, under the risk-neutral measure Q, the dynamics of S t , V t and r t are given by the SDEs:

       dS t = S t r + bσ s ϕ ε t,s dt + σ s χ H 1 ,ε d B t,s , dV t = V t r + bσ v ϕ ε t,v dt + σ v χ H 2 ,ε d B t,v , dr t = κ * (θ * -r t ) + σ r ϕ ε t,r dt + σ r χ H 3 ,ε d B t,r . (31)
where B t,s , B t,v and B t,r are three Brownian motions defined by Eq.( 26), Eq.( 28) and Eq.( 30) respectively; with

κ * = κ + λ σr √ χ H 3 ,ε , θ * = κθ κ+ λ σr √ χ H 3 ,ε
, are the risk-neutral parameters and λ the interest rate risk parameter.

Simulation of MMFHWV Model

We often face some difficult problems which can't be solved analytically. In these cases, simulation methods are appropriate and really helpful. In order to simulate the asset price S t at time t, the volatility of asset price V t at time t and the interest rate denote by r t at time t, we need sample paths of asset price movement, sample paths of volatility of asset price movement and sample paths of interest rates evolution.

The Milstein discretization can be used to approximate the asset path of the asset, the volatility path of the volatility and the interest path of the interest rate on a discrete time grid, see Glasserman [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF], Kloeden and Platen [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] and Kahl and Jackel [START_REF] Kahl | Fast Strong Approximation Monte-Carlo Schemes for Stochastic Volatility Models[END_REF].

Let Υ = {t 1 , t 2 , • • • , t N } be a partition on time interval [0, T ]; i.e. 0 = t 0 < t 1 < • • • < t N = T.
Then, we have

         S j+1 = S j + r j + bσ s ϕ ε j,s S j ∆ t + σ s χ H 1 ,ε ∆B i,j V j+1 = V j + r j + bσ v ϕ ε j,v V j ∆ t + σ v χ H 2 ,ε ∆B i,j r j+1 = r j+1 + κ * (θ * -r j ) + σ r ϕ ε j,r ∆ t + σ r χ H 3 ,ε ∆B i,j , (32)
where B i,j is a Brownian motion and we have the following formula:

B i,j = B i (t j+1 ) -B i (t j ), 0 ⩽ j ⩽ N -1 and i = 1, 2, 3.
Moreover, ∆B i,j ∼ N (0, ∆ t ) and t j = j∆ t .

According to the Central Limit Theorem (CLT), we have ∆B i,j = Z i √ ∆ t , Z i ∼ N (0, 1) and we obtain from Eq.( 32)

         S j+1 = S j + r t + bσ s ϕ ε j,s S j ∆ t + S j σ s χ H 1 ,ε √ ∆ t Z 1 + 1 2 σ 2 s χ H 1 ,ε S j ∆ t Z 2 1 -1 V j+1 = V j + r j + bσ s ϕ ε j,s V j ∆ t + V j σ s χ H 2 ,ε √ ∆ t Ψ 1 + 1 2 σ 2 v χ H 2 ,ε V j ∆ t Ψ 2 1 -1 r j+1 = r j + κ * (θ * -r j ) + bσ r ϕ ε j,r dt + σ r χ H 3 ,ε √ ∆ t Ψ 2 . (33)
Where

Ψ 1 = ρ sv Z 1 + 1 -ρ 2 sv Z 2 and Ψ 2 = ρ sr Z 1 + 1 -ρ 2
sr Z 3 are two correlated normal variables generated by the Cholesky decomposition. We have

         S j+1 = S j + r t + bσ s ϕ ε j,s S j ∆ t + S j σ s χ H 1 ,ε √ ∆ t Z 1 + 1 2 σ 2 s χ H 1 ,ε S j ∆ t Z 2 1 -1 , V j+1 = V j + r j + bσ s ϕ ε j,s V j ∆ t + V j σ s χ H 2 ,ε √ ∆ t Ψ 1 + 1 2 σ 2 v χ H 2 ,ε V j ∆ t Ψ 2 1 -1 , r j+1 = r j + κ * (θ * -r j ) + bσ r ϕ ε j,r ∆ t + σ r χ H 3 ,ε √ ∆ t Ψ 2 . ( 34
)
The simulation of the MMFHWV model is given by the following algorithm.

Algorithm. MMFHWV model simulation.

1. Set ∆t = t N . 2. For i = 1 to number of simulation.

3. Generate independent standard normal variables, Z j ∼ N (0, 1), j = 1, . . . , N .

Set

V j+1 = V j + r j + bσ s ϕ ε j,s V j ∆ t + V j σ s χ H2,ε √ ∆ t Ψ 1 + 1 2 σ 2 v χ H2,ε V j ∆ t Ψ 2 1 -1 . 5. Set S j+1 = S j + r t + bσ s ϕ ε j,s S j ∆ t + S j σ s χ H1,ε √ ∆ t Z 1 + 1 2 σ 2 s χ H1,ε S j ∆ t Z 2 1 -1 . 6. Set r j+1 = r j + κ * (θ * -r j ) + bσ r ϕ ε j,r ∆ t + σ r χ H3,ε √ ∆ t Ψ 2 .
7. End For.

In fig. 1, 2 and 3, we illustrate changes in the value of stock path, volatility of the stock's path and interest rate path under modified Hurst parameter. The result indicate that:

• By increasing the value of H 1 ,(see Fig. 1), the value of the stock price is reduced • By increasing the value of H 2 ,(see Fig. 2),the value of the volatility decrease and more, quickly convergence to zero.

• By increasing the value of H 3 ,(see Fig. 3, the value of the Interest rate increase. 5 Value of the call option under the MMFHWV model using the double Mellin transform

Bond price formula

In this subsection,we first give the pricing formula of a zero coupon bond in the mixed modified fractional Hull-White Vasicek interest rate model. Let P (t, r, T ) denote the price of a default-free zero-coupon bond with the nominal value 1 unit of money at time T. By the Feynman-Kac formula, we have the following pricing formula.

Theorem 5.1. In the fixed MMFHWV model, the price of a zero-coupon bond at time t is given by:

P (τ, r, T ) = A(τ )e -B(τ )τ , (35)
with τ = T -t.A(τ ) and B(τ ) are given by

   A(τ ) = exp χ H 3 ,ε σ 2 r 2κ 2 (τ+ 2 κ * (e -κ * τ -1)-1 2κ * (e -κ * τ -1)) - τ 0 κ * θ * (T -z)B(z)dz B(τ ) = 1-e -κ * τ κ * (36)
Proof of Theorem 5.1. We recall that the no-arbitrage price at time t of the bond with the nominal value 1 is given by

P (t, r, T ) = E e -T t rsds |r t = r) , t < T. ( 37
)
with the final condition P (T, r, T ) = 1. Then by using the Mixed modified Fractional Feynman-Kac formula, the solution of P (t, r, T ) satisfies the following PDE problem

∂P ∂t + 1 2 χ H3,ε σ 2 r ∂ 2 P ∂r 2 + κ * (θ * -r) ∂P ∂r -rP = 0. ( 38
)
The affine term structure solution of Eq.( 38) is given by

P (τ, r, T ) = A(τ )e -B(τ )τ , τ = T -t (39)
with If we substitute Eq.(39) into Eq.(38), we derive the following ordinary differential equation:

dA(τ ) dτ -κ * θ * B(τ ) + 1 2 χ H3,ε σ 2 r B 2 (τ ) = 0 (40) and dB(τ ) dτ -κ * B(τ ) + 1 = 0 (41)
The resolution of Eq.(41) gives

B(τ ) = 1 -e -κ * τ κ * (42)
and the resolution of Eq.( 40) can be expressed in terms of B(τ ) to give the value of A(τ )

Pricing for Bull Spread Option under MMFHWV model

We recall that, if we use the traditional probabilistic method such as the Monte Carlo method to solve system given by Eq.( 31), for Bull Spread options, that method produce interesting result, but it could requires a long time, depending on the complexity of calculations. However, in this section, using the double Mellin transform, we obtain a closed form formula for the option price, that is credible for both the accuracy and efficiency.

The bull spread option is an optimistic option strategy designed to take advantage of a moderate rise in the price of a security or asset. Its payoff function is defined as follows:

l(S T ) = (S T -K 1 ) + -(S T -K 2 ) + . ( 43
)
The vulnerable option is a kind of option with credit risk which refers to a risk that a borrower will default on any type of debt by failing to make the required payments.Its payoff function is defined as follows:

m(S T ) = (S T -K) + 1 {V T ⩾γ * } + 1 {V T <γ * } (1 -α)V T γ , (44)
where α denotes the dead weight costs associated with bankruptcy expressed as a percentage of the value of the assets of the counterparty, the constant parameter D denotes the default boundary. The parameter γ * may be less than D due to the possibility of a counterparty continuing its operation even while V T < γ. In the event of a credit loss, only the proportion (1-α)V T γ of the nominal claim is paid out by the counterparty. As aforementioned case for vulnerable bull spread option, expiring at time T , with the strike price K 1 and K 2 , the payoff function is given by

h(S T , V T )) = (S T -K 1 ) + -(S T -K 2 ) + 1 {V T ⩾γ * } + 1 {V T <γ * } (1 -α)V T γ , (45)
with γ * is a critical value such that a credit loss occurs if the value of the option writer's assets and the no-arbitrage price of the option is given by

U (t, s, v, r) = E * e -T t rsds h(S T , V T )|S t = s, V t = v, r t = r) . (46)
By the non-arbitrage assumptions, the corresponding MMFHWV PDE of U (t, s, v, r) according to the Feynman-Kac formula of bull spread option.

In this paper, we consider U = U (S, V, r, T -t) as the common arguments in financial mathematics imply that if the option value function U is sufficiently smooth then under the MMFHWV model, it satisfies the partial differential equation given in the following theorem.

Theorem 5.2. Under the MMFHWV model, the option value function U satisfies in the domain:

Ω = {(S, V, r, t) : 0 < S < +∞, 0 < V < +∞, -∞ < r < +∞, 0 < t < T } , (47)
the following partial differential equation is given by:

∂U ∂t + 1 2 χ H1,ε σ 2 s s 2 ∂ 2 U ∂s 2 + 1 2 χ H2,ε σ 2 v v 2 t ∂ 2 U ∂v 2 + 1 2 χ H3,ε σ 2 r ∂ 2 U ∂r 2 + ρ sv σ s σ v sv χ H1,ε χ H2,ε ∂ 2 U ∂s∂v + ρ sr σ s σ r s χ H1,ε χ H3,ε ∂ 2 U ∂s∂r + ρ vr σ v σ r v χ H2,ε χ H3,ε ∂ 2 U ∂v∂r + r s ∂U ∂s + v ∂U ∂v + κ * (θ * -r) ∂U ∂r -rU = 0 (48)
with the terminal condition U (T, s, v, r) = h(s, v) and χ Hi,ε = (a + bε Hi-1 2 ) 2 ,i=1.2.3.

Proof of Theorem 5.2. The proof of this result is found in Appendix A.2

Option price formula

In this subsection, we use the approach of Ji-Hun and Jeong-Hoon [START_REF] Ji-Hun | The pricing of vulnerable options with double Mellin transforms[END_REF] to derive a price formula of a vulnerable bull Spread option under the Mixed-Modified-Fractional-Hull-White-Vasicek model by using the double mellin transform.

We observe that the payoff function h of a vulnerable bull spread call option is not bounded, we require to modify h by define a sequence of function h n such that each h n is bounded and lim n-→∞ h n (s, v, r) = h(s, v, r). They are given by

             h n (s, v, r) = h 1 n (s, v, r)h 2 n (s, v, r) h 1 n (s, v, r) = (s -K 1 ) -(s -K 2 ), if K 2 ⩽ K 1 < s < n 0, otherwise h 2 n (s, v, r) = 1, if γ * ⩽ v < n (1-α)v γ , otherwise (49)
Setting U n (T, s, v, r) by

U n (T, s, v, r) = E * e -T t rsds h(S T , V T )|S t = s, V t = v, r t = r) . (50)
We observe that, U n (t, s, v, r) satisfies Eq.( 48) with the terminal condition U n (T, s, v, r) = h n (s, v).

The price U (t, s, v, r) is given by taking the limit

U (t, s, v, r) = lim n-→∞ U n (t, s, v, r) Let U n (t, w 1 , w 2 , r)
be the double Mellin transform of U n (t, s, v, r). Using Eq.( 2),the expression of the inverse double Mellin transform of U n (t, s, v, r) is given by

U n (t, s, v, r) = 1 (2π) 2 a1+i∞ a1-i∞ a2+i∞ a2-i∞ U n (t, w 1 , w 2 , r)s -w1 v -w2 dw 1 dw 2 (51)
By applying Eq.(51) into Eq.( 48), we obtain

- ∂ U n ∂τ + 1 2 χ H1,ε σ 2 s w 1 (w 1 + 1) + ρ sv χ H1,ε χ H2,ε σ v σ s w 1 w 2 + 1 2 χ H2,ε σ 2 v w 2 (w 2 + 1) -r(w 1 + w 2 + 1) U n + 1 2 χ H3,ε σ 2 r ∂ 2 U n ∂r 2 + κ * θ * -κ * r -ρ sr χ H1,ε χ H3,ε σ s σ r w 1 -ρ vr χ H2,ε χ H3,ε σ r w 2 ∂ U n ∂r = 0 (52)
with τ = T -t and U n (0, w 1 , w 2 , r) = h(w 1 , w 2 ) which is the Mellin transform of h n (s, v). Let

U n (τ, w 1 , w 2 , r) = exp 1 2 χ H1,ε σ 2 s w 1 (w 1 + 1) + ρ sv χ H1,ε χ H2,ε σ v σ s w 1 w 2 + 1 2 χ H2,ε σ 2 v w 2 (w 2 + 1) τ × f n (τ, w 1 , w 2 , r). (53)
Then, Eq.( 52) is transformed into the following PDE for f n .

-

∂ f n ∂τ + 1 2 χ H3,ε σ 2 r ∂ 2 f n ∂r 2 + κ * θ * -κ * r -ρ sr χ H1,ε χ H3,ε σ s σ r w 1 -ρ vr χ H2,ε χ H3,ε σ v σ r w 2 ∂ f n ∂r -r(w 1 + w 2 + 1) f n = 0 (54) with f n (0, w 1 , w 2 , r) = U n (0, w 1 , w 2 , r) = h(w 1 , w 2 ).
To solve Eq.( 54), we let

f n (τ, w 1 , w 2 , r) = h(w 1 , w 2 )Q(τ, w 1 , w 2 )e -L(τ )(w 1 + w 2 + 1)r. ( 55
)
We have the terminal condition f n (0, w 1 , w 2 , r) = h n (w 1 , w 2 ). Substituting Eq.(55) into Eq.(54), we obtain the following ODEs:

∂Q ∂τ + κθ(T -τ ) -ρ sr χ H1,ε χ H3,ε σ s σ r w 1 -ρ vr χ H2,ε χ H3,ε σ v σ r w 2 L(τ )(w 1 + w 2 + 1) Q - 1 2 χ H,ε σ 2 r L 2 (τ )(w 1 + w 2 + 1) 2 Q = 0 (56) ∂L ∂τ + κL -1 = 0 (57)
where Q(0, w 1 , w 2 ) = 1 and L(0) = 0. Then, by solving Eq.( 56), we have

Q(τ, w 1 , w 2 ) = exp χ H3,ε σ 2 r 2κ 2 τ + 2 κ * e -κ * τ -1 - 1 2κ * e -2κ * τ -1 (w 1 + w 2 + 1) 2 + exp ρ sr χ H1,ε χ H3,ε σ s σ r κ * w 1 (w 1 + w 2 + 1) + ρ vr χ H2,ε χ H3,ε σ v σ r κ * w 2 (w 1 + w 2 + 1) (τ -L(τ )) -exp τ 0 (w 1 + w 2 + 1)κ * θ * (T -ν)L(ν)dν and
L(τ ) = 1 -e -κ * τ κ * . ( 58
)
Then, from Eq.( 53), we obtain

U n (τ, w 1 , w 2 , r) = h n (w 1 , w 2 ) exp N (τ )(w 1 + w 2 + 1) 2 + exp 1 2 χ H1,ε σ 2 s w 1 (w 1 + 1) + ρ sv χ H1,ε χ H2,ε + 1 2 χ H2,ε σ 2 v w 2 (w 2 + 1) + exp ρ sr χ H1,ε χ H3,ε σ s σ r κ * w 1 (w 1 + w 2 + 1) + ρ vr χ H2,ε χ H3,ε σ v σ r κ * w 2 (w 1 + w 2 + 1) (τ -L(τ )) -exp {M (τ ) + L(τ )r(w 1 + w 2 + 1)} (59)
with N (τ ) and M (τ ) are given by

N (τ ) = χ H 3 ,ε σ 2 r 2(κ * ) 2 τ + 2 κ * e -κ * τ -1 -1 2κ * e -2κ * τ -1 M (τ ) = τ 0 κθ(T -ν)L(ν)dν (60)
We have

U n (τ, s, v, r) = 1 (2π) 2 a1+i∞ a1-i∞ a2+i∞ a2-i∞
h n (w 1 , w 2 )e J(τ,w1,w2,r) s -w1 v -w2 dw 1 dw 2 (61)

where J(τ, w 1 , w 2 , r) is given by

J(τ, w 1 , w 2 , r) = N (τ )(w 1 + w 2 + 1) 2 + 1 2 χ H2,ε σ 2 s w 1 (w 1 + 1) + ρ sv χ H1,ε χ H2,ε + 1 2 χ H2,ε σ 2 v w 2 (w 2 + 1) τ + ρ sr χ H2,ε χ H3,ε σ s σ r κ * w 1 (w 1 + w 2 + 1) + ρ vr χ H2,ε χ H3,ε σ v σ r κ * w 2 (w 1 + w 2 + 1) (τ -L(τ )) -(M (τ ) + L(τ )r) (w 1 + w 2 + 1). (62)
We define the following expression to use to compute U n (τ, s, v, r)

Λ(τ, s, v, r) = 1 (2π) 2 a1+i∞ a1-i∞ e C(τ,w2,r) a2+i∞ a2-i∞ h n (w 1 , w 2 )e {A(τ)ξ 2 (τ )} s -w1 dw 1 v -w2 dw 2 . ( 63
)
Now, the term J(τ, w 1 , w 2 , r) have the following quadratic form of w 1 and we obtain

Λ(τ, s, v, r) = 1 (2π) 2 a1+i∞ a1-i∞ e C(τ,w2,r) a2+i∞ a2-i∞ h n (w 1 , w 2 )e {A(τ)ξ 2 (τ )} s -w1 dw 1 v -w2 dw 2 (64) with ξ(τ ) = w 1 + A(τ )+B(τ )+B(τ )w2 2A(τ )
,where

         A(τ ) = N (τ ) + 1 2 χ H1,ε σ 2 s τ + ρsr √ χ H 1 ,ε χ H 3 ,ε σsσr κ * (τ -L(τ )) B(τ ) = 2N (τ ) + ρ sv χ H1,ε χ H2,ε + ρsr √ χ H 1 ,ε χ H 3 ,ε σsσr κ * + ρvr √ χ H 2 ,ε χ H 3 ,ε σvσr κ * (τ -L(τ )) B(τ ) = N (τ ) -M (τ ) -L(τ )r.
and

C(τ, w 2 , r) = - 1 4 A(τ ) + B(τ ) + B(τ )w 2 2 A(τ ) + N (τ )(w 2 + 1) 2 + 1 2 χ H2,ε σ 2 v (w 2 2 + w 2 )τ + ρ vr χ H2,ε χ H3,ε σ v σ r κ * (τ -L(τ )) -M (τ )(w 2 + 1) -L(τ )(w 2 + 1)r.
We observe that A(τ ) ⩾ 1 2 τ 0 (σ s -σ r L(ν)) 2 dν, we deduce that A(τ ) > 0, for all τ > 0, according Lemma 2.1, we have

1 2πi a2+i∞ a2-i∞ e {A(τ)ξ 2 (τ )} s -w1 dw 1 = 1 πA(τ ) s A(τ )+B(τ )+B(τ )w 2 2A(τ ) e -1 4A(τ ) (ln s) 2 .
Therefore, Eq.(64) becomes

Λ(τ, s, v, r) = e -1 4A(τ ) (ln s) 2 2 πA(τ ) s A(τ )+B(τ ) 2A(τ ) 1 2πi a1+i∞ a1-i∞ e C(τ,w2,r) s B(τ )w 2 2A(τ ) v -w2 dw 2 . (65)
In the Eq.( 65), the integrand have a quadratic form of w 2 as follows

1 2πi a1+i∞ a1-i∞ e C(τ,w2,r) s B(τ )w 2 2A(τ ) v -w2 dw 2 = exp      - B 3 (τ ) + B(τ ) 2A(τ ) ln s 2 4B 1 (τ ) -B 2 (τ ) 4A - A(τ ) -B 2 2 4A(τ )      × 1 2πi a1+i∞ a1-i∞ exp      B 1 (τ ) - B 2 (τ ) 4A(τ )   w 2 2 + B 3 (τ ) + B(τ ) 2A(τ ) ln s 2B 1 (τ ) -B 2 (τ ) 2A(τ )   2      v -w2 dw 2 , (66)
where

   B 1 (τ ) = N (τ ) + 1 2 χ H2,ε σ 2 v τ + ρvr √ χ H 2 ,ε χ H 3 ,ε σvσr κ * (τ -L(τ )) B 3 (τ ) = -B(τ ) 2A(τ ) A(τ ) + B 2 (τ ) + B 1 (τ ) + B 2 (τ ). By setting σ 2 s (t) = χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(t) + χ H3,ε σ 2 r L 2 (t) σ 2 v (t) = χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(t) + χ H3,ε σ 2 r L 2 (t). (67)
We have

     A(τ ) = 1 2 τ 0 σ 2 s (ν)dν B 1 (τ ) = 1 2 τ 0 σ 2 v (ν)dν B(τ ) = τ 0 ρ sv χ H1,ε χ H2,ε σ s σ v + σ r L(ν)(ρ sr χ H1,ε χ H3,ε σ s + ρ vr χ H2,ε χ H3,ε σ r ) + χ H3 σ 2 r L 2 (ν) dν. ( 68
)
While the correlation coefficients ρ sr and ρ vr are constants satisfying Eq.( 20), then B 1 (τ ) -B 2 (τ ) 4A(τ ) > 0 is satisfied for τ > 0. From Lemma 2.1 applied to Eq.(66), we have

1 2πi a1+i∞ a1-i∞ exp      B 1 (τ ) - B 2 (τ ) 4A(τ )   w 2 2 + B 3 (τ ) + B(τ ) 2A(τ ) ln s 2B 1 (τ ) -B 2 (τ ) 2A(τ )   2      v -w2 dw 2 = 1 2 π B 1 (τ ) -B 2 (τ ) 4A(τ ) exp      - 1 4 B 1 (τ ) -B 2 (τ ) 4A(τ ) (ln v) 2      v B 3 (τ )+ B(τ ) 2A(τ ) ln s 2B 1 (τ )- B 2 (τ ) 2A(τ ) (69)
Therefore, using Eq.(65), Eq.(66) and Eq.( 69), we obtain the following explicit of Λ(τ, s, v, r) which is given by

Λ(τ, s, v, r) = exp      - B 3 (τ ) + B(τ ) 2A(τ ) ln s 2 4 B 1 (τ ) -B 2 (τ ) 4A(τ ) - A(τ ) -B(τ ) 2 4A(τ )      × s A(τ )+B(τ )+B(τ )w 2 2A(τ ) v B 3 (τ )+ B(τ ) 2A(τ ) ln s 2 B 1 (τ )- B 2 (τ ) 4A(τ ) × e -1 4A(τ ) (ln s) 2 2 πA(τ ) e - 1 4 B 1 (τ )- B 2 (τ ) 4A(τ ) (ln v) 2 2 π B 1 (τ ) -B 2 (τ ) 4A(τ ) . ( 70
)
Now the expression of U n (τ, s, v, r) is obtain by using the Mellin convolution result given by Lemma 2.2, since h(w 1 , w 2 ) and e J(τ,s,v,r) are the double Mellin transforms of h n (s, v) and Λ(τ, s, v, r), respectively. By using Eq.( 45), Eq.( 61) becomes

U n (τ, s, v, r) = ∞ 0 ∞ 0 h n (u, w)Λ(τ, s u , v w , r)u -1 w -1 dudw = e k n γ * n K1 e ϕ1(u) (u -K 1 ) ϕ 2 (u, w) 1 u 1 w dudw -e k n γ * n K2 e ϕ1(u) (u -K 2 ) ϕ 2 (u, w) 1 u 1 w dudw + e k δ γ * 0 ∞ K1 e ϕ1(u) w (u -K 1 ) ϕ 2 (u, w) 1 u 1 w dudw -e k δ γ * 0 ∞ K2 e ϕ1(u) w (u -K 2 ) ϕ 2 (u, w) 1 u 1 w dudw. with k = - 1 4 A(τ ) -B(τ ) 2 A(τ ) , where      δ = 1-α γ ϕ 1 (u) = -1 4
B3(τ )+ B(τ ) 2A(τ ) ln s 2 B1(τ )- B 2 (τ ) 4A(τ )
, and

ϕ 2 (u, w) = s u A(τ )+B(τ ) 2A(τ ) v w B 3 (τ )+ B(τ ) 2A(τ ) ln ( s u ) B 1 (τ )- B 2 (τ ) 4A(τ ) × e -1 4A(τ ) (ln( s u )) 2 2 πA(τ ) e - 1 4 B 1 (τ )- B 2 (τ ) 4A(τ ) (ln( v w )) 2 2 π B 1 (τ ) -B 2 (τ ) 4A (τ )
.

By setting U (τ, s, v, r) = lim n-→∞ U n (τ, s, v, r), it follow that U (τ, s, v, r) = ∞ 0 ∞ 0 h n (u, w)Λ(τ, s u , v w , r)u -1 w -1 dudw = e k n γ * n K1 e ϕ1(u) (u -K 1 ) ϕ 2 (u, w) 1 u 1 w dudw -e k n γ * n K2 e ϕ1(u) (u -K 2 ) ϕ 2 (u, w) 1 u 1 w dudw + e k δ γ * 0 ∞ K1 e ϕ1(u) w (u -K 1 ) ϕ 2 (u, w) 1 u 1 w dudw -e k δ γ * 0 ∞ K2 e ϕ1(u) w (u -K 2 ) ϕ 2 (u, w) 1 u 1 w dudw = U 1 (τ, s, v, r) + V 1 (τ, s, v, r). (71)
We deduce that

U (τ, s, v, r) = U 1 (τ, s, v, r) + V 1 (τ, s, v, r), (72)
with

U 1 (τ, s, v, r) = U 1 (τ, s, v, r) -U 1 (τ, s, v, r) (73) V 1 (τ, s, v, r) = V 1 (τ, s, v, r) -V 1 (τ, s, v, r), (74)
where

U 1 (τ, s, v, r) = e k n γ * n K1 e ϕ1(u) (u -K 1 ) ϕ 2 (u, w) 1 u 1 w dudw (75) U 1 (τ, s, v, r) = e k n γ * n K2 e ϕ1(u) (u -K 2 ) ϕ 2 (u, w) 1 u 1 w dudw. ( 76
)
and

V 1 (τ, s, v, r) = e k δ γ * 0 ∞ K1 e ϕ1(u) w (u -K 1 ) ϕ 2 (u, w) 1 u 1 w dudw (77) V 1 (τ, s, v, r) = e k δ γ * 0 ∞ K2 e ϕ1(u) w (u -K 2 ) ϕ 2 (u, w) 1 u 1 w dudw. ( 78
)
The explicit calculation of U 1 (τ, s, v, r) and V 1 (τ, s, v, r) of Eq.( 72) leads to the following analytic closed form formula of the price of vulnerable bull spread call option. By letting U (t, s, v, r) = CallBullSpread(t, s, v, r), we have:

Theorem 5.3. The price of Bull Spread vulnerable call option denote by CallBullSpread(t, s, v, r) with the Mixed-Modified-Fractional-Hull-White-Vasicek model,defined by Eq.( 46), is explicitly given by

CallBullSpread(t, s, v, r) = sN ( a 1 , a 2 , ξ) -K 1 P (t, r, T ))N ( b 1 , b 2 , ξ) (79) -sN ( a 1 , a 2 , ξ) + K 2 B(t, r, T ))N ( b 1 , b 2 , ξ) + δv s p(t, T ) P (t, r, T ) N ( c 1 , c 2 , -ξ) -K 1 N ( d 1 , d 2 , -ξ) (80)
-δv s p(t, T ) P (t, r, T )

N ( c 1 , c 2 , -ξ) + K 2 N ( d 1 , d 2 , -ξ) ( 81
)
with ξ = B(τ ) 2 √ A(τ )B1(τ )
where

a 1 = ln s K1 -ln P (t, r, T ) + 1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (82) b 1 = ln v γ * -ln P (t, r, T ) -1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H2,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (83) a 2 = ln s K1 -ln P (t, r, T ) -1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H2,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (84) b 2 = ln v γ * -ln P (t, r, T ) -1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (85)
and

a 1 = ln s K2 -ln P (t, r, T ) + 1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (86) b 1 = ln v γ * -ln P (t, r, T ) -1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H2,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (87)
a 2 = ln s K2 -ln P (t, r, T ) -1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (88) b 2 = ln v γ * -ln P (t, r, T ) -1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (89) c 1 = ln s K1 -ln P (t, r, T ) + 1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (90) c 2 = - ln v γ * -ln P (t, r, T ) + 1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν - 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (91)
d 1 = ln s K1 -ln P (t, r, T ) -1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (92)
d 2 = - ln v γ * -ln P (t, r, T ) + 1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (93) c 1 = ln s K2 -ln P (t, r, T ) + 1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (94) c 2 = - ln v γ * -ln P (t, r, T ) + 1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν - 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (95)
d 1 = ln s K2 -ln P (t, r, T ) -1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (96)
d 2 = - ln v γ * -ln P (t, r, T ) + 1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (97)
Proof of Theorem 5.3. The proof of this result is found in Appendix A.3

In the case of a vulnerable bull spread put option with Hull-White-Vasicek volatility and interest rates model using Mixed Modified Fractional Brownian Motion, the option pricing formula is given by the following result:

Theorem 5.4. The price of Bull Spread vulnerable put option denote by P utBullSpread(t, s, v, r), with the Mixed-Modified-Fractional-Hull-White-Vasicek volatility and interest rates model, defined by Eq.( 46), is explicitly given by

P utBullSpread(t, s, v, r) = -sN (-a 1 , a 2 , ξ) + K 1 P (t, r, T ))N (-b 1 , b 2 , ξ) (98)
+ sN ( -a 1 , a 2 , ξ) -K 2 P (t, r, T ))N (-b 1 , b 2 , ξ) + δv -s p(t, T ) P (t, r, T ) N (-c 1 , c 2 , -ξ) + K 1 N (-d 1 , d 2 , -ξ) (99)
-δv -s p(t, T ) P (t, r, T ) N (-c 1 , c 2 , -ξ) -K 2 N (-d 1 , d 2 , -ξ) (100)
where

a 1 , a 2 , b 1 , b 2 , c 1 , c 2 , d 1 , d 2 , c 1 , c 2 , d 1 , d 2 ,
ξ and p(t, T ) are given by Theorem 5.3.

We deduce that the bull Call spread reduces the cost of the call option, but it comes with a trade-off. The gains in the stok's price are also capped, creating a limited range where the investor can make a profit. Traders will use the bull call spread if they believe an asset will moderately rise in value.

Conclusion

As an extension of Hull and White [START_REF] Hull | The pricing of options with stochastic volatilities[END_REF] and Vasicek [START_REF] Vasicek | An Equilibrium Characterisation of the Term Structure[END_REF] models, this paper proposes a new three-factor pricing model for vulnerable European options incorporating spreads and interest rates whose movements are governed by mixed modified fractional stochastic differential equations. The dynamics of the short-term interest rate is represented implicitly by a fractional stochastic process of the price of zero-coupon bond. Under our original assumption, we have obtain the bond price closed form formula. We have discretize the stochastic process with Milstein discretization scheme and we have showed the sample paths of asset price movement, the sample paths of volatility of asset price movement and the sample paths of interest rates of asset price movement. Furthermore, we exploit Mellin's transform techniques to derive a closed-form solution for vulnerable Bull Spread options model, which is simply calculated using the standard normal cumulative distribution function so that pricing and coverage of vulnerable European options can be calculated accurately and quickly.

The applications of our results are numerous, in particular to for a realistic credit risk measurement (see Johnson and Stuly [START_REF] Johnson | The pricing of options with default risk[END_REF]) and Klein [START_REF] Klein | Pricing Black-Scholes options with correlated credit risk[END_REF]) and Amman [START_REF] Ammann | Credit Risk Valuation: Methods, Models, and Applications[END_REF].

A Appendix

A.1 Notations

• µ the drift of the process for the stock and for the variance;

• σ s is the volatility of the stock price;

• σ v is the volatility of the variance;

• κ the mean reversion speed for the interest rate;

• θ the mean reversion level for the interest rate;

• σ r the volatility of the interest rate;

• µ the drift of the process for the stock and for the variance;

• S 0 > 0 is the spot asset price;

• V 0 > 0 is the initial (time zero)level of the variance;

• r 0 > 0 is the initial (time zero)level of the interest rates;

• Q is a Risk neutral probability measure and P is a probability measure;

• S t is the dynamic stock price,V t is its volatility and r t is its interest rate;

• P (t, r, T ) is the default price of a zero-coupon bond maturing at time T with the nominal value 1;

• U is the price of the bull spread vulnerable option at time t;

• N is the central distribution function of the normal distribution,N (0, 1);

• M H1,ε
.,s and M H2,ε .,v ,M H3,ε .,r are three Mixed Modified Fractional Brownian motions with Hurst parameters

H 1 , H 2 , H 3 ∈] 1 2 , 1[; • ε

is the adjustment coefficient between the times of cotations on the financial market;

• γ * is a critical value such that a credit loss occurs if the value of the option writer's assets;

• K 1 and K 2 are the strike of the bull spread option;

• p(t, T ) is the price at time 0 < t < T of a zero coupon with maturity T .

A.2 Proof of Theorem 5.2

Let U = U (s, v, r, t) be the price of a bull spread option and Π denotes the value of the portfolio at time t. We have

Π = U -n 1 s -n 2 v -n 3 P (101)
The variation of the portfolio is defined by

dΠ = dU -n 1 ds -n 2 v -n 3 dP (102)
By applying the fractional Ito's lemma to Eq.(102) as well as choosing

n 1 = ∂U ∂s ; n 2 = ∂U ∂v ; n 3 = ∂U ∂r ∂P ∂r , (103)
to eliminate the stochastic noise contained in the dS, dV and dr respectively. We obtain

dΠ = ∂U ∂t dt+ 1 2 ∂ 2 U ∂s 2 (ds) 2 + 1 2 ∂ 2 U ∂v 2 (dV ) 2 + ∂ 2 U ∂s∂v (dsdV )+ ∂ 2 U ∂s∂r (dsdr)+ ∂ 2 U ∂v∂r (dV dr)-n 3 ∂U ∂t dt + 1 2
∂ 2 U ∂r 2 (dr) 2 (104) with (dS t ) 2 , (dV t ) 2 (dr) 2 , dS t dV t , dS t dr dV t dr t are given respectively by

               (dS t ) 2 = χ H1,ε σ 2 s s 2 dt (dV t ) 2 = χ H2,ε σ 2 v v 2 dt (dr t ) 2 = χ H3,ε σ 2
r dt dS t t = χ H1,ε χ H2,ε σ s σ v svρ sv dt dS t dr t = χ H1,ε χ H3,ε σ s σ r sρ sr dt dV t dr t = χ H2,ε χ H3,ε σ v σ r vρ vr dt.

(105) Substituting Eq.(105) into Eq.( 104), we have

         dΠ = ∂U ∂t + 1 2 χ H1,ε σ 2 s s 2 ∂ 2 U ∂s 2 + 1 2 χ H2,ε σ 2 v ∂ 2 U ∂v 2 dt + χ H1,ε χ H2,ε σ s σ v svρ sv ∂ 2 U ∂s∂v + χ H1,ε χ H3,ε σ s σ r sρ sr ∂ 2 U ∂s∂r + χ H2,ε χ H3,ε σ v σ r vρ vr ∂ 2 U ∂v∂r dt -n 3 ∂U ∂t + 1 2 χ H3,ε σ 2 r ∂ 2 U ∂r 2 . ( 106
)
On the other hand, it is known that the return for any risk-free portfolio should be r. Hence, the variation value of the portofolio is

dΠ t = r U -s ∂U ∂s -v ∂U ∂v -n 3 P dt (107)
By substituting Eq.(106) and Eq.(38) into Eq.(107) gives

         ∂U ∂t + 1 2 χ H1,ε σ 2 s s 2 ∂ 2 U ∂s 2 + 1 2 χ H2,ε σ 2 v ∂ 2 U ∂v 2 dt + χ H1,ε χ H2,ε σ s σ v svρ sv ∂ 2 U ∂s∂v + χ H1,ε χ H3,ε σ s σ r sρ sr ∂ 2 U ∂s∂r + χ H2,ε χ H3,ε σ v σ r vρ vr ∂ 2 U ∂v∂r dt -n 3 ∂U ∂t + 1 2 χ H3,ε σ 2 r ∂ 2 U ∂r 2
. = rU -rs ∂U ∂s -rv ∂U ∂v -rn 3 P .

(108)

In the above statement Eq.(109), we use Eq.( 38) and the value of n 3 given by Eq.( 103) to obtain

∂U ∂t + 1 2 χ H1,ε σ 2 s s 2 ∂ 2 U ∂s 2 + 1 2 χ H2,ε σ 2 v v 2 t ∂ 2 U ∂v 2 + 1 2 χ H3,ε σ 2 r ∂ 2 U ∂r 2
+ ρ sv σ s σ v sv χ H1,ε χ H2,ε ∂ 2 U ∂s∂v + ρ sr σ s σ r s χ H1,ε χ H3,ε ∂ 2 U ∂s∂r

+ ρ vr σ v σ r v χ H2,ε χ H3,ε ∂ 2 U ∂v∂r + r s ∂U ∂s + v ∂U ∂v + κ * (θ * -r) ∂U ∂r -rU = 0 (109)
A.3 Proof of Theorem 5.3

A.3.1 Expression of U 1 (τ, s, v, r)
• Expression of U 1 (τ, s, v, r).We begin to take changes of variables

     x = ln( s u ) √ 2A(τ ) y = ln( s w ) √ 2B(τ ) . (110) and      η 1 (x) = e -xσs √ τ -K 1 Σ(τ ) = B 2 3 (τ ) 4 B1(τ )- B 2 (τ ) 4A(τ )
.

(111)

We have:

U 1 (τ, s, v, r) (112) = 1 2π B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) ∞ ln ( v γ * ) √ 2B 1 (τ ) ∞ ln ( s K 1 ) √ 2A(τ ) e k exp          - 1 4 B 3 (τ ) + B(τ ) √ 2A(τ ) x 2 B 1 (τ ) -B 2 (τ ) 4A(τ )          η 1 (x) × exp        A(τ ) + B(τ ) 2A(τ ) x + y B 1 (τ ) B 3 (τ ) + B(τ ) √ 2A(τ ) x √ 2 B 1 (τ ) -B 2 (τ ) 4A(τ )        e -1 2 x 2 e -1 2 B 1 (τ ) B 1 (τ )- B 2 (τ ) 4A(τ ) y 2 dxdy = s 2π B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) ∞ ln ( v γ * ) √ 2B 1 (τ ) ∞ ln ( s K 1 ) √ 2A(τ ) exp    k - 1 2 B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) (x 2 + y 2 )    × exp            -2A(τ ) + A(τ ) + B(τ ) 2A(τ ) - 1 2
B 3 (τ ) + B(τ ) √ 2A(τ ) B 1 (τ ) -B 2 (τ ) 4A(τ )     x -Σ(τ )        × exp        y B 1 (τ )B 3 (τ ) √ 2 B 1 (τ ) -B 2 (τ ) 4A(τ ) + B1(τ ) B1(τ )- B 2 (τ ) 2A(τ ) × B 1 (τ ) - B 2 (τ ) 4A(τ ) × B(τ ) 2 A(τ )B 1 (τ ) xy        dxdy × K 1 2π B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) ∞ ln ( v γ * ) √ 2B 1 (τ ) ∞ ln ( s K 1 ) √ 2A(τ ) exp    k - 1 2 B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) (x 2 + y 2 )    × exp            -2A(τ ) + A(τ ) + B(τ ) 2A(τ ) - 1 2
B 3 (τ ) + B(τ ) √ 2A(τ ) B 1 (τ ) -B 2 (τ ) 4A(τ )     x -Σ(τ )        × exp        y B 1 (τ )B 3 (τ ) √ 2 B 1 (τ ) -B 2 (τ ) 4A(τ ) + B1(τ ) B1(τ )- B 2 (τ ) 2A(τ ) × B 1 (τ ) - B 2 (τ ) 4A(τ ) × B(τ ) 2 A(τ )B 1 (τ ) xy        dxdy. (113)
According to Lemma 2.2, From Eq.(72), we let

U 1 (τ, s, v, r) = s 2π B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) ∞ ln ( v γ * ) √ 2B 1 (τ ) ∞ ln ( s K 1 ) √ 2A(τ ) exp    k - 1 2 B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) (x 2 + y 2 )    × exp            -2A(τ ) + A(τ ) + B(τ ) 2A(τ ) - 1 2
B 3 (τ ) + B(τ ) √ 2A(τ ) B 1 (τ ) -B 2 (τ ) 4A(τ )     x -Σ(τ )        × exp        y B 1 (τ )B 3 (τ ) √ 2 B 1 (τ ) -B 2 (τ ) 4A(τ ) + B1(τ ) B1(τ )- B 2 (τ ) 2A(τ ) × B 1 (τ ) - B 2 (τ ) 4A(τ ) × B(τ ) 2 A(τ )B 1 (τ ) xy        dxdy = s        1 2π √ Σ * ∞ ln ( v γ * ) √ 2B 1 (τ ) ∞ ln ( s K 1 ) √ 2A(τ ) exp - 1 2Σ * (x + a) 2 + (y + b) 2 + ξ 2Σ * (x + a)(y + b)        . where              Σ * = 1 -ξ 2 , a = A(τ )-B(τ ) √ 2A(τ ) , b = B(τ )-B1(τ )-B(τ ) √ 2B1(τ ) , -1 < ξ < 1.
Using the following transformation of variables,

                   x 1 = x + A(τ )-B(τ ) √ 2A(τ )
,

y 1 = y + B(τ )-B1(τ )-B(τ ) √ 2B1(τ ) , x 2 = x -A(τ )+B(τ ) √ 2A(τ ) , y 2 = y -B1(τ )+B(τ ) √ 2B1(τ ) . (114)
We obtain

U 1 (τ, s, v, r) = s 1 2π 1 -ξ 2 a1 -∞ a2 -∞ e 1 2(1-ξ 2 ) (x 2 1 -2ξx2y2+y 2 1 ) dx 1 dy 1 , (115)
where

â1 = ln s K1 -ln P (t, r, T ) + 1 2 τ 0 χ H2,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν â2 = ln s K1 -ln P (t, r, T ) -1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν
and according always to Lemma 2.2, we also let

U 2 (τ, s, v, r) = K 1 2π B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) ∞ ln ( v γ * ) √ 2B 1 (τ ) ∞ ln ( s K 1 ) √ 2A(τ ) exp    k - 1 2 B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) (x 2 + y 2 )    × exp            -2A(τ ) + A(τ ) + B(τ ) 2A(τ ) - 1 2
B 3 (τ ) + B(τ ) √ 2A(τ ) B 1 (τ ) -B 2 (τ ) 4A(τ )     x -Σ(τ )        × exp        y B 1 (τ )B 3 (τ ) √ 2 B 1 (τ ) -B 2 (τ ) 4A(τ ) + B1(τ ) B1(τ )- B 2 (τ ) 2A(τ ) × B 1 (τ ) - B 2 (τ ) 4A(τ ) × B(τ ) 2 A(τ )B 1 (τ ) xy        dxdy = K 1 P (t, r, T ) 2π √ Σ *        ∞ ln ( v γ * ) √ 2B 1 (τ ) ∞ ln ( s K 1 ) √ 2A(τ ) exp P (t, r, T ) - 1 2Σ * (x + c) 2 + (y + d) 2 + q(τ )        . where          q(τ ) = ξ 2Σ * (x + c)(y + d) , c = -A(τ )+B(τ ) √ 2A(τ ) , d = -B1(τ )+B(τ ) √ 2B1 (τ )
.

By also using the transformation of variables given by Eq.( 114), we have

U 2 (τ, s, v, r) = K 1 P (t, r, T ) 1 2π 1 -ξ 2 b1 -∞ b2 -∞ e 1 2(1-ξ 2 ) (x 2 2 -2ξx2y2+y 2
2 ) dx 2 dy 2

where b1 = ln v γ * -ln P (t, r, T ) -

1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν b2 = ln v γ * -ln P (t, r, T ) -1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν
Combining, Eq.( 115) and Eq.( 116), the value of U 1 (t, s, v, r) given by Eq.(75) becomes

U 1 (t, s, v, r) = U 1 (t, s, v, r) + U 2 (t, s, v, r) = sN ( a 1 , a 2 , ξ) -K 1 P (t, r, T ))N ( b 1 , b 2 , ξ).
Therefore,we have:

U 1 (t, s, v, r) = sN ( a 1 , a 2 , ξ) -K 1 P (t, r, T ))N ( b 1 , b 2 , ξ) (117)
• Expression of U 1 (t, s, v, r). We make the same transformation used to have the value of U 1 (t, s, v, r)

and then, we changes the value of η 1 (x) by η 2 (x) = e -xσs √ τ -K 2 . We deduce that, the expression of second term U 1 (t, s, v, r) given by Eq.( 76) is defined as follows:

U 1 (t, s, v, r) = sN ( a 1 , a 2 , ξ) -K 2 P (t, r, T ))N ( b 1 , b 2 , ξ) (118)
with

â1 = ln s K2 -ln P (t, r, T ) + 1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (119) b1 = ln v γ * -ln P (t, r, T ) -1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (120) â2 = ln s K2 -ln P (t, r, T ) -1 2 τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (121) b2 = ln v γ * -ln P (t, r, T ) -1 2 τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (122)
Therefore, combining Eq.( 117) and Eq.( 118), Eq.( 73) gives:

U 1 (t, s, v, r) = sN ( a 1 , a 2 , ξ)-K 2 P (t, r, T ))N ( b 1 , b 2 , ξ)-sN ( a 1 , a 2 , ξ)-K 2 P (t, r, T ))N ( b 1 , b 2 , ξ) (123) A.3.2 Expression of V 1 (τ, s, v, r) • Expression of V 1 (τ, s, v, r) The expression of V 1 (τ, s, v, r) of Eq.(77) becomes V 1 (τ, s, v, r) = δvs 2π B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) ln ( v γ * ) √ 2B 1 (τ ) ∞ ∞ ln ( s K 1 ) √ 2A(τ ) exp    k - 1 2 B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) (x 2 + y 2 )    × exp            -2A(τ ) + A(τ ) + B(τ ) 2A(τ ) - 1 2
B 3 (τ ) + B(τ ) √ 2A(τ ) B 1 (τ ) -B 2 (τ ) 4A(τ )     x -Σ(τ )        × exp        y B 1 (τ )B 3 (τ ) √ 2 B 1 (τ ) -B 2 (τ ) 4A(τ ) + B1(τ ) B1(τ )- B 2 (τ ) 2A(τ ) × B 1 (τ ) - B 2 (τ ) 4A(τ ) × B(τ ) 2 A(τ )B 1 (τ ) xy        dxdy × δvK 1 2π B 1 (τ ) B 1 (τ ) -B 2 (τ ) 2A(τ ) ln ( v γ * ) √ 2B 1 (τ ) ∞ ∞ ln ( s K 1 ) √ 2A(τ ) exp    k - 1 2 B 1 (τ )
Then, following the same procedure as for U 1 (τ, s, v, r), and using the following transformation defined by:

                  
x 3 = x + -B(τ )+A+B √ 2A y 3 = y + -B(τ )+B1(τ )+B √ 2B1(τ )

x 4 = x + -B(τ )-A+B √ 2A y 4 = y + -B(τ )+B1(τ ) √

2B1(τ ) , (125)
the expression of V 1 (τ, s, v, r) is given by: χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (130)

• Expression of V 1 (τ, s, v, r). The expression of V 1 (τ, s, v, r) of Eq.( 78) is obtained by following the same procedure as for V 1 (τ, s, v, r). Then, the expression of V 1 (τ, s, v, r) is given by:

By using Eq.( 123) and Eq.(136), Eq.(72) gives the expression of U (τ, s, v, r).

Definition 3 . 1 .

 31 (MMFHWV model) In the Mixed Modified Fractional Hull-White model Eq.(

  (a) H1 = 0.75 (b) H1 = 0.80 (c) H1 = 0.85 (d) H1 = 0.90

Figure 1 :

 1 Figure 1: The Stock price paths under MMFHWV model with S 0 = 115, V 0 = 0.01, r 0 = 0.9, κ = 0.2, θ = 0.2, σ = 0.8.

Figure 2 :

 2 Figure 2: The volatility paths under MMFHWV model with S 0 = 115, V 0 = 0.01, r 0 = 0.9, κ = 0.2, θ = 0.2, σ = 0.8.

Figure 3 :

 3 Figure 3: The Interest rate paths under MMFHWV model with S 0 = 115, V 0 = 0.01, r 0 = 0.9, κ = 0.2, θ = 0.2, σ = 0.8.

B 1 (B 1 y B 1

 111 (τ ) -B 2 (τ ) 4A(τ ) (τ )B 3 (τ ) √ 2 B 1 (τ ) -B 2 (τ ) 4A(τ ) + B1(τ ) B1(τ )-B 2 (τ ) 2A(τ ) × B 1 (τ ) -B 2 (τ ) 4A(τ ) × B(τ ) 2 A(τ )B 1 (τ )

V 1 (( x 2 3 -2ξ x3 y3+ y 2 3 ) d x 3 d y 3 - 2 τ 0 χχ 2 τ 0 χχ 2 τ 0 χ 2 (χ

 13232020202 τ, s, v, r) = δvse -B(τ )+B(τ ) δvsK 1 e -B(τ )+B(τ ) , T ) P (t, r, T )N ( c 1 , -c 2 , ξ) -K 1 N ( d 1 , -d 2 , -ξ) . -ln P (t, r, T ) + 1 H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν(127)c 2 = -ln v γ * -ln P (t, r, T ) + 1 H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν -H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν -ln P (t, r, T ) -1 H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν(129) d 2 = -ln v γ * -ln P (t, r, T ) + 1 2 τ 0

V 1 ( 2 τ 0 χ 2 τ 0 χχ 2 τ 0 χχ 2 τ 0 χ

 120202020 τ, s, v, r) = -δvse -B(τ )+B(τ ) T ) P (t, r, T )N ( c 1 , -c 2 , ξ) -K 2 N ( d 1 , -d 2 , -ξ) . -ln P (t, r, T ) + 1 H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν + 2ξ τ 0 σ 2 s (ν)dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν (132) c 2 = -ln v γ * -ln P (t, r, T ) + 1 H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν -H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν(133)andd 1 = ln s K2 -ln P (t, r, T ) -1 H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν(134)d 2 = -ln v γ * -ln P (t, r, T ) + 1 H2,ε σ 2 v + 2ρ vr χ H2,ε χ H3,ε σ v σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν τ 0 χ H1,ε σ 2 s + 2ρ sr χ H1,ε χ H3,ε σ s σ r L(ν) + χ H3,ε σ 2 r L 2 (ν) dν(135) By combining Eq.(126) and Eq.(131), Eq.(74) gives: V 1 (τ, s, v, r) = δv p(t, T ) P (t, r, T ) N ( c 1 , -c 2 , ξ) -K 1 N ( d 1 , -d 2 , -ξ) -δv p(t, T ) P (t, r, T ) N ( c 1 , -c 2 , ξ) -K 2 N ( d 1 , -d 2 , -ξ) .

  t,s dB t,v ] = ρ sv dt. B .,s and B .,v are the wiener processes; ρ sv ∈ [-1; 1] is the correlation between W t,s and W t,v .

	The parameters of the model are given in Appendix A.1
	Definition 2.3. Vasicek[33] The Vasicek model is a stochastic interest rate model and the corre-
	sponding Stochastic Differential Equation (SDE) as follows,

The hurst parameters H 1 , H 2 and H 3 can be verified H 3 > H 2 > H 1 > 1

	1 ,ε .,s	and M H 1 ,ε .,v	are two Mixed Modified Fractional Brownian motions with Hurst
	parameters H 1 , H 2 ∈] 1 2 , 1[.
	(c) The evolution of interest rate r t is given by the Mixed Modified Fractional Vasicek no-arbitrage
	model, namely :
			dr t = κ(θ -r t )dt + r t σ r dM H 3 ,ε t,r .	(17)
	where M H 3 ,ε .,r	is the Mixed Modified Fractional Brownian motion with Hurst parameter H 3 ∈
	] 1 2 , 1[.
	(d) 2 .
	(e) The correlation of the Mixed Modified Fractional Brownian motions M H 1 ,ε .,s	and M H 1 ,ε .,v