# European Option Pricing under Wishart Processes

Authors: Raphael Naryongo; Philip Ngare; Anthony Waititu
Journal: Journal of Mathematics
 (2021)

---

## Abstract

This study deals with a single risky asset pricing model whose volatility is described by Wishart affine processes. This multifactor model with two dependency matrices describing the correlation between the asset dynamic and Wishart processes makes it more flexible enough to fit the market data for short or long maturities. The aim of the study is to derive and solve the call option pricing problem under the double Wishart stochastic volatility model. The Fourier transform techniques combined with perturbation methods are employed in order to price the European call options. The numerical illustrations on pricing predictions show similar behavior of price movements under the double Wishart model with respect to the market price.


---

## Body

## 1. Introduction

The limitation of the standard Black and Scholes [1] model of not accounting for the fact that implied volatility of derivative products varies by strike and maturity; this makes it less flexible to reproduce certain market conditions observed on derivative prices. This caused the introduction of the Heston [2] model, which has been remarkably famous and heavily applied in financial markets due to its flexibility, financial interpretation of parameters, and analytical tractability property of the model since it belongs to the class of the affine process. This helps it to obtain the call price of a European option by inverting the Fourier transform and forming a closed-form solution of the log-price characteristic function.However, despite the Heston model popularity, the existing literature documented the limitation of the model. Da Fonseca et al. [3], Christoffersen et al. [4], Ahdida and Alfonsi [5], Kang et al. [6], La Bua and Marazzina [7], and Gourieroux [8] pointed out that the major deficiency of the model does not generate the realistic term structure of the volatility smiles. The Heston model provides too flat implied volatility surface to attain reality. Yet, the implied volatility curve generally has steep slope, and it is convex in short maturity and tends to be linear for long maturity.The literature recommends that the problem of not capturing such stylized facts can be solved by generalizing the Heston model using two approaches: by adding jump in the stock dynamic or volatility and investigating the multifactor nature of the implied volatility as emphasized in the studies by Benabid et al. [9], Da Fonseca et al. [3], La Bua and Marazzina [7], and Kang and Kang [10]. It is well understood that the multifactor approach is capable of handling the pricing problem of derivative products and volatility skew. This usher in the application of the stochastic matrix defined process, that is, the Wishart multidimensional stochastic volatility model. The Wishart process was introduced in Bru [11] and is defined as a positive semidefinite matrix-valued generalization of the square root process. The square root process is a Bessel process which is a generalization of the multiple chi-squared distribution [12, 13]. The Wishart stochastic volatility model found its application in finance by Gourieroux and Sufana [14], and it has been applied in financial markets because of its nature of matrix specification which makes it flexible.The aim of this study is to price European option under two Wishart variance processes, through generalization of the Heston model into the multifactor nature of volatility modeling with two dependence matrices for single asset, that is, the asset dynamics depends on two Wishart volatility processes (double Wishart volatility processes). Following Kang et al. [6] and Benabid et al. [9], the matrix specification of the Wishart model makes it flexible enough to describe the market prices. The numerical illustrations show that price predictions under the double Wishart model have a similar behavior like the market prices.The European call option pricing problem is handled through transform methods of Duffie et al. [15] and fast Fourier transform by Carr and Madan [16] combined with perturbation techniques of Benabid et al. [9] and Fouque et al. [17], in order to obtain the European call option price formula.The study is structured as follows: Section2 presents the definition of the Wishart process, and we introduce the Wishart stochastic volatility model and double Wishart model, correlation structures, and the infinitesimal generator of the Bi-Wishart model. Section 3 provides the solution for European call option pricing problem through Fourier transform methods combined with perturbation techniques. Section 4 presents the numerical illustrations to compare the price predictions under the double Wishart model with the market prices. Section 5 provides conclusion and further research work, together with technical proofs in the Appendix.


## 2. Wishart Process

Definition 1.
LetWt,t≥0 under the probability measure ℚ be a n×n matrix-valued Brownian motion. The Wishart matrix process Σ satisfies the following equation:(1)dΣt=βQQT+MΣt+MΣtTdt+ΣtdWtQ+QTdWtTΣt,where Q∈GLnℝ is the invertible matrix, M∈Mn is the nonpositive matrix, and Σ0∈S˜n+ is the nonnegative symmetric matrix, while β is a real parameter. The condition β>n−1 is considered to ensure existence and uniqueness of the Σt∈S˜n+ solution for equation (1), and Eigen values of the solution are nonnegative for all t≥0 a.s Σt∈S˜n+. Like in the study by Benabid et al. [9], the probability measure ℚ corresponds to a risk-neutral measure (Appendix: change of the probability measure).


### 2.1. Wishart Volatility Model in the Stock Market

In Da Fonseca et al. [3], under the risk-neutral probability measure, the risky asset price dynamic and its quadratic variation are follows:(2)dStSt=rdt+TrΣtdZt,S0=s,dΣt=βQQT+MΣt+MΣtTdt+ΣtdWtQ+QTdWtTΣt,Σ0=Σ,where r denotes risk-less interest rate, Tr is the trace, Z∈Mn is the Brownian matrix, and Σt is in a set of symmetric n×n positive-definite matrices. We observe that asset volatility is a trace of Σt matrix, which is a multidimensional process with M,Q∈Mn, and Wt∈Mn is the Brownian matrix.In the study by Bru [11], the Wishart process provides a matrix analogue of the square root mean-reverting process and M is considered negative to ensure positivity and mean-reverting property of the volatility with parameter β>n−1 for uniqueness and existence of the solution.


### 2.2. Correlation Structure

The Brownian matricesWt and Zt are correlated, in such a way that it gives a constant correlated matrix R∈Mn, like in the study by Da Fonseca et al. [3], which describes the correlation structure, in such a way that Zt can be presented as follows:(3)Zt≔WtRT+BtI−RRT,where I denotes the identity matrix, T is the transpose, and Bt is an independent matrix Brownian motion from Wt. The correlation structure is a Brownian motion (Appendix).


### 2.3. Bivariate Wishart Stochastic Processes in Stock Market

In this section, we present a proposed novel model, the multifactor model with two Wishart variance processes or double Wishart stochastic volatility model with two dependence matrices. The model takes two volatility components which is the trace of Wishart whose diagonal components will be the controlling factors for the dynamics of volatilities.Under arbitrage-free financial market and probability measure, we consider risky asset dynamic as follows:(4)dStSt=rdt+TrΣtdZt+Σ¯tdZ¯t,S0=s,dΣt=βQQT+MΣt+ΣtMTdt+ΣtdWtQ+QTdWtTΣt,Σ0=Σ,dΣ¯t=β¯Q¯Q¯T+M¯Σ¯t+Σ¯tM¯Tdt+Σ¯tdW¯tQ¯+Q¯TdW¯tTΣ¯t,Σ¯0=Σ¯,where βandβ¯ are the real parameters, such that β,β¯>n−1, Q,Q¯,M,M¯∈Mn, Q is the invertible matrix, and Wt,W¯t∈Mn are the matrices and Brownian motions, also Zt,Z¯t∈Mn.Lemma 1.
The correlations between the Brownian matrices of the stock price dynamic and the Brownian matrices of the Wishart processes in equation (4) are given by(5)ρt=TrRTQΣtTrΣtTrQTQΣt,ρt¯=TrR¯TQ¯Σ¯tTrΣ¯tTrQ¯TQ¯Σ¯t.Proof.
We derive the correlations as follows:(6)dStSt=rdt+TrΣtdZt+Σ¯tdZ¯t=rdt+TrΣtTrΣtdZtTrΣt+TrΣ¯tTrΣ¯tdZ¯tTrΣ¯t=rdt+TrΣtdXt+TrΣ¯tdX¯t,where Xt and X¯t are the standard Brownian motions (Appendix), and also taking the trace of the Wishart volatility dynamics (4), we have(7)dTrΣt=βTrQTQ+2TrMΣtdt+2TrQdWtΣt,dTrΣ¯t=β¯TrQ¯TQ¯+2TrM¯Σ¯tdt+2TrQ¯dW¯tΣ¯t.
These processes can still be written in the following form:(8)dTrΣt=βTrQTQ+2TrMΣtdt+2TrQTQΣtTrQdWtΣtTrQTQΣt,dTrΣ¯t=β¯TrQ¯TQ¯+2TrM¯Σ¯tdt+2TrQ¯TQ¯Σ¯tTrQ¯dW¯tΣ¯tTrQ¯TQ¯Σ¯t.
Then,(9)dTrΣt=βTrQTQ+2TrMΣtdt+2TrQTQΣtdξt,dTrΣ¯t=β¯TrQ¯TQ¯+2TrM¯Σ¯tdt+2TrQ¯TQ¯Σ¯tdηt,where ξt and ηt are the Brownian motions (Appendix).
We now determine the covariation of the stock price and Wishart processes:(10)CovtdXt,dξt=CovtTrΣtdZtTrΣt,TrQdWtΣtTrQTQΣt=EtTrΣtdWtRTTrΣtTrQdWtΣtTrQTQΣt=∑ijCovteiTΣtdWtRTei,ejTQdWtΣtejTrΣtTrQTQΣt=∑ijEteiTΣtdWtRTeiejTQTdWtTΣtejTrΣtTrQTQΣt=∑ijeiTΣtTrRTeiejTQTΣtejdtTrΣtTrQTQΣt=∑ijTrQReiejTeiTΣtejdtTrΣtTrQTQΣt=∑ijejTQReieiTΣtejdtTrΣtTrQTQΣt=∑ijejTQRΣtejdtTrΣtTrQTQΣtCovtdXt,dξt=TrRTQΣtTrΣtTrQTQΣt.
Similarly, for the second stochastic differential equation of the Wishart processes, the determination of covariation follows the same procedures, that is,(11)CovtdX¯t,dηt=TrR¯TQ¯Σ¯tTrΣ¯tTrQ¯TQ¯Σ¯t.


### 2.4. The Correlation Structure of the Model

The correlated Brownian motionsWt, Zt and W¯t, Z¯t, respectively, result in constant correlations R,R¯∈Mn, which describe the two respective correlation structures, such that Zt and Z¯t can be presented as follows:(12)Zt≔WtRT+BtI−RRT,Z¯t≔W¯tR¯T+B¯tI−R¯R¯T,where I is the identity matrix, T is the transposition, and Bt and B¯t are the Brownian matrices independent of Wt and W¯t, respectively.


