2
$\begingroup$

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?

$\endgroup$

1 Answer 1

1
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Alper
    Oct 10, 2023 at 7:20
  • $\begingroup$ @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 $\endgroup$
    – SimoPape
    Oct 10, 2023 at 7:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.