Evaluating the Wishart Heston Characteristic Function Across Fourier Inputs

# Characteristic Function for Wishart Heston Model


# 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)

\begin{equation} A(\tau) = F(\tau)^{-1}G(\tau) = A_{22}(\tau)^{-1}A_{21}(\tau) \end{equation} where $F(\tau), G(\tau)$ solve a system of $(2n)$ linear ODE whose solution is given by \begin{equation} \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] \end{equation} instead the solution for $C(\tau)$ is given by \begin{equation} C(\tau) = -\frac{\beta}{2} Tr\Big[ \log F(\tau) + \tau(M^T + 2i\gamma RQ) \Big] + i\gamma(r-q)\tau \end{equation}

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?

## Answer by SimoPape (score 1, accepted)

https://quant.stackexchange.com/a/76940

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
```