### 2.5. The Dynamic of Log-Price for the Double Wishart Model

The matricesR and R¯ describe the correlations between the Brownian of the asset price and those of the Wishart processes.Lemma 2.
The log-price dynamicYt=logSt under the double Wishart model is given as(13)dYt=r−12TrΣt+Σt¯dt+TrΣtdZt+Σt¯dZ¯t.Proof.
LetYt=logSt, and the asset dynamic is given as(14)dStSt=rdt+TrΣtdZt+Σt¯dZ¯t.
By applying Ito’s formula onYt, referring the studies by Björk [18] and Shreve [19], we get(15)dYt=dlogSt=dStSt−12dSt2St2.
Substituting the asset process (14) in the derivative expression (15) of Yt,(16)dYt=rdt+TrΣtdZt+Σt¯dZ¯t−12TrΣt+Σ¯tdt,dYt=r−12TrΣt+Σt¯dt+TrΣtdZt+Σt¯dZ¯t,which can be written as(17)dYt=r−12TrΣt+Σ¯tdt+TrΣtdWtRT+dBtI−RRT+Σt¯dW¯tR¯T+dB¯tI−R¯R¯T,Y0=y.


## 2.1. Wishart Volatility Model in the Stock Market

In Da Fonseca et al. [3], under the risk-neutral probability measure, the risky asset price dynamic and its quadratic variation are follows:(2)dStSt=rdt+TrΣtdZt,S0=s,dΣt=βQQT+MΣt+MΣtTdt+ΣtdWtQ+QTdWtTΣt,Σ0=Σ,where r denotes risk-less interest rate, Tr is the trace, Z∈Mn is the Brownian matrix, and Σt is in a set of symmetric n×n positive-definite matrices. We observe that asset volatility is a trace of Σt matrix, which is a multidimensional process with M,Q∈Mn, and Wt∈Mn is the Brownian matrix.In the study by Bru [11], the Wishart process provides a matrix analogue of the square root mean-reverting process and M is considered negative to ensure positivity and mean-reverting property of the volatility with parameter β>n−1 for uniqueness and existence of the solution.


## 2.2. Correlation Structure

The Brownian matricesWt and Zt are correlated, in such a way that it gives a constant correlated matrix R∈Mn, like in the study by Da Fonseca et al. [3], which describes the correlation structure, in such a way that Zt can be presented as follows:(3)Zt≔WtRT+BtI−RRT,where I denotes the identity matrix, T is the transpose, and Bt is an independent matrix Brownian motion from Wt. The correlation structure is a Brownian motion (Appendix).


## 2.3. Bivariate Wishart Stochastic Processes in Stock Market

In this section, we present a proposed novel model, the multifactor model with two Wishart variance processes or double Wishart stochastic volatility model with two dependence matrices. The model takes two volatility components which is the trace of Wishart whose diagonal components will be the controlling factors for the dynamics of volatilities.Under arbitrage-free financial market and probability measure, we consider risky asset dynamic as follows:(4)dStSt=rdt+TrΣtdZt+Σ¯tdZ¯t,S0=s,dΣt=βQQT+MΣt+ΣtMTdt+ΣtdWtQ+QTdWtTΣt,Σ0=Σ,dΣ¯t=β¯Q¯Q¯T+M¯Σ¯t+Σ¯tM¯Tdt+Σ¯tdW¯tQ¯+Q¯TdW¯tTΣ¯t,Σ¯0=Σ¯,where βandβ¯ are the real parameters, such that β,β¯>n−1, Q,Q¯,M,M¯∈Mn, Q is the invertible matrix, and Wt,W¯t∈Mn are the matrices and Brownian motions, also Zt,Z¯t∈Mn.Lemma 1.
The correlations between the Brownian matrices of the stock price dynamic and the Brownian matrices of the Wishart processes in equation (4) are given by(5)ρt=TrRTQΣtTrΣtTrQTQΣt,ρt¯=TrR¯TQ¯Σ¯tTrΣ¯tTrQ¯TQ¯Σ¯t.Proof.
We derive the correlations as follows:(6)dStSt=rdt+TrΣtdZt+Σ¯tdZ¯t=rdt+TrΣtTrΣtdZtTrΣt+TrΣ¯tTrΣ¯tdZ¯tTrΣ¯t=rdt+TrΣtdXt+TrΣ¯tdX¯t,where Xt and X¯t are the standard Brownian motions (Appendix), and also taking the trace of the Wishart volatility dynamics (4), we have(7)dTrΣt=βTrQTQ+2TrMΣtdt+2TrQdWtΣt,dTrΣ¯t=β¯TrQ¯TQ¯+2TrM¯Σ¯tdt+2TrQ¯dW¯tΣ¯t.
These processes can still be written in the following form:(8)dTrΣt=βTrQTQ+2TrMΣtdt+2TrQTQΣtTrQdWtΣtTrQTQΣt,dTrΣ¯t=β¯TrQ¯TQ¯+2TrM¯Σ¯tdt+2TrQ¯TQ¯Σ¯tTrQ¯dW¯tΣ¯tTrQ¯TQ¯Σ¯t.
Then,(9)dTrΣt=βTrQTQ+2TrMΣtdt+2TrQTQΣtdξt,dTrΣ¯t=β¯TrQ¯TQ¯+2TrM¯Σ¯tdt+2TrQ¯TQ¯Σ¯tdηt,where ξt and ηt are the Brownian motions (Appendix).
We now determine the covariation of the stock price and Wishart processes:(10)CovtdXt,dξt=CovtTrΣtdZtTrΣt,TrQdWtΣtTrQTQΣt=EtTrΣtdWtRTTrΣtTrQdWtΣtTrQTQΣt=∑ijCovteiTΣtdWtRTei,ejTQdWtΣtejTrΣtTrQTQΣt=∑ijEteiTΣtdWtRTeiejTQTdWtTΣtejTrΣtTrQTQΣt=∑ijeiTΣtTrRTeiejTQTΣtejdtTrΣtTrQTQΣt=∑ijTrQReiejTeiTΣtejdtTrΣtTrQTQΣt=∑ijejTQReieiTΣtejdtTrΣtTrQTQΣt=∑ijejTQRΣtejdtTrΣtTrQTQΣtCovtdXt,dξt=TrRTQΣtTrΣtTrQTQΣt.
Similarly, for the second stochastic differential equation of the Wishart processes, the determination of covariation follows the same procedures, that is,(11)CovtdX¯t,dηt=TrR¯TQ¯Σ¯tTrΣ¯tTrQ¯TQ¯Σ¯t.


## 2.4. The Correlation Structure of the Model

The correlated Brownian motionsWt, Zt and W¯t, Z¯t, respectively, result in constant correlations R,R¯∈Mn, which describe the two respective correlation structures, such that Zt and Z¯t can be presented as follows:(12)Zt≔WtRT+BtI−RRT,Z¯t≔W¯tR¯T+B¯tI−R¯R¯T,where I is the identity matrix, T is the transposition, and Bt and B¯t are the Brownian matrices independent of Wt and W¯t, respectively.


## 2.5. The Dynamic of Log-Price for the Double Wishart Model

The matricesR and R¯ describe the correlations between the Brownian of the asset price and those of the Wishart processes.Lemma 2.
The log-price dynamicYt=logSt under the double Wishart model is given as(13)dYt=r−12TrΣt+Σt¯dt+TrΣtdZt+Σt¯dZ¯t.Proof.
LetYt=logSt, and the asset dynamic is given as(14)dStSt=rdt+TrΣtdZt+Σt¯dZ¯t.
By applying Ito’s formula onYt, referring the studies by Björk [18] and Shreve [19], we get(15)dYt=dlogSt=dStSt−12dSt2St2.
Substituting the asset process (14) in the derivative expression (15) of Yt,(16)dYt=rdt+TrΣtdZt+Σt¯dZ¯t−12TrΣt+Σ¯tdt,dYt=r−12TrΣt+Σt¯dt+TrΣtdZt+Σt¯dZ¯t,which can be written as(17)dYt=r−12TrΣt+Σ¯tdt+TrΣtdWtRT+dBtI−RRT+Σt¯dW¯tR¯T+dB¯tI−R¯R¯T,Y0=y.


## 3. Option Pricing Problem

This part deals with the European call option pricing problem, given its payoff as(18)ST−K+.To handle this pricing problem, the infinitesimal generator of the Wishart processes should be obtained, in order to employ the conditional characteristic function on the log asset return. The Riccati ordinary differential equations are linearized to give the closed-form solution to the pricing problem through fast Fourier transforms of Duffie et al. [15]. This is possible because the Wishart processes preserve the property of analytical tractability, since it belongs to the class of affine [20].


### 3.1. Infinitesimal Generator

