Función característica para el modelo Wishart Heston

Question

No sé si este es el lugar correcto (como mucho cerrarán el post). De todos modos, estoy tratando de implementar la función característica del modelo de volatilidad estocástica Heston Wishart ilustrado en el trabajo de Da Fonseca, Grasselli y Tebaldi [2008] para usar métodos de transformada de Fourier.

La función característica es

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

donde $Y_{t} := \log S_{t}$ y $\Sigma_{t}$ sigue el SDE de wishart.

Utilizando el procedimiento de linearización (como se explica en el papel), podríamos resolver la ODE de Riccati de la matriz $A(\tau)$. La solución para $A(\tau)$ está dada por la ecuación (12)

\begin{equation} A(\tau) = F(\tau)^{-1}G(\tau) = A_{22}(\tau)^{-1}A_{21}(\tau) \end{equation} donde $F(\tau), G(\tau)$ resuelven un sistema de ODE lineales de $(2n)$ cuya solución está dada por \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} en lugar de la solución para $C(\tau)$ que es \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}

donde $F(\tau) = A_{22}(\tau)$. Observar que Da Fonsenca, Grasselli y Tebaldi derivaron la función generadora de momentos del logaritmo del precio. Por lo tanto, para obtener la función característica debemos evaluar $\gamma := i\gamma$ donde $\gamma \in \mathbb{R}$.

Ahora, intenté implementar la función característica del modelo WMSV en MatLab pero tengo algunos problemas. Publicaré un esbozo de código

clear variables; close all; clc

% -----------------------------
% Parámetros
% -----------------------------
S0 = 100; % Precio inicial de la acción
r = 0.05; % Tasa de interés libre de riesgo
q = 0; % Rendimiento de dividendos
t = 1.0; % Tiempo hasta el vencimiento
beta = 3; % Rol de la condición de Feller

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];  % Varianza inicial 

% u = linspace(-10,10,1000);
u = 3.5; % por ejemplo

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

% ----------------------------------------------------------
% Aquí se calcula la función característica del modelo WMSV
% ----------------------------------------------------------
function resu = cf_wmsv(u,t,r,q,M,R,Q,Sigma0,S0,beta)
    i = complex(0,1);
    % Exponencial de la matriz (14) en 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)];
    % Cálculo de la función de matriz A(tau) ecuación (12-15) en Da Fonseca et al.
    A = A22\A21;
    % Cálculo de la función escalar C(tau) 
    C = -0.5*beta*trace(log(A22) + t*(M' + 2*R*Q)) + i*u*(r-q);
    % Función Característica para el modelo WMSV ecuación (5) en Da Fonseca et al.
    resu = exp(trace(A*Sigma0) + C + i*u*log(S0));
end

En el código el "rol" de $\gamma$ está siendo reemplazado por $u$. La salida debería ser un vector doble complejo $1\times \mathtt{length(u)}$. Si considero una entrada escalar $u$ no tengo problema, pero si considero (Aun si reemplazo el producto matricial $*$ con el producto elemento por elemento $.*$) no obtengo el resultado deseado.

Consideré que podría usar un bucle for y evaluar la matriz compleja 4*4 elemento por elemento, pero no sé cómo ajustar el resto. ¿Alguien puede darme una sugerencia o solución?

Answer 1

He resuelto el problema por mí mismo. En resumen, uno tiene que evaluar la matriz compleja para cada argumento de la transformada de Fourier/Laplace $u_1, \dots, u_{1000}$. Esto no se puede hacer con el producto elemento por elemento de MatLab. Así que utilicé un bucle for que multiplicaría por cada argumento la matriz por el escalar $u_1, \dots, u_{1000}$.

Aquí está el código revisado:

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)
        % Exponencial de la matriz (14) en 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)];
        % Cálculo de la función de la matriz A(tau)
        A = A22\A21;
        % Cálculo de la función escalar C(tau)
        C = -0.5*beta * trace(log(A22) + t*(M' + 2*R*Q)) + i*w*(r-q);
        % Función Característica para el modelo WMSV
        resu(1,idx) = exp(trace(A*Sigma0) + C + i*w*log(S0));
    end
end