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?