The log-price process and Wishart processes, together with corresponding pair of two correlated Brownian motions,ZsY,ZsΣ and Z¯sY,Z¯sΣ¯, can be written in the form as in the study by Benabid et al. [9], for easy handling of the complexity of the dynamics:(19)dYt=r−12TrΣt+Σt¯dt+TrΣtdZtY+Σt¯dZ¯tY,dTrΣt=βTrQTQ+2TrMΣtdt+2TrQTQΣtρtdZtY+I−ρt2dZtΣ,dTrΣ¯t=β¯TrQ¯TQ¯+2TrM¯Σ¯tdt+2TrQ¯TQ¯Σ¯tρ¯tdZ¯tY+I−ρ¯t2dZ¯tΣ¯,where(20)dZY,ZΣt=ρt=TrRTQΣtTrΣtTrQTQΣt,dZ¯Y,Z¯Σ¯t=ρ¯t=TrR¯TQ¯Σ¯tTrΣ¯tTrQ¯TQ¯Σ¯t.Proposition 1.
The infinitesimal generator of the double Wishart stochastic volatility model forYt,Σt,Σ¯t is given by(21)ℒY,Σ,Σ¯=r−TrΣ+Σ¯2∂∂y+TrΣ+Σ¯2∂2∂y2+βTrQTQ+2TrMΣ∂∂Σ+2TrΣ∂∂ΣQTQ∂∂Σ+β¯TrQ¯TQ¯+2TrM¯Σ¯∂∂Σ¯+2TrΣ¯∂∂Σ¯Q¯TQ¯∂∂Σ¯+2TrΣRQ∂∂Σ∂∂y+2TrΣ¯R¯Q¯∂∂Σ¯∂∂y.Proof.
The infinitesimal generator has a nontrivial term which arises from covariationdΣθij,Y corresponding to the coefficients of the term(22)∂2∂xθ;ij∂y=∂∂xθ;ij∂∂y,i,j=1,…,n,θ=1,2.
LetVθ;t≔Σθ;t be the square root matrix, such that(23)Σθ;tij=∑t=1nVθ;tilVθ;tlj=∑t=1nVθ;tilVθ;tjl.
Due to the fact thatVθ;t is symmetric, we find the covariation terms matching with ∂2/∂xθ;ij∂y coefficients:(24)dΣθij,Y=Et∑i,k=1nVθ;tildWlkθQkjθ+∑l,k=1nVθ;tjldWlkθQkiθ∑l,k,h=1nVθ;tlkdWkhθRlhθ=∑l,k,h=1nVθ;tilQkjθ+Vθ;tjlQkiθVθ;thlRhkθdt=∑k,h=1n∑l=1nVθ;tilVθ;thlQkjθ+∑l=1nVθ;tjlVθ;thlQkiθRhkθdt=∑k,h=1nΣθ;tihQkjθ+Σθ;tjhQkiθRhkθdt.
This gives the corresponding coefficients of the term since(25)2TrΣθRθQθDθ∂∂y=2∑i,j,k,h=1nDθijΣθjhRhkθQkiθ∂∂y.
The notation whenθ=1,Σ1=Σ,R1=R,Q1=Q,D1=∂/∂xij for θ=2,Σ2=Σ¯,R2=R¯,Q2=Q¯,D2=∂/∂x¯ij.


### 3.2. The Laplace Transform of Asset Returns

Following Duffie et al. [15] and Da Fonseca et al. [3], in order to solve the European options pricing problem, we consider the Laplace transform of the process (17). The Laplace transform of the Wishart process is exponentially affine as in the study by Bru [11]; hence, the conditional moment generating function of the log-asset returns is the exponential of an affine combination of Y and Wishart components, so we provide for the deterministic functions λ1t,λ2t∈Mn, and δt,εt∈ℝ, as parameters of the Laplace transform:(26)ψγ,tτ=EeγYt+τ=expTrλ1τΣt+Trλ2τΣ¯t+δτYt+ετ,where γ∈ℝ and Et denote the conditional expected value with respect to the probability measure.Using Feynman–Kac argument helps us to get matrix Riccati equations:(27)∂ψγ,t∂τ=ℒY,Σ,Σ¯ψγ,t,ψγ,t0=expγYt.Proposition 2.
The Laplace transform of the asset returns(28)ψγ,tτ=expTrλ1τΣt+λ2τΣ¯t+δτYt+ετ,where λ1,λ2,andε are the solutions of the differential equations through linearization methods:(29)λ1′τ=Mλ1τ+MT+2γRQλ1τ+2λ1τQTQλ1τ+γγ−12In,λ2′τ=M¯λ2τ+M¯T+2γR¯Q¯λ2τ+2λ2τQ¯TQ¯λ2τ+γγ−12In,ε′τ=rγ+βTrQTQλ1τ+β¯TrQ¯TQ¯λ2τ,with boundary conditions λ10=0,λ20=0,ε0=0,andδτ=γ=C0.The solutionsλ1,λ2,andε are obtained as follows:(30)λ1τ=H1τ−1H2τ,λ2τ=I1−1τI2τ,ετ=−β2TrlogH1τ+MT+2γRQτ−β¯2TrlogI1τ+M¯T+2γR¯Q¯τ+rγτ,with(31)H2τH1τ=H20H10expτM−2QTQγγ−12In−MT+2γRQ,I2τI1τ=I20I10expτM¯−2Q¯TQ¯γγ−12In−M¯T+2γR¯Q¯.Proof.
From equation (27) and the Laplace transform (26), we solve the problem by considering Proposition 1 for the double Wishart stochastic volatility model as follows:(32)∂ψγ,t∂τ=r−TrΣ+Σ¯2∂ψγ,t∂y+TrΣ+Σ¯2∂2ψγ,t∂y2+βTrQTQ+2TrMΣ∂ψγ,t∂Σ+2TrΣ∂∂ΣQTQ∂∂Σψγ,t+β¯TrQ¯TQ¯+2TrM¯Σ¯∂ψγ,t∂Σ¯+2TrΣ¯∂∂Σ¯Q¯TQ¯∂∂Σ¯ψγ,t+2TrΣRQ∂∂Σ∂ψγ,t∂y+2TrΣ¯R¯Q¯∂∂Σ¯∂ψγ,t∂y.
From equation (27), we have(33)λ10=0,λ20=0,δ0=γ,ε0=0.
Then,(34)∂ψγ,tτ∂τ=Trdλ1τdτΣt+dλ2τ∂τΣ¯t+ddτδtYt+ddτετψγ,tτ=r−Σ+Σ¯2δτψγ,t+TrΣ+Σ¯2δ2τψγ,t+βTrQTQ+2TrMΣλ1τ+β¯TrQ¯TQ¯+2TrM¯Σ¯λ2τψγ,t+2TrΣλ1τQTQλ1τψγ,tτ+2TrΣ¯λ1τQ¯TQ¯λ1τψγ,t+2TrΣRQδτλ1τψγ,t+2TrΣ¯R¯Q¯δτλ1τψγ,t,−Trλ´1τΣt+λ´2τΣ¯t+δ´τYt+ε´τ+r−TrΣ+Σ¯2δτ+TrΣ+Σ¯2δ2τ+TrβQTQ+MΣ+ΣMTλ1τ+β¯Q¯TQ¯+M¯Σ¯+Σ¯M¯Tλ2τ+2Σλ1τQTQλ1τ+2Σ¯λ2τQ¯TQ¯λ2τ+2ΣRQδτλ1τ+2Σ¯R¯Q¯δτλ2τ=0,−Trddτλ1τΣ−12TrΣδτ+12TrΣδ2τ+TrMΣ+ΣMTλ1τ+2Σλ1τQTQλ1τ+2ΣRQδτλ1τ=0,−Trddτλ2τΣ¯−12TrΣ¯δτ+12TrΣ¯δ2τ+TrM¯Σ¯+Σ¯M¯Tλ2τ+2Σ¯λ2τQ¯TQ¯λ2τ+2Σ¯R¯Q¯δτλ2τ=0,−δ´τYt−ε´τ+μδτ+TrβQTQλ1τ+Trβ¯Q¯TQ¯λ2τ=0.
Now, we need to identify the coefficients for above equations and obtain the matrix Riccati ordinary differential equations:(35)ddτλ1τ=Mλ1τ+MT+2γRQλ1τ+2λ1τQTQλ1τ+γγ−12In,λ10=0,(36)ddτλ2τ=M¯λ2τ+M¯T+2γR¯Q¯λ2τ+2λ2τQ¯TQ¯λ2τ+γγ−12In,λ20=0.
For the constantε, the matrix Riccati ordinary differential is obtained as follows:(37)−δ′τYt−ε′τ+rδτ+TrβQTQλ1τ+Trβ¯Q¯TQ¯λ2τ=0,δτ=C0=γ,sinceδ′τ=0so⇒δτ=γ,dετdτ=rγ+βTrQTQλ1τ+β¯TrQ¯TQ¯λ2τ,ε0=0.
Then,ετ is obtained by integrating directly(38)ετ=∫0τrγds+∫0τTrβQTQλ1s+β¯Q¯TQ¯λ2sds,with λ1,λ2∈Mnℝ, and δτ,ετ∈ℝ.
The matrix Riccati equations above are linearized to provide the closed-form solution, following the techniques of Benabid et al. [9] and Da Fonseca et al. [3]. So from equations (35) and (36), let(39)λ1τ=H1τ−1H2τ,where H1τ∈GLnℝ, H2τ∈Mnℝ, and therefore,(40)ddτH1τλ1τ=dH1τdτλ1τ+H1τdλ1τdτ,H1τdλ1τdτ=ddτH1τλ1τ−dH1τdτλ1τ,H1τλ1τM+H1τMT+2γRQλ1τ+2H1τλ1τQTQλ1τ+H1τγγ−12=H1τdλ1τdτ,H2τM+H1τMT+2γRQλ1τ+2H2τQTQλ2τ+H1τγγ−12=H1τdλ1τdτ.
SinceH2τ=H1τλ1τ, now we have(41)dH2τdτ=dH1τdτλ1τ+H1τdλ1τdτ.
This implies that(42)H1τdλ1τdτ=dH2τdτ−dH1τdτλ1τ.
Then, we get expressions(43)dH2τdτ−dH1τdτλ1τ=H2τM+H1τMT+2γRQλ1τ+2H2τQTQλ1τ+H1τγγ−12=H1τγγ−12+H2τM+H1τMT+2γRQ+2H2τQTQλ1τ,(44)dH2τdτ=H1τγγ−12+H2τM,(45)dH1τdτ=−H1τMT+2γRQ−2H2τQTQ,(46)ddτH2τH1τ=H2τH1τM−2QTQγγ−12In−MT+2γRQ.
The solution of differential equation (46) above is obtained through exponentiation as follows:(47)H2τH1τ=H20H10eQτ,with conditions H10=KIn and H1−10=K−1In and for K=1,H1−10=In, such that(48)H2τH1τ=λ10IneτQ,where(49)eτQ=λ111τλ112τλ121τλ122τ.
Then,(50)H2τH1τ=λ10Inλ111τλ112τλ121τλ122τ=λ10λ111τ+λ121τλ10λ112τ+λ122τ,and since λ10=0, then(51)H2τH1τ=λ121τλ122τ.
In conclusion, we get(52)λ1τ=λ122τ−1λ121τ.
This represents the closed-form solution of the Riccati equation (35).
Now, we look at the solution of the second Riccati equation (36). Let(53)λ2τ=I1−1τI2τ,then,(54)dI2τdτ=I1τγγ−12+I2τM¯,(55)dI1τdτ=−I1τM¯T+2γR¯Q¯−2I2τQ¯TQ¯,(56)ddτI2τI1τ=I2τI1τM¯−2Q¯TQ¯γγ−12In−M¯T+2γR¯Q¯,such that(57)I2τI1τ=I20I10eQ¯τ.
From equation (40), we go through the same procedures as above, with conditions(58)I20=I10λ20=λ20,I10=K¯I=K¯=I,therefore,(59)I2τI1τ=λ20Inλ211τλ212τλ221τλ222τ,I2τI1τ=λ20λ211τ+λ221τλ20λ212τ+λ222τ,since λ20=0, then(60)I2τI1τ=λ221τλ222τ.
Thus, it provides the expression for a closed-form solution of the Riccati equation:(61)λ2τ=λ222τ−1λ221τ.
Let compute the last Riccati equation (37) for constant ε,(62)dετdτ=rγ+βTrQTQλ1τ+β¯TrQ¯TQ¯λ2τ,ε0=0.
From Riccati equations (35) and (36), we have(63)dH2τdτ=H1τγγ−12+H2τM,dH1τdτ=−H1τMT+2γRQ−2H2τQTQ,(64)dI2τdτ=I1τγγ−12+I2τM¯,dI1τdτ=−I1τM¯T+2γR¯Q¯−2I2τQ¯TQ¯,(65)dετdτ=rγ+βTrQTQH1−1τH2τ+β¯TrQ¯TQ¯I1−1τI2τ.
From equation (63), we get(66)H2τ=−12dH1τdτ+H1τMT+2γRQQTQ−1,and as well as from equation (64),(67)I2τ=−12dI1τdτ+I1τM¯T+2γR¯Q¯Q¯TQ¯−1,and then, substituting equations (66) and (67) in (65) provides(68)dετdτ=−β2TrH1−1τdH1τdτ+MT+2γRQ−β¯2TrI1−1τdI1τdτ+M¯T+2γR¯Q¯+rγ,ε0=0,(69)dετ=−β2TrdH1τH1τ+MT+2γRQdτ−β¯2TrdI1τI1τ+M¯T+2γR¯Q¯dτ+rγdτ.
By integrating (69), we obtain(70)ετ=−β2TrlogH1τ+MT+2γRQτ−β¯2TrlogI1τ+M¯T+2γR¯Q¯τ+rγτ.


