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