Home > code > core > cf > cf_wasc_v.m

cf_wasc_v

PURPOSE ^

CF_WASC_V Calculates characteristic function values for a WASC model.

SYNOPSIS ^

function phi = cf_wasc_v(mu, Sigma_0, M, Q, rho, beta, omega, y_t, tau)

DESCRIPTION ^

 CF_WASC_V Calculates characteristic function values for a WASC model.

  phi = cf_wasc_v(mu, Sigma_0, M, Q, rho, beta, omega, y_t, tau)
    calculates the vectorized conditional characteristic function
    of continuous returns of a p-dimensional WASC model as described
    in da Fonseca et al. 2007 (all references to this paper).
    Note that the cf in eq. (21) does not depend on y_t, but it is included
    in the argument list for compatibility reasons (an arbitrary value
    can be passed).
    Assuming that there are T observations (=time points),
    N grid points for the cf evaluation,
    p the dimension of the problem,
    the parameters accepted will be:

    INPUT      mu: a p-dimensional mean vector mu
          Sigma_0: a p x p covariance matrix
                M: a p x p matrix
                Q: an orthogonal p x p matrix
              rho: a p-dimensional correlation vector
             beta: a real value
            omega: a p x N matrix containing the evaluation grid
              y_t: a T x p matrix containing the observations (only
                   the length is used to determin T)
              tau: a real valued time difference

    OUTPUT    phi: a (T-1) x N matrix epresenting the cf value
                   for time (row) and grid point (column).

 See also CF_PCSV_V, CF_GBM_V.

 created by Benedikt Rudolph
 DATE: 20-Aug-2012

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function phi = cf_wasc_v(mu, Sigma_0, M, Q, rho, beta ...
0002                               , omega, y_t, tau)
0003 % CF_WASC_V Calculates characteristic function values for a WASC model.
0004 %
0005 %  phi = cf_wasc_v(mu, Sigma_0, M, Q, rho, beta, omega, y_t, tau)
0006 %    calculates the vectorized conditional characteristic function
0007 %    of continuous returns of a p-dimensional WASC model as described
0008 %    in da Fonseca et al. 2007 (all references to this paper).
0009 %    Note that the cf in eq. (21) does not depend on y_t, but it is included
0010 %    in the argument list for compatibility reasons (an arbitrary value
0011 %    can be passed).
0012 %    Assuming that there are T observations (=time points),
0013 %    N grid points for the cf evaluation,
0014 %    p the dimension of the problem,
0015 %    the parameters accepted will be:
0016 %
0017 %    INPUT      mu: a p-dimensional mean vector mu
0018 %          Sigma_0: a p x p covariance matrix
0019 %                M: a p x p matrix
0020 %                Q: an orthogonal p x p matrix
0021 %              rho: a p-dimensional correlation vector
0022 %             beta: a real value
0023 %            omega: a p x N matrix containing the evaluation grid
0024 %              y_t: a T x p matrix containing the observations (only
0025 %                   the length is used to determin T)
0026 %              tau: a real valued time difference
0027 %
0028 %    OUTPUT    phi: a (T-1) x N matrix epresenting the cf value
0029 %                   for time (row) and grid point (column).
0030 %
0031 % See also CF_PCSV_V, CF_GBM_V.
0032 %
0033 % created by Benedikt Rudolph
0034 % DATE: 20-Aug-2012
0035   
0036   % p: number of assets
0037   % N: number of grid points
0038   [p, N] = size(omega);
0039   T = size(y_t,1);
0040   t = cumsum(repmat(tau,T,1))-tau;
0041   
0042   phi = zeros(T-1, N);% pre-allocate phi
0043   
0044   A_large = compute_a_large(M, Q, rho, omega, tau);
0045   B_large = compute_b_large(M, Q, t(1:end-1)); % note call with t for tau! eq. (21)
0046   trace_M = trace(M);
0047   t_trace_M = trace_M*t(1:end-1);
0048   
0049   i_1_to_p = 1:p;
0050   parfor k=1:N
0051     A_21 = A_large(p+(1:p),2*p*(k-1)+(1:p)); % eq. (14)
0052     A_22 = A_large(p+(1:p),2*p*(k-1)+p+(1:p)); % eq. (14)
0053     A = A_22 \ A_21; % eq. (13)
0054     omega_k = omega(:,k);
0055     c = -0.5*beta*(real(log(det(A_22)))+tau*trace_M ...
0056       + tau*1i*sum(sum((omega_k*rho').*Q'))) ...
0057       + tau*1i*mu'*omega_k;  % eq. (15)
0058     Delta = -1i*A;  % eq. (21)
0059     phi_latest = 1;
0060     for j=1:(T-1)
0061       i_2pj_1 = 2*p*(j-1);
0062       B_11 = B_large(i_2pj_1+i_1_to_p,i_1_to_p);
0063       B_12 = B_large(i_2pj_1+i_1_to_p,p+i_1_to_p);
0064       B_21 = B_large(i_2pj_1+p+i_1_to_p,i_1_to_p);
0065       B_22 = B_large(i_2pj_1+p+i_1_to_p,p+i_1_to_p);
0066       to_inv = 1i*Delta*B_12+B_22;
0067       if ~( rcond(to_inv)<sqrt(eps) )
0068         B = to_inv \ (1i*Delta*B_11+B_21); % eq. (18)
0069         trace_B_Sigma_0 = sum(sum(B.*Sigma_0'));
0070         trace_log_D_B = real(log(det(to_inv)));
0071         % note call with t for tau (!) eq. (21)
0072         C = -0.5*beta*(trace_log_D_B+t_trace_M(j)); % eq. (18)
0073         phi_latest = c+trace_B_Sigma_0+C; % eq. (21)/(17)
0074       end
0075       phi(j,k) = phi_latest;
0076     end
0077   end
0078   phi = exp(phi);
0079 end
0080 
0081 function A_large = compute_a_large(M, Q, rho, omega, tau)
0082   % Calculates the 2p x 2p matrix described in eq. (14).
0083   % That matrix is referred to as "A_large" here in the code.
0084   % The vectorized version of A_large is a 2pT x 2pN matrix.
0085   %
0086   % A_large = compute_a_large(M, Q, rho, omega, tau)
0087   
0088   p = size(omega, 1); % number of assets == dimension of the problem
0089   N = size(omega, 2);
0090   
0091   % pre-calulate matrix elements of the matrix to exponentiated in eq. (14)
0092   
0093   omega_2 = zeros(p, p*N); % pre-allocate
0094   % compute matrix omega_k * omega_k' for each grid point
0095   for k=1:N
0096     omega_2(:,(k-1)*p+1:k*p) = omega(:,k) * omega(:,k)';
0097   end
0098   
0099   omega_flat = reshape(omega, 1, p*N); % flatten omega to a single row vector
0100   
0101   ul = repmat(M, 1, N) + 1i*Q' * rho * omega_flat; % p x p*N matrix
0102   ll = 0.5*(-omega_2 - 1i*repmat(eye(p),1,N).*repmat(omega_flat,p,1)); % p x p*N matrix
0103   ur = -2*(Q'*Q); % p x p matrix
0104   lr = -(repmat(M',1,N) + 1i*kron(omega,rho'*Q)); % p x p*N matrix
0105     
0106   % distribute the four parts of the matrix to be exponentiated
0107   % on a 2*p x 2*p*N matrix
0108   a = zeros(2*p, 2*p*N);
0109   a(repmat(logical([ones(p),zeros(p);zeros(p,2*p)]),1,N)) = ul;
0110   a(repmat(logical([zeros(p),ones(p);zeros(p,2*p)]),1,N))=repmat(ur,1,N);
0111   a(repmat(logical([zeros(p,2*p);ones(p),zeros(p)]),1,N)) = ll;
0112   a(repmat(logical([zeros(p,2*p);zeros(p),ones(p)]),1,N)) = lr;
0113   
0114   % exponentiate
0115   p2 = 2*p;
0116   A_large_cell = cell(1,N);
0117   parfor j=1:N
0118     A_large_cell{j} = expm(tau*a(:,(j-1)*p2+1:j*p2));
0119   end
0120   
0121   A_large = [A_large_cell{:}]; % 2*p x 2*p*N matrix
0122 end
0123 
0124 function B_large = compute_b_large(M, Q, tau)
0125   % Calculates the 2p x 2p matrix described in eq. (18) in a vectorized way.
0126   % That matrix is referred to as "B_large" here in the code.
0127   % The vectorized version of B_large is a 2p(T-1) x 2p matrix.
0128   %
0129   % B_large = compute_b_large(M, Q, tau)
0130   
0131   p = size(M, 1); % number of assets == dimension of the problem
0132   T = length(tau);
0133   
0134   B_log = [M,-2*(Q'*Q);zeros(p),-M']; % matrix of which exp is to be computed
0135   
0136   p2 = 2*p;
0137   B_large = zeros(p2*(T-1), p2);
0138   
0139   for j=1:T
0140     B_large((j-1)*p2+1:j*p2,:) = expm(tau(j)*B_log);
0141   end
0142   
0143 end
0144

Generated on Mon 29-Apr-2013 19:29:13 by m2html © 2005