### 3.3. Fast Fourier Transform Method and Characteristic Function

In this section, we consider the fast Fourier transform (FFT) method in the study by Carr and Madan [16] to price European call option with α>0, at time t, strike k=logK, and time to maturity T, given as(71)CtT,K=e−rT−tEXT−K+|ℱt,CtT,K=CtT,k=e−rT−tEexpYt−expk+|ℱt.The modified priceCtα in the study by Carr and Madan [16] is considered, with α=1.1 as a good empirical value for the Heston model. The square integrable function is obtained through the modified price in order to use inverse Fourier transform:(72)CtαT,k=expαkCtT,K.We introduce the Fourier transform of the modified price together with the application of the Fubini integration theorem:(73)ψtαT,θ=∫−∞+∞CtαT,keiθkdk=∫−∞+∞eαkCtT,Keiθkdk=∫−∞+∞eαke−rT−tEexpYt−expk+|ℱteiθkdk=e−rT−t∫−∞+∞expα+iθkEexpYt−expk+|ℱtdk,(74)ψtαT,θ=e−rT−tϕtT,θ−1+αiα+iθα+iθ+1.The call price can be provided through inversion of Fourier transform givenψtαT,θ function having both odd imaginary and even real parts; so from equation (73), we have(75)∫−∞+∞CtαT,keiθkdk=e−rT−tϕtT,θ−1+αiα+iθα+iθ+1,CtαT,k=12π∫−∞∞ψtαT,θe−iθkdθ.Then, we have(76)CtT,K=e−αk2π∫−∞∞ψtαT,θe−iθkdθ,CtT,K=e−αk2π∫−∞∞e−rT−tϕ˜tT,θ−1+αiα+iθα+iθ+1e−iθkdθ.This is a Fourier transform:(77)CtT,K=e−αk2πRe∫−∞∞e−rT−tϕ˜tT,θ−1+αiα+iθα+iθ+1e−iθkdθ.Corollary 1.
LetD be the symmetric matrix, and it is sufficient to find the conditional characteristic function of the double Wishart Σt and Σ¯t given by (Appendix)(78)ϕ˜Σt,Σ¯tD1D2=EexpiTrD1Σt+D2Σ¯t=expTrA1τΣt+A2τΣ¯t+Cτ,where A1τ,A2τ∈Mn, and C∈ℂ verify the following dynamics:(79)dA1τdτ=A1τM+MTA1τ+2A1τQTQA1τ,A10=iD1,dA2τdτ=A2τM¯+M¯TA2τ+2A2τQ¯TQ¯A2τ,A20=iD2,Cτ=∫0τTrβQTQA1τdu+β¯Q¯TQ¯A2τdu.Proposition 3.
The call price under double Wishart is given as(80)CtT,K=e−αk2πRe∫−∞∞e−rT−tϕ˜Σ,Σ¯tϕθe−iθkdθ,where(81)ϕ˜Σ,Σ¯t=eTrA1tΣ0+A2tΣ¯0+Ct+εT−t,ϕθ=1α+iθα+iθ+1.Proof.
Letϕ˜t,T be the characteristic function of the log-price Yt. We have,(82)Yt,T=lnXTXt=lnXT−lnXt=YT−Yt,ϕ˜γ,0t,T=EexpiγYt,T=EEexpiγYT−Yt=Eexp−iγYtEtexpiγYT=Eexp−iγYtexpTrλ1T−tΣt+λ2T−tΣ¯t+iγYt+εT−t=expεT−tEexpTrλ1T−tΣt+λ2T−tΣ¯t=expεT−texpTrA1τΣ0+A2τΣ¯0+Cτ=expεT−texpTrA1tΣ0+A2tΣ¯0+Ctϕ˜γ,0t,T=expTrA1tΣ0+A2tΣ¯0+Ct+εT−t,where Ajt is obtained from equation (79) with τ=t provided Aj0=λjT−t,j=1,2. Hence,(83)CtT,K=e−αk2πRe∫−∞∞e−rT−teTrA1tΣ0+A2tΣ¯0+Ct+εT−tα+iθα+iθ+1e−iθkdθ,where(84)ετ=Trlog1H1τβI1τβ¯−12βMT+β¯M¯T+rγ.


### 3.4. Perturbation Techniques of the Riccati Differential Equations

