# Characteristic Function for Wishart Heston Model

I don't know if this is the right place (at most they will close the post). Anyway, I am trying to implement the characteristic function of the Heston Wishart Stochastic Volatility model illustrated in the work of Da Fonseca, Grasselli and Tebaldi [2008] in order to use Fourier transform methods

The characteristic function is

$$\Psi_{WMSV}(u;\tau) = \exp\Big\{ Tr[A(\tau)\Sigma_{t}] + i\gamma Y_{t} + C(\tau) \Big\}$$

where $$Y_{t} := \log S_{t}$$ and $$\Sigma_{t}$$ follows the wishart SDE.

Using the linearization procedure (as explained in the paper) we could solve the matrix Riccati ODE $$A(\tau)$$. The solution for $$A(\tau)$$ is given by equation (12)

$$$$A(\tau) = F(\tau)^{-1}G(\tau) = A_{22}(\tau)^{-1}A_{21}(\tau)$$$$ where $$F(\tau), G(\tau)$$ solve a system of $$(2n)$$ linear ODE whose solution is given by $$$$\begin{pmatrix} A_{11}(\tau) & A_{12}(\tau) \\ A_{21}(\tau) & A_{22}(\tau) \end{pmatrix} = \exp\left[ \tau \begin{pmatrix} M & -2Q^TQ \\ \frac{\gamma(\gamma-1)}{2}\mathbb{I}_{n} &-(M^T + 2\gamma RQ) \end{pmatrix} \right]$$$$ instead the solution for $$C(\tau)$$ is given by $$$$C(\tau) = -\frac{\beta}{2} Tr\Big[ \log F(\tau) + \tau(M^T + 2i\gamma RQ) \Big] + i\gamma(r-q)\tau$$$$

where $$F(\tau) = A_{22}(\tau)$$. Remark that Da Fonsenca, Grasselli and Tebaldi derived the moment generating function of log-price. So in order to obtain the characteristic function we have to valuate $$\gamma := i\gamma$$ where $$\gamma \in \mathbb{R}$$.

Now, I tried to implement the characteristic function of WMSV model in MatLab but i have some problems. I poste a sketch of code

clear variables; close all; clc

% -----------------------------
% Parameters
% -----------------------------
S0 = 100; % Initial stock price
r = 0.05; % Risk free rate
q = 0; % Dividend yield
t = 1.0; % Time to maturity
beta = 3; % Role of Feller condition

M = [-3, 0.0; 0.0, -3];
R = [-0.7, 0.0; 0.0, -0.7];
Q = [0.25, 0.0; 0.0, 0.25];
Sigma0 = [0.01, 0.0; 0.0, 0.01];  % Initial variance

% u = linspace(-10,10,1000);
u = 3.5; % for example

cf = cf_wmsv(u,t,r,q,M,R,Q,Sigma0,S0,beta);

% ----------------------------------------------------------
% Here compute the WMSV model characteristic function
% ----------------------------------------------------------
function resu = cf_wmsv(u,t,r,q,M,R,Q,Sigma0,S0,beta)
i = complex(0,1);
% Exponential of the matrix (14) in Da Fonseca et al. [2008]
MATEXP = expm(t*[M -2.0*(Q'*Q); 0.5*i*u*(i*u - 1)*eye(2) -1.0*(M' + 2.0*i*u*R*Q)]);
A11 = [MATEXP(1:2,1:2)];
A12 = [MATEXP(1:2,3:4)];
A21 = [MATEXP(3:4,1:2)];
A22 = [MATEXP(3:4,3:4)];
% Computation of matrix function A(tau) equation (12-15) in Da Fonseca et al.
A = A22\A21;
% Computation of scalar function C(tau)
C = -0.5*beta*trace(log(A22) + t*(M' + 2*R*Q)) + i*u*(r-q);
% Characteristic Function for WMSV model equation (5) in Da Fonseca et al.
resu = exp(trace(A*Sigma0) + C + i*u*log(S0));
end



In the code the "role" of $$\gamma$$ is replaced by $$u$$. The output should be a complex double vector $$1\times \mathtt{length(u)}$$. If I consider a scalar input $$u$$ i have no problem but if I consider (Even if I replace the matrix product $$*$$ with the element-by-element product $$.*$$) I do not get the desired result.

I considered that I could use a for loop and evaluate the complex 4*4 matrix element by element, but I do not know how to adjust the remaining. Can anyone give me a suggestion or solution?

I have solved the problem myself. In summary, one has to evaluate the complex matrix for each argument of the Fourier/Laplace transform $$u_1, \dots, u_{1000}$$. This cannot be done with MatLab's element-by-element product. So I used a for loop that would multiply for each argument the matrix by the scalar $$u_1, \dots, u_{1000}$$.

Here is the revised code:

function resu = chf_wmsv(u,t,r,q,M,R,Q,Sigma0,S0,beta)
i = complex(0,1);
resu = zeros(1,length(u));
for idx = 1:length(u)
% Exponential of the matrix (14) in Da Fonseca et al. [2008]
w = u(1,idx);
MATEXP = expm(t*[M, -2.0*(Q'*Q); ...
0.5*i*w*(i*w - 1)*eye(2), -1.0*(M' + 2.0*i*w*R*Q)]);
% A11 = [MATEXP(1:2,1:2)];
% A12 = [MATEXP(1:2,3:4)];
A21 = [MATEXP(3:4,1:2)];
A22 = [MATEXP(3:4,3:4)];
% Computation of matrix function A(tau)
A = A22\A21;
% Computation of scalar function C(tau)
C = -0.5*beta * trace(log(A22) + t*(M' + 2*R*Q)) + i*w*(r-q);
% Characteristic Function for WMSV model
resu(1,idx) = exp(trace(A*Sigma0) + C + i*w*log(S0));
end
end

• Can you pls add a bit of detail to your answer? What was wrong with the original formulation or code and how did you change it? Oct 10, 2023 at 7:20
• @Alper In a nutshell, you have to evaluate the complex matrix for each argument of the Fourier/Laplace transform u_1,u_1000. This cannot be done with MatLab's element-by-element product so I used a for-loop that would multiply for each argument the matrix by the scalar u_1,....,u_1000 Oct 10, 2023 at 7:28