This section deals with implementation of perturbation techniques since the system of the Riccati differential equations does not allow an analytical closed-form solution. This is due to noncommutativity of matrix multiplication, in order to approximate the call option price.The method procedures retain the affine properties and higher orders comparable to the standard perturbation scheme on partial differential equations that is complex after the first order [9, 17]. Consider the differential Riccati equations of the double Wishart stochastic volatility model and take dimension n=2. The two characteristics orders in perturbation p and q are considered. The solution of Aτ is given in the following form:(85)Aτ=∑i,jpi/2qj/2Ai,jτ.The developments in the perturbed differential equations are done by comparing the coefficients and identifying the terms inp and q to give expected approximations. Let p=M1 and q=M2 are small while νi quantities remain constant. Then, we consider approximation at order one p,q and order two p,q together with the following notations.(86)M=−p00−q=−p1000−q0001=−pM1−qM2,M¯=−p¯00−q¯=−p¯1000−q¯0001=−p¯M1−q¯M2,and also noting that Q is the volatility,(87)Q=pν1M1+qν2M2,Q2=pν12M1+qν22M2,which can be written as(88)Q=pQ1+qQ2.We can now rewrite the Riccati equations as in the form(89)A1′τ=p−A1τM1−M1A1τ+2A1τQ12A1τ+q−A1τM2−M2A1τ+2A1τQ22A1τ,and for(90)λ1′τ=γγ−12I2+p−λ1τM1−M1λ1τ+2λ1τQ12λ1τ+2pγRQ1λ1τ+2qγRQ2λ1τ+q−λ1τM2−M2λ1τ+2λ1τQ22λ1τ,(91)λ2′τ=γγ−12I2+p¯−λ2τM1−M1λ2τ+2λ2τQ¯12λ2τ+2p¯γR¯Q¯1λ2τ+2q¯γR¯Q¯2λ2τ+q¯−λ2τM2−M2λ2τ+2λ2τQ¯22λ2τ.Development of the Riccati functions: by dealing with the development of the Riccati functionsA1,A2,λ1,λ2,C, and ε, the perturbation of the Riccati equations at order two (order 2) is considered to obtain the European call option pricing formula as follows:(92)Akτ=Ak0,0τ+pAk1,0τ+qAk0,1τ+pAk2,0τ+qAk0,2τ+pqAk1,1τ+omaxp,q,k=1,2.Let us determineAk0,0,…,Ak1,1.By identifying the terms in respective orders, we have(93)A10,0τ′=0,A10,0τ∈M2ℂ.SinceA10=iD1,(94)D1=d100d1.Then, we obtain(95)A10,00+pA11,00+qA10,10+pA12,00+qA10,20+pqA11,10=iD1.Then,(96)A10,00=iD1=id1I2,A11,00=A10,10=A12,00=A10,20=A11,10=0,A12,0τ=−2d1i+d1ν12τM1,A10,2τ=−2d1i+d1ν22τM2,A11,1τ′=0,⇒A11,1τ=C,since A11,10=0, for all τ∈ℝ+, C=0.For functionA2τ, we go through the same procedures:(97)A2′τ=p¯−A2τM1−M1A2τ+2A2τQ¯12A2τ+q¯−A2τM2−M2A2τ+2A2τQ¯22A2τ,A2′τ=A20,0τ′+p¯A21,0τ′+q¯A20,1τ′+p¯A22,0τ′+q¯A20,2τ′+p¯q¯A21,1τ′,then,(98)A20,0τ′=0,⇒A20,0τ∈M2ℂ.SinceA20=iD2,(99)D2=d200d2.Then, we have(100)A20,00′+p¯A21,00′+q¯A20,10′+p¯A22,00′+q¯A20,20′+p¯q¯A21,10′=iD2,A20,00=iD2=id2I2,A21,00=A20,10=A22,00=A20,20=A21,10=0.We proceed by obtaining(101)A22,0τ=−2d2i+d2ν¯12τM1,A20,2τ=−2d2i+d2ν¯22τM2,and A21,1τ′=0,A21,1τ=0; since A21,10=0, for all τ∈ℝ+, C=0 is constant. So now, we can get the functions of λ1 and λ2 from the ordinary differential equations (90) and (91) as follows:(102)λ1′τ=γγ−12I2+p−λ1τM1−M1λ1τ+2λ1τQ12λ1τ+2pγRQ1λ1τ+2qγRQ2λ1τ+q−λ1τM2−M2λ1τ+2λ1τQ22λ1τ,λ1τ=λ10,0τ+pλ11,0τ+qλ10,1τ+pλ12,0τ+qλ10,2τ+pqλ11,1τ+omaxp,q,such that(103)λ10,0τ=γγ−12τI2,λ11,0τ=γ2γ−12ν1τ2RM1,λ10,1τ=γ2γ−12ν2τ2RM2,λ12,0τ=−γγ−12τ2M1+γ2γ−126ν12τ3M1+γ3γ−13ν12τ3RM12,λ10,2τ=−γγ−12τ2M2+γ2γ−16ν22τ3M2+γ3γ−13ν22τ3RM22,λ11,1τ=γ3γ−13ν1ν2τ3RM1RM2+RM2RM1.Then, similar forλ2, we have the following approximations:(104)λ20,0τ=γγ−12τI2,λ21,0τ=γ2γ−12ν¯1τ2R¯M1,λ20,1τ=γ2γ−12ν¯2τ2R¯M2,λ22,0τ=−γγ−12τ2M1+γ2γ−126ν¯12τ3M1+γ3γ−13ν¯12τ3R¯M12,λ20,2τ=−γγ−12τ2M2+γ2γ−16ν¯22τ3M2+γ3γ−13ν¯22τ3R¯M22,λ21,1τ=γ3γ−13ν¯1ν¯2τ3R¯M1R¯M2+R¯M2R¯M1.Then, we look at the Riccati differential equation forCτ:(105)C′τ=βTrpν12M1A1τ+qν22M2A1τ+β¯Trpν¯12M1A2τ+qν¯22M2A2τ,C0,0τ=0.Since C0=0,(106)C1,0τ=C0,1τ=0,C2,0τ=iβν12d1τ+iβ¯ν¯12d2τ,C0,2τ=iβν22d1τ+iβ¯ν¯22d2τ,C1,1τ=0.For theε differential equation,(107)ε′τ=rγ+βTrpν12M1λ1τ+qν22M2λ1τ+β¯Trpν¯12M1λ2τ+qν¯22M2λ2τ,we have,(108)ε0,0τ=rγτ,ε0,1τ=0,ε2,0τ=βν12γγ−14τ2+β¯ν¯12γγ−14τ2,ε0,2τ=βν22γγ−14τ2+β¯ν¯22γγ−14τ2,ε1,1τ=0.Substituting back in the perturbation function, we obtain(109)A1τ=A10,0τ+pA11,0τ+qA10,1τ+pA12,0τ+qA10,2τ+pqA11,1τ+omaxp,q,A1τ=id100id1+−2pd1i+d1ν12τ000+000−2qd1i+d1ν22τ+omaxp,q,A1τ=−2pd12ν12τ+id11−2pτ00−2qd12ν22τ+id11−2qτ,A2τ=−2pd22ν¯12τ+id21−2pτ00−2qd22ν¯22τ+id21−2qτ.The variance processesΣ0 and Σ¯0 are given as(110)Σ0=uvvw,Σ¯0=u¯v¯v¯w¯,A1τΣ0=−2pd12ν12τ+id11−2pτ00−2qd12ν22τ+id11−2qτuvvw,A1τΣ0=−2upd12ν12τ+id11−2pτu−2vpd12ν12τ+id11−2pτv−2qvd12ν22τ+id11−2qτv−2qwd12ν22τ+id11−2qτw,A2τΣ¯0=−2u¯pd22ν¯12τ+id21−2pτu¯−2pv¯d22ν¯12τ+id21−2pτv¯−2qv¯d22ν¯22τ+id21−2qτv¯−2qw¯d22ν¯22τ+id21−2qτw¯.Taking the trace,(111)TrA1τΣ0=−2d12upν12+wqν22τ+i1−2pτu+1−2qτwd1,TrA2τΣ¯0=−2d22u¯pν¯12+w¯qν¯22τ+i1−2pτu¯+1−2qτw¯d2.For functionsC and ε,(112)Cτ=pC2,0τ+qC0,2τ+omaxp,q=piβν12d1τ+iβ¯ν¯12d2τ+qiβν22d1τ+iβ¯ν¯22d2τ+omaxp,q=iβd1pν12+qν22τ+iβ¯d2pν¯12+qν¯22τ+omaxp,q.Then,(113)ετ=ε0,0τ+pε2,0τ+qε0,2τ+omaxp,q=rγτ+γγ−14βpν12+qν22τ2+γγ−14β¯pν¯12+qν¯22τ2+omaxp,q.Substituting back in equation (80), we have(114)ϕ˜ΣΣ¯τ=exprγT−t+p−2d12uν12t−2d22u¯ν¯12t+γγ−14βν12T−t2+γγ−14β¯ν¯12T−t2+q−2d12wν22t−2d22w¯ν¯22t+γγ−14βν22T−t2+γγ−14β¯ν¯22T−t2×expi1−2ptud1+1−2qtwd1+1−2ptu¯d2+1−2qtw¯d2+βd1pν12+qν22t+β¯d2pν¯12+qν¯22t+omaxp,q.This can still be written as(115)ϕ˜ΣΣ¯τ=eΔ1teiΔ2t,where(116)Δ1t=rγT−t+p−2d12uν12t−2d22u¯ν¯12t+γγ−14βν12T−t2+γγ−14β¯ν¯12T−t2+q−2d12wν22t−2d22w¯ν¯22t+γγ−14βν22T−t2+γγ−14β¯ν¯22T−t2,Δ2t=1−2ptud1+1−2qtwd1+1−2ptu¯d2+1−2qtw¯d2+βd1pν12+qν22t+β¯d2pν¯12+qν¯22t+omaxp,q,such that the call price is as follows:(117)CtT,K=e−αk2πRe∫−∞∞e−rT−tϕ˜ΣΣ¯tϕθe−iθkdθ≈e−αk2πRe∫−∞∞e−rT−teΔ1teiΔ2te−iθkα+iθα+1+iθdθ,since(118)∫−∞∞e−rT−teΔ1teiΔ2te−iθkα+iθα+1+iθdθ=e−rT−teΔ1teiΔ2t∫−∞∞e−iθkα+iθα+1+iθdθ=e−rT−teΔ1teiΔ2t2πeαk−ekα+1.Finally, we obtain the approximated European call option price as follows:(119)CtT,K≈e−rT−t+Δ1t1−ekcosΔ2t.


## 3.1. Infinitesimal Generator

The log-price process and Wishart processes, together with corresponding pair of two correlated Brownian motions,ZsY,ZsΣ and Z¯sY,Z¯sΣ¯, can be written in the form as in the study by Benabid et al. [9], for easy handling of the complexity of the dynamics:(19)dYt=r−12TrΣt+Σt¯dt+TrΣtdZtY+Σt¯dZ¯tY,dTrΣt=βTrQTQ+2TrMΣtdt+2TrQTQΣtρtdZtY+I−ρt2dZtΣ,dTrΣ¯t=β¯TrQ¯TQ¯+2TrM¯Σ¯tdt+2TrQ¯TQ¯Σ¯tρ¯tdZ¯tY+I−ρ¯t2dZ¯tΣ¯,where(20)dZY,ZΣt=ρt=TrRTQΣtTrΣtTrQTQΣt,dZ¯Y,Z¯Σ¯t=ρ¯t=TrR¯TQ¯Σ¯tTrΣ¯tTrQ¯TQ¯Σ¯t.Proposition 1.
The infinitesimal generator of the double Wishart stochastic volatility model forYt,Σt,Σ¯t is given by(21)ℒY,Σ,Σ¯=r−TrΣ+Σ¯2∂∂y+TrΣ+Σ¯2∂2∂y2+βTrQTQ+2TrMΣ∂∂Σ+2TrΣ∂∂ΣQTQ∂∂Σ+β¯TrQ¯TQ¯+2TrM¯Σ¯∂∂Σ¯+2TrΣ¯∂∂Σ¯Q¯TQ¯∂∂Σ¯+2TrΣRQ∂∂Σ∂∂y+2TrΣ¯R¯Q¯∂∂Σ¯∂∂y.Proof.
The infinitesimal generator has a nontrivial term which arises from covariationdΣθij,Y corresponding to the coefficients of the term(22)∂2∂xθ;ij∂y=∂∂xθ;ij∂∂y,i,j=1,…,n,θ=1,2.
LetVθ;t≔Σθ;t be the square root matrix, such that(23)Σθ;tij=∑t=1nVθ;tilVθ;tlj=∑t=1nVθ;tilVθ;tjl.
Due to the fact thatVθ;t is symmetric, we find the covariation terms matching with ∂2/∂xθ;ij∂y coefficients:(24)dΣθij,Y=Et∑i,k=1nVθ;tildWlkθQkjθ+∑l,k=1nVθ;tjldWlkθQkiθ∑l,k,h=1nVθ;tlkdWkhθRlhθ=∑l,k,h=1nVθ;tilQkjθ+Vθ;tjlQkiθVθ;thlRhkθdt=∑k,h=1n∑l=1nVθ;tilVθ;thlQkjθ+∑l=1nVθ;tjlVθ;thlQkiθRhkθdt=∑k,h=1nΣθ;tihQkjθ+Σθ;tjhQkiθRhkθdt.
This gives the corresponding coefficients of the term since(25)2TrΣθRθQθDθ∂∂y=2∑i,j,k,h=1nDθijΣθjhRhkθQkiθ∂∂y.
The notation whenθ=1,Σ1=Σ,R1=R,Q1=Q,D1=∂/∂xij for θ=2,Σ2=Σ¯,R2=R¯,Q2=Q¯,D2=∂/∂x¯ij.


## 3.2. The Laplace Transform of Asset Returns

Following Duffie et al. [15] and Da Fonseca et al. [3], in order to solve the European options pricing problem, we consider the Laplace transform of the process (17). The Laplace transform of the Wishart process is exponentially affine as in the study by Bru [11]; hence, the conditional moment generating function of the log-asset returns is the exponential of an affine combination of Y and Wishart components, so we provide for the deterministic functions λ1t,λ2t∈Mn, and δt,εt∈ℝ, as parameters of the Laplace transform:(26)ψγ,tτ=EeγYt+τ=expTrλ1τΣt+Trλ2τΣ¯t+δτYt+ετ,where γ∈ℝ and Et denote the conditional expected value with respect to the probability measure.Using Feynman–Kac argument helps us to get matrix Riccati equations:(27)∂ψγ,t∂τ=ℒY,Σ,Σ¯ψγ,t,ψγ,t0=expγYt.Proposition 2.
The Laplace transform of the asset returns(28)ψγ,tτ=expTrλ1τΣt+λ2τΣ¯t+δτYt+ετ,where λ1,λ2,andε are the solutions of the differential equations through linearization methods:(29)λ1′τ=Mλ1τ+MT+2γRQλ1τ+2λ1τQTQλ1τ+γγ−12In,λ2′τ=M¯λ2τ+M¯T+2γR¯Q¯λ2τ+2λ2τQ¯TQ¯λ2τ+γγ−12In,ε′τ=rγ+βTrQTQλ1τ+β¯TrQ¯TQ¯λ2τ,with boundary conditions λ10=0,λ20=0,ε0=0,andδτ=γ=C0.The solutionsλ1,λ2,andε are obtained as follows:(30)λ1τ=H1τ−1H2τ,λ2τ=I1−1τI2τ,ετ=−β2TrlogH1τ+MT+2γRQτ−β¯2TrlogI1τ+M¯T+2γR¯Q¯τ+rγτ,with(31)H2τH1τ=H20H10expτM−2QTQγγ−12In−MT+2γRQ,I2τI1τ=I20I10expτM¯−2Q¯TQ¯γγ−12In−M¯T+2γR¯Q¯.Proof.
From equation (27) and the Laplace transform (26), we solve the problem by considering Proposition 1 for the double Wishart stochastic volatility model as follows:(32)∂ψγ,t∂τ=r−TrΣ+Σ¯2∂ψγ,t∂y+TrΣ+Σ¯2∂2ψγ,t∂y2+βTrQTQ+2TrMΣ∂ψγ,t∂Σ+2TrΣ∂∂ΣQTQ∂∂Σψγ,t+β¯TrQ¯TQ¯+2TrM¯Σ¯∂ψγ,t∂Σ¯+2TrΣ¯∂∂Σ¯Q¯TQ¯∂∂Σ¯ψγ,t+2TrΣRQ∂∂Σ∂ψγ,t∂y+2TrΣ¯R¯Q¯∂∂Σ¯∂ψγ,t∂y.
From equation (27), we have(33)λ10=0,λ20=0,δ0=γ,ε0=0.
Then,(34)∂ψγ,tτ∂τ=Trdλ1τdτΣt+dλ2τ∂τΣ¯t+ddτδtYt+ddτετψγ,tτ=r−Σ+Σ¯2δτψγ,t+TrΣ+Σ¯2δ2τψγ,t+βTrQTQ+2TrMΣλ1τ+β¯TrQ¯TQ¯+2TrM¯Σ¯λ2τψγ,t+2TrΣλ1τQTQλ1τψγ,tτ+2TrΣ¯λ1τQ¯TQ¯λ1τψγ,t+2TrΣRQδτλ1τψγ,t+2TrΣ¯R¯Q¯δτλ1τψγ,t,−Trλ´1τΣt+λ´2τΣ¯t+δ´τYt+ε´τ+r−TrΣ+Σ¯2δτ+TrΣ+Σ¯2δ2τ+TrβQTQ+MΣ+ΣMTλ1τ+β¯Q¯TQ¯+M¯Σ¯+Σ¯M¯Tλ2τ+2Σλ1τQTQλ1τ+2Σ¯λ2τQ¯TQ¯λ2τ+2ΣRQδτλ1τ+2Σ¯R¯Q¯δτλ2τ=0,−Trddτλ1τΣ−12TrΣδτ+12TrΣδ2τ+TrMΣ+ΣMTλ1τ+2Σλ1τQTQλ1τ+2ΣRQδτλ1τ=0,−Trddτλ2τΣ¯−12TrΣ¯δτ+12TrΣ¯δ2τ+TrM¯Σ¯+Σ¯M¯Tλ2τ+2Σ¯λ2τQ¯TQ¯λ2τ+2Σ¯R¯Q¯δτλ2τ=0,−δ´τYt−ε´τ+μδτ+TrβQTQλ1τ+Trβ¯Q¯TQ¯λ2τ=0.
Now, we need to identify the coefficients for above equations and obtain the matrix Riccati ordinary differential equations:(35)ddτλ1τ=Mλ1τ+MT+2γRQλ1τ+2λ1τQTQλ1τ+γγ−12In,λ10=0,(36)ddτλ2τ=M¯λ2τ+M¯T+2γR¯Q¯λ2τ+2λ2τQ¯TQ¯λ2τ+γγ−12In,λ20=0.
For the constantε, the matrix Riccati ordinary differential is obtained as follows:(37)−δ′τYt−ε′τ+rδτ+TrβQTQλ1τ+Trβ¯Q¯TQ¯λ2τ=0,δτ=C0=γ,sinceδ′τ=0so⇒δτ=γ,dετdτ=rγ+βTrQTQλ1τ+β¯TrQ¯TQ¯λ2τ,ε0=0.
Then,ετ is obtained by integrating directly(38)ετ=∫0τrγds+∫0τTrβQTQλ1s+β¯Q¯TQ¯λ2sds,with λ1,λ2∈Mnℝ, and δτ,ετ∈ℝ.
The matrix Riccati equations above are linearized to provide the closed-form solution, following the techniques of Benabid et al. [9] and Da Fonseca et al. [3]. So from equations (35) and (36), let(39)λ1τ=H1τ−1H2τ,where H1τ∈GLnℝ, H2τ∈Mnℝ, and therefore,(40)ddτH1τλ1τ=dH1τdτλ1τ+H1τdλ1τdτ,H1τdλ1τdτ=ddτH1τλ1τ−dH1τdτλ1τ,H1τλ1τM+H1τMT+2γRQλ1τ+2H1τλ1τQTQλ1τ+H1τγγ−12=H1τdλ1τdτ,H2τM+H1τMT+2γRQλ1τ+2H2τQTQλ2τ+H1τγγ−12=H1τdλ1τdτ.
SinceH2τ=H1τλ1τ, now we have(41)dH2τdτ=dH1τdτλ1τ+H1τdλ1τdτ.
This implies that(42)H1τdλ1τdτ=dH2τdτ−dH1τdτλ1τ.
Then, we get expressions(43)dH2τdτ−dH1τdτλ1τ=H2τM+H1τMT+2γRQλ1τ+2H2τQTQλ1τ+H1τγγ−12=H1τγγ−12+H2τM+H1τMT+2γRQ+2H2τQTQλ1τ,(44)dH2τdτ=H1τγγ−12+H2τM,(45)dH1τdτ=−H1τMT+2γRQ−2H2τQTQ,(46)ddτH2τH1τ=H2τH1τM−2QTQγγ−12In−MT+2γRQ.
The solution of differential equation (46) above is obtained through exponentiation as follows:(47)H2τH1τ=H20H10eQτ,with conditions H10=KIn and H1−10=K−1In and for K=1,H1−10=In, such that(48)H2τH1τ=λ10IneτQ,where(49)eτQ=λ111τλ112τλ121τλ122τ.
Then,(50)H2τH1τ=λ10Inλ111τλ112τλ121τλ122τ=λ10λ111τ+λ121τλ10λ112τ+λ122τ,and since λ10=0, then(51)H2τH1τ=λ121τλ122τ.
In conclusion, we get(52)λ1τ=λ122τ−1λ121τ.
This represents the closed-form solution of the Riccati equation (35).
Now, we look at the solution of the second Riccati equation (36). Let(53)λ2τ=I1−1τI2τ,then,(54)dI2τdτ=I1τγγ−12+I2τM¯,(55)dI1τdτ=−I1τM¯T+2γR¯Q¯−2I2τQ¯TQ¯,(56)ddτI2τI1τ=I2τI1τM¯−2Q¯TQ¯γγ−12In−M¯T+2γR¯Q¯,such that(57)I2τI1τ=I20I10eQ¯τ.
From equation (40), we go through the same procedures as above, with conditions(58)I20=I10λ20=λ20,I10=K¯I=K¯=I,therefore,(59)I2τI1τ=λ20Inλ211τλ212τλ221τλ222τ,I2τI1τ=λ20λ211τ+λ221τλ20λ212τ+λ222τ,since λ20=0, then(60)I2τI1τ=λ221τλ222τ.
Thus, it provides the expression for a closed-form solution of the Riccati equation:(61)λ2τ=λ222τ−1λ221τ.
Let compute the last Riccati equation (37) for constant ε,(62)dετdτ=rγ+βTrQTQλ1τ+β¯TrQ¯TQ¯λ2τ,ε0=0.
From Riccati equations (35) and (36), we have(63)dH2τdτ=H1τγγ−12+H2τM,dH1τdτ=−H1τMT+2γRQ−2H2τQTQ,(64)dI2τdτ=I1τγγ−12+I2τM¯,dI1τdτ=−I1τM¯T+2γR¯Q¯−2I2τQ¯TQ¯,(65)dετdτ=rγ+βTrQTQH1−1τH2τ+β¯TrQ¯TQ¯I1−1τI2τ.
From equation (63), we get(66)H2τ=−12dH1τdτ+H1τMT+2γRQQTQ−1,and as well as from equation (64),(67)I2τ=−12dI1τdτ+I1τM¯T+2γR¯Q¯Q¯TQ¯−1,and then, substituting equations (66) and (67) in (65) provides(68)dετdτ=−β2TrH1−1τdH1τdτ+MT+2γRQ−β¯2TrI1−1τdI1τdτ+M¯T+2γR¯Q¯+rγ,ε0=0,(69)dετ=−β2TrdH1τH1τ+MT+2γRQdτ−β¯2TrdI1τI1τ+M¯T+2γR¯Q¯dτ+rγdτ.
By integrating (69), we obtain(70)ετ=−β2TrlogH1τ+MT+2γRQτ−β¯2TrlogI1τ+M¯T+2γR¯Q¯τ+rγτ.


## 3.3. Fast Fourier Transform Method and Characteristic Function

In this section, we consider the fast Fourier transform (FFT) method in the study by Carr and Madan [16] to price European call option with α>0, at time t, strike k=logK, and time to maturity T, given as(71)CtT,K=e−rT−tEXT−K+|ℱt,CtT,K=CtT,k=e−rT−tEexpYt−expk+|ℱt.The modified priceCtα in the study by Carr and Madan [16] is considered, with α=1.1 as a good empirical value for the Heston model. The square integrable function is obtained through the modified price in order to use inverse Fourier transform:(72)CtαT,k=expαkCtT,K.We introduce the Fourier transform of the modified price together with the application of the Fubini integration theorem:(73)ψtαT,θ=∫−∞+∞CtαT,keiθkdk=∫−∞+∞eαkCtT,Keiθkdk=∫−∞+∞eαke−rT−tEexpYt−expk+|ℱteiθkdk=e−rT−t∫−∞+∞expα+iθkEexpYt−expk+|ℱtdk,(74)ψtαT,θ=e−rT−tϕtT,θ−1+αiα+iθα+iθ+1.The call price can be provided through inversion of Fourier transform givenψtαT,θ function having both odd imaginary and even real parts; so from equation (73), we have(75)∫−∞+∞CtαT,keiθkdk=e−rT−tϕtT,θ−1+αiα+iθα+iθ+1,CtαT,k=12π∫−∞∞ψtαT,θe−iθkdθ.Then, we have(76)CtT,K=e−αk2π∫−∞∞ψtαT,θe−iθkdθ,CtT,K=e−αk2π∫−∞∞e−rT−tϕ˜tT,θ−1+αiα+iθα+iθ+1e−iθkdθ.This is a Fourier transform:(77)CtT,K=e−αk2πRe∫−∞∞e−rT−tϕ˜tT,θ−1+αiα+iθα+iθ+1e−iθkdθ.Corollary 1.
LetD be the symmetric matrix, and it is sufficient to find the conditional characteristic function of the double Wishart Σt and Σ¯t given by (Appendix)(78)ϕ˜Σt,Σ¯tD1D2=EexpiTrD1Σt+D2Σ¯t=expTrA1τΣt+A2τΣ¯t+Cτ,where A1τ,A2τ∈Mn, and C∈ℂ verify the following dynamics:(79)dA1τdτ=A1τM+MTA1τ+2A1τQTQA1τ,A10=iD1,dA2τdτ=A2τM¯+M¯TA2τ+2A2τQ¯TQ¯A2τ,A20=iD2,Cτ=∫0τTrβQTQA1τdu+β¯Q¯TQ¯A2τdu.Proposition 3.
The call price under double Wishart is given as(80)CtT,K=e−αk2πRe∫−∞∞e−rT−tϕ˜Σ,Σ¯tϕθe−iθkdθ,where(81)ϕ˜Σ,Σ¯t=eTrA1tΣ0+A2tΣ¯0+Ct+εT−t,ϕθ=1α+iθα+iθ+1.Proof.
Letϕ˜t,T be the characteristic function of the log-price Yt. We have,(82)Yt,T=lnXTXt=lnXT−lnXt=YT−Yt,ϕ˜γ,0t,T=EexpiγYt,T=EEexpiγYT−Yt=Eexp−iγYtEtexpiγYT=Eexp−iγYtexpTrλ1T−tΣt+λ2T−tΣ¯t+iγYt+εT−t=expεT−tEexpTrλ1T−tΣt+λ2T−tΣ¯t=expεT−texpTrA1τΣ0+A2τΣ¯0+Cτ=expεT−texpTrA1tΣ0+A2tΣ¯0+Ctϕ˜γ,0t,T=expTrA1tΣ0+A2tΣ¯0+Ct+εT−t,where Ajt is obtained from equation (79) with τ=t provided Aj0=λjT−t,j=1,2. Hence,(83)CtT,K=e−αk2πRe∫−∞∞e−rT−teTrA1tΣ0+A2tΣ¯0+Ct+εT−tα+iθα+iθ+1e−iθkdθ,where(84)ετ=Trlog1H1τβI1τβ¯−12βMT+β¯M¯T+rγ.


## 3.4. Perturbation Techniques of the Riccati Differential Equations

This section deals with implementation of perturbation techniques since the system of the Riccati differential equations does not allow an analytical closed-form solution. This is due to noncommutativity of matrix multiplication, in order to approximate the call option price.The method procedures retain the affine properties and higher orders comparable to the standard perturbation scheme on partial differential equations that is complex after the first order [9, 17]. Consider the differential Riccati equations of the double Wishart stochastic volatility model and take dimension n=2. The two characteristics orders in perturbation p and q are considered. The solution of Aτ is given in the following form:(85)Aτ=∑i,jpi/2qj/2Ai,jτ.The developments in the perturbed differential equations are done by comparing the coefficients and identifying the terms inp and q to give expected approximations. Let p=M1 and q=M2 are small while νi quantities remain constant. Then, we consider approximation at order one p,q and order two p,q together with the following notations.(86)M=−p00−q=−p1000−q0001=−pM1−qM2,M¯=−p¯00−q¯=−p¯1000−q¯0001=−p¯M1−q¯M2,and also noting that Q is the volatility,(87)Q=pν1M1+qν2M2,Q2=pν12M1+qν22M2,which can be written as(88)Q=pQ1+qQ2.We can now rewrite the Riccati equations as in the form(89)A1′τ=p−A1τM1−M1A1τ+2A1τQ12A1τ+q−A1τM2−M2A1τ+2A1τQ22A1τ,and for(90)λ1′τ=γγ−12I2+p−λ1τM1−M1λ1τ+2λ1τQ12λ1τ+2pγRQ1λ1τ+2qγRQ2λ1τ+q−λ1τM2−M2λ1τ+2λ1τQ22λ1τ,(91)λ2′τ=γγ−12I2+p¯−λ2τM1−M1λ2τ+2λ2τQ¯12λ2τ+2p¯γR¯Q¯1λ2τ+2q¯γR¯Q¯2λ2τ+q¯−λ2τM2−M2λ2τ+2λ2τQ¯22λ2τ.Development of the Riccati functions: by dealing with the development of the Riccati functionsA1,A2,λ1,λ2,C, and ε, the perturbation of the Riccati equations at order two (order 2) is considered to obtain the European call option pricing formula as follows:(92)Akτ=Ak0,0τ+pAk1,0τ+qAk0,1τ+pAk2,0τ+qAk0,2τ+pqAk1,1τ+omaxp,q,k=1,2.Let us determineAk0,0,…,Ak1,1.By identifying the terms in respective orders, we have(93)A10,0τ′=0,A10,0τ∈M2ℂ.SinceA10=iD1,(94)D1=d100d1.Then, we obtain(95)A10,00+pA11,00+qA10,10+pA12,00+qA10,20+pqA11,10=iD1.Then,(96)A10,00=iD1=id1I2,A11,00=A10,10=A12,00=A10,20=A11,10=0,A12,0τ=−2d1i+d1ν12τM1,A10,2τ=−2d1i+d1ν22τM2,A11,1τ′=0,⇒A11,1τ=C,since A11,10=0, for all τ∈ℝ+, C=0.For functionA2τ, we go through the same procedures:(97)A2′τ=p¯−A2τM1−M1A2τ+2A2τQ¯12A2τ+q¯−A2τM2−M2A2τ+2A2τQ¯22A2τ,A2′τ=A20,0τ′+p¯A21,0τ′+q¯A20,1τ′+p¯A22,0τ′+q¯A20,2τ′+p¯q¯A21,1τ′,then,(98)A20,0τ′=0,⇒A20,0τ∈M2ℂ.SinceA20=iD2,(99)D2=d200d2.Then, we have(100)A20,00′+p¯A21,00′+q¯A20,10′+p¯A22,00′+q¯A20,20′+p¯q¯A21,10′=iD2,A20,00=iD2=id2I2,A21,00=A20,10=A22,00=A20,20=A21,10=0.We proceed by obtaining(101)A22,0τ=−2d2i+d2ν¯12τM1,A20,2τ=−2d2i+d2ν¯22τM2,and A21,1τ′=0,A21,1τ=0; since A21,10=0, for all τ∈ℝ+, C=0 is constant. So now, we can get the functions of λ1 and λ2 from the ordinary differential equations (90) and (91) as follows:(102)λ1′τ=γγ−12I2+p−λ1τM1−M1λ1τ+2λ1τQ12λ1τ+2pγRQ1λ1τ+2qγRQ2λ1τ+q−λ1τM2−M2λ1τ+2λ1τQ22λ1τ,λ1τ=λ10,0τ+pλ11,0τ+qλ10,1τ+pλ12,0τ+qλ10,2τ+pqλ11,1τ+omaxp,q,such that(103)λ10,0τ=γγ−12τI2,λ11,0τ=γ2γ−12ν1τ2RM1,λ10,1τ=γ2γ−12ν2τ2RM2,λ12,0τ=−γγ−12τ2M1+γ2γ−126ν12τ3M1+γ3γ−13ν12τ3RM12,λ10,2τ=−γγ−12τ2M2+γ2γ−16ν22τ3M2+γ3γ−13ν22τ3RM22,λ11,1τ=γ3γ−13ν1ν2τ3RM1RM2+RM2RM1.Then, similar forλ2, we have the following approximations:(104)λ20,0τ=γγ−12τI2,λ21,0τ=γ2γ−12ν¯1τ2R¯M1,λ20,1τ=γ2γ−12ν¯2τ2R¯M2,λ22,0τ=−γγ−12τ2M1+γ2γ−126ν¯12τ3M1+γ3γ−13ν¯12τ3R¯M12,λ20,2τ=−γγ−12τ2M2+γ2γ−16ν¯22τ3M2+γ3γ−13ν¯22τ3R¯M22,λ21,1τ=γ3γ−13ν¯1ν¯2τ3R¯M1R¯M2+R¯M2R¯M1.Then, we look at the Riccati differential equation forCτ:(105)C′τ=βTrpν12M1A1τ+qν22M2A1τ+β¯Trpν¯12M1A2τ+qν¯22M2A2τ,C0,0τ=0.Since C0=0,(106)C1,0τ=C0,1τ=0,C2,0τ=iβν12d1τ+iβ¯ν¯12d2τ,C0,2τ=iβν22d1τ+iβ¯ν¯22d2τ,C1,1τ=0.For theε differential equation,(107)ε′τ=rγ+βTrpν12M1λ1τ+qν22M2λ1τ+β¯Trpν¯12M1λ2τ+qν¯22M2λ2τ,we have,(108)ε0,0τ=rγτ,ε0,1τ=0,ε2,0τ=βν12γγ−14τ2+β¯ν¯12γγ−14τ2,ε0,2τ=βν22γγ−14τ2+β¯ν¯22γγ−14τ2,ε1,1τ=0.Substituting back in the perturbation function, we obtain(109)A1τ=A10,0τ+pA11,0τ+qA10,1τ+pA12,0τ+qA10,2τ+pqA11,1τ+omaxp,q,A1τ=id100id1+−2pd1i+d1ν12τ000+000−2qd1i+d1ν22τ+omaxp,q,A1τ=−2pd12ν12τ+id11−2pτ00−2qd12ν22τ+id11−2qτ,A2τ=−2pd22ν¯12τ+id21−2pτ00−2qd22ν¯22τ+id21−2qτ.The variance processesΣ0 and Σ¯0 are given as(110)Σ0=uvvw,Σ¯0=u¯v¯v¯w¯,A1τΣ0=−2pd12ν12τ+id11−2pτ00−2qd12ν22τ+id11−2qτuvvw,A1τΣ0=−2upd12ν12τ+id11−2pτu−2vpd12ν12τ+id11−2pτv−2qvd12ν22τ+id11−2qτv−2qwd12ν22τ+id11−2qτw,A2τΣ¯0=−2u¯pd22ν¯12τ+id21−2pτu¯−2pv¯d22ν¯12τ+id21−2pτv¯−2qv¯d22ν¯22τ+id21−2qτv¯−2qw¯d22ν¯22τ+id21−2qτw¯.Taking the trace,(111)TrA1τΣ0=−2d12upν12+wqν22τ+i1−2pτu+1−2qτwd1,TrA2τΣ¯0=−2d22u¯pν¯12+w¯qν¯22τ+i1−2pτu¯+1−2qτw¯d2.For functionsC and ε,(112)Cτ=pC2,0τ+qC0,2τ+omaxp,q=piβν12d1τ+iβ¯ν¯12d2τ+qiβν22d1τ+iβ¯ν¯22d2τ+omaxp,q=iβd1pν12+qν22τ+iβ¯d2pν¯12+qν¯22τ+omaxp,q.Then,(113)ετ=ε0,0τ+pε2,0τ+qε0,2τ+omaxp,q=rγτ+γγ−14βpν12+qν22τ2+γγ−14β¯pν¯12+qν¯22τ2+omaxp,q.Substituting back in equation (80), we have(114)ϕ˜ΣΣ¯τ=exprγT−t+p−2d12uν12t−2d22u¯ν¯12t+γγ−14βν12T−t2+γγ−14β¯ν¯12T−t2+q−2d12wν22t−2d22w¯ν¯22t+γγ−14βν22T−t2+γγ−14β¯ν¯22T−t2×expi1−2ptud1+1−2qtwd1+1−2ptu¯d2+1−2qtw¯d2+βd1pν12+qν22t+β¯d2pν¯12+qν¯22t+omaxp,q.This can still be written as(115)ϕ˜ΣΣ¯τ=eΔ1teiΔ2t,where(116)Δ1t=rγT−t+p−2d12uν12t−2d22u¯ν¯12t+γγ−14βν12T−t2+γγ−14β¯ν¯12T−t2+q−2d12wν22t−2d22w¯ν¯22t+γγ−14βν22T−t2+γγ−14β¯ν¯22T−t2,Δ2t=1−2ptud1+1−2qtwd1+1−2ptu¯d2+1−2qtw¯d2+βd1pν12+qν22t+β¯d2pν¯12+qν¯22t+omaxp,q,such that the call price is as follows:(117)CtT,K=e−αk2πRe∫−∞∞e−rT−tϕ˜ΣΣ¯tϕθe−iθkdθ≈e−αk2πRe∫−∞∞e−rT−teΔ1teiΔ2te−iθkα+iθα+1+iθdθ,since(118)∫−∞∞e−rT−teΔ1teiΔ2te−iθkα+iθα+1+iθdθ=e−rT−teΔ1teiΔ2t∫−∞∞e−iθkα+iθα+1+iθdθ=e−rT−teΔ1teiΔ2t2πeαk−ekα+1.Finally, we obtain the approximated European call option price as follows:(119)CtT,K≈e−rT−t+Δ1t1−ekcosΔ2t.


## 4. Numerical Illustration

In this section, we deal with illustrative numerical examples to examine the implications of volatility in option pricing under the double Wishart stochastic volatility model.


### 4.1. Comparison between the Double Wishart Model Price and the Market Price

The volatility specification of the double Wishart model remarkably makes it flexible to produce price predictions which exhibit a similar price behavior with the market price. We quote market data drawn from QQQ (a fund by Invesco that tracks the performance of the stocks listed under the NASDAQ Index) options, April 2020.Let us consider the estimated parameter values from the market data; the variance matricesΣ and Σ¯, with volatility of volatility matrices Q and Q¯, are considered since they are very important parameters in a stochastic volatility model.(120)Σ0=0.00008000.000095,Σ¯0=0.00012000.0001,Q=0.10000.23,Q¯=0.20000.25.While, the strike price isK=189 and interest rate is r=0.05.Figure1 illustrates the effect of implied volatilities on the call price under the double Wishart model. Given the parameter values, we observe that the price predictions under the double Wishart model exhibits the same behavior with respect to the market price in short maturity time. This depends also on the choice of the model parameters, such as γ=0.5,β=4, β¯=3, d1=0.45,andd2=0.479. This example proves that the double Wishart volatility model has greater flexibility.Figure 1
The price predictions at 6 months maturity time and the effect of model parameters on the call option price.Figure2 shows that a slight change in model parameters leads to the changes in the model call option price predictions. We observe that as time tends to maturity time of six months, d2=0.478, the model price movement behaves almost the same with the market price. The parameters γ=0.5,β=4,β¯=3, and d1=0.45 are maintained.Figure 2
The call price predictions under the double Wishart model and market prices at 6 months maturity time.Figure3 demonstrates the effect of parameters on the call European option price under the double Wishart model. The parameters γ=1,β=4,β¯=3, d1=0.45, and d2=0.4968. The model shows price predictions exhibiting the same market price behavior, in long maturity time of 1 year. We note that the model parameters influence the prices greatly and help long position holders not to experience arbitrage profits.Figure 3
The double Wishart model price predictions and market price at 1 year maturity time.


## 4.1. Comparison between the Double Wishart Model Price and the Market Price

The volatility specification of the double Wishart model remarkably makes it flexible to produce price predictions which exhibit a similar price behavior with the market price. We quote market data drawn from QQQ (a fund by Invesco that tracks the performance of the stocks listed under the NASDAQ Index) options, April 2020.Let us consider the estimated parameter values from the market data; the variance matricesΣ and Σ¯, with volatility of volatility matrices Q and Q¯, are considered since they are very important parameters in a stochastic volatility model.(120)Σ0=0.00008000.000095,Σ¯0=0.00012000.0001,Q=0.10000.23,Q¯=0.20000.25.While, the strike price isK=189 and interest rate is r=0.05.Figure1 illustrates the effect of implied volatilities on the call price under the double Wishart model. Given the parameter values, we observe that the price predictions under the double Wishart model exhibits the same behavior with respect to the market price in short maturity time. This depends also on the choice of the model parameters, such as γ=0.5,β=4, β¯=3, d1=0.45,andd2=0.479. This example proves that the double Wishart volatility model has greater flexibility.Figure 1
The price predictions at 6 months maturity time and the effect of model parameters on the call option price.Figure2 shows that a slight change in model parameters leads to the changes in the model call option price predictions. We observe that as time tends to maturity time of six months, d2=0.478, the model price movement behaves almost the same with the market price. The parameters γ=0.5,β=4,β¯=3, and d1=0.45 are maintained.Figure 2
The call price predictions under the double Wishart model and market prices at 6 months maturity time.Figure3 demonstrates the effect of parameters on the call European option price under the double Wishart model. The parameters γ=1,β=4,β¯=3, d1=0.45, and d2=0.4968. The model shows price predictions exhibiting the same market price behavior, in long maturity time of 1 year. We note that the model parameters influence the prices greatly and help long position holders not to experience arbitrage profits.Figure 3
The double Wishart model price predictions and market price at 1 year maturity time.


## 5. Conclusion

The generalization of the Heston model into a multifactor form with two dependence matrices efficiently solves the problem of pricing European call option and financial market data fitting for short or long maturities. The new Wishart affine model retains the remarkable property of analytical tractability, and this allows the model to form a closed-form solution of the conditional characteristic function for the call option price expression, which is obtained through Fourier transforms together with perturbation methods.The numerical results show that the double Wishart stochastic volatility model produces price predictions similar to the market prices. Notably the effect of parameters on our model makes it flexible enough for short or long maturities. We recommend future work on discretization schemes and also to investigate the behavior of nondiagonal matrix components in the model.



---

*Source: 7411885-2021-07-12.xml*