| gradient(A,B,C,V,W,P,R,S1,S2,Q1,Q2,met,met_sylv)
|
function [Pdif,Rdif,S1dif,S2dif,Q1dif,Q2dif,INFO] = gradient(A,B,C,V,W,P,R,S1,S2,Q1,Q2,met,met_sylv)
% PURPOSE: finds first differential of matrices representing solution with
% respect to parameter.
% Model takes the form
%
% 0 = A(p) y_t + B(p) y_t+1 + C(p) E_t{y_t+1} + V(p) epsilon_t + W(p) epsilon_t+1
%
% where $p$ is a scalar, and its solution takes the form
%
% y_t = R(p) u_t + S_1(p) epsilon_t + S_2(p) omega_t
% u_t = P(p) u_t-1 + Q_2(p) epsilon_t + Q_2(p) omega_t
%
% where u_t is a state variable, omega_t is i.i.d. random
% variable with mean zero. This routine calculates R'(p), P'(p),
% S_1'(p), S_2'(p), Q_1'(p), Q_2'(p) around p = 0.
%
% USAGE: [Pdif,Rdif,S1dif,S2dif,Q1dif,Q2dif,INFO] = gradient(A,B,C,V,W,P,R,S1,S2,Q1,Q2,met,met_sylv)
% where:
% A,B,C,V,W cells 1xk. The first element is value of given
% matrix for p = 0. The i+1-th element represents
% a differential of a given matrix with respect
% to i-th parameter
% P,R,S1,S2,Q1,Q2 value of matrices forming solution to
% the model at p=0
% met method of calculating null spaces
% (optional)
% 1 - qr with pivoting (default)
% 2 - svd
% met_sylv method of calculating generalized
% Sylvester equation (optional)
% 1 - using gsylvester (default)
% 2 - vectorization
%
% Pdif,Rdif,S1dif,S2dif,Q1dif,Q2dif
% cells 1xk. The i-th element represents
% a differential of a given matrix with respect
% to i-th parameter
% INFO cell 1xk. The i-th element concerns
% i-th parameter
% 0: one of the differentian does not
% exist
% 1: success
% 2: solving generalized Sylvester
% equation failed
%
% COMMENTS:
%
% Copyright (c) Pawel Kowal (2006)
% All rights reserved
% LREM_SOLVE toolbox is available free for noncommercial academic use only.
% pkowal3@sgh.waw.pl
if nargin<12
met = 1;
end
if nargin<13
met_sylv = 1;
end
Pdif = {};
Rdif = {};
S1dif = {};
S2dif = {};
Q1dif = {};
Q2dif = {};
info = 0;
%---------------------------- PREPARING MATRICES --------------------------
tol = 300*100*eps;
%------------- DETERMINISTIC --------------------
[BC_bar,U_BC,k] = putv(B{1}+C{1},met,tol);
[BCR_bar,U_BCR,k2] = putv(BC_bar*R,met,tol);
M_BC = U_BC(:,1:k);
N_BC = U_BC(:,k+1:end);
M_BCR = U_BCR(:,1:k2);
N_BCR = U_BCR(:,k2+1:end);
[R_bar,U_R,k] = putv(R,met,tol);
M_R = U_R(:,1:k);
N_R = U_R(:,k+1:end);
H = null2(N_BC'*A{1}*N_R,met,tol);
Y_B1 = M_BCR'*M_BC'*A{1}*N_R*H;
Y_B2 = M_BCR'*BC_bar*N_R*H;
Y_A = BCR_bar;
H2 = N_BCR'*M_BC'*A{1}*N_R*H;
H3 = N_BCR'*BC_bar*N_R*H;
n = size(P,1);
if met_sylv==2
VEC = kron(eye(n),H2)+kron(P',H3);
end
G_base = -(N_BC'*A{1}*N_R)\(N_BC');
G_util1 = A{1}*N_R*G_base;
G_util2 = BC_bar*N_R*G_base;
%------------- STOCHASTIC --------------------
[A_bar,U_A,k] = putv(A{1},met,tol);
M_A = U_A(:,1:k);
N_A = U_A(:,k+1:end);
PSI = null2(A{1},met,tol);
[B_bar,U_B,k] = putv(B{1}*[R PSI],met,tol);
M_B = U_B(:,1:k);
N_B = U_B(:,k+1:end);
%---------------------------- LOOP OVER PARAMETERS ------------------------
n = length(A)-1;
INFO = {};
for i=1:1:n
Rdif{i} = 0*R;
Pdif{i} = 0*P;
Q1dif{i} = 0*Q1;
Q2dif{i} = 0*Q2;
S1dif{i} = 0*S1;
S2dif{i} = 0*S2;
%------------- DETERMINISTIC PART --------------
G1 = A{i+1}*R+(B{i+1}+C{i+1})*R*P;
G = G_base*G1;
G2 = M_BC'*(G1+G_util1*G1)+G_util2*G1*P;
G22 = N_BCR'*G2;
if met_sylv==1
[Z,dum,info] = gsylvester(H2,P,-G22,H3,eye(size(P)),0,2);
else
mz = size(R,2);
nz = size(H,2);
Z = reshape(-VEC\G22(:),nz,mz);
end
RR = N_R*G+N_R*H*Z;
Rdif{i} = RR;
Y_sum = M_BCR'*G2+Y_B1*Z+Y_B2*Z*P;
Y = -Y_A\Y_sum;
Pdif{i} = Y;
%------------- STOCHASTIC PART ------------------
STOCH_1 = A{i+1}*S1+V{i+1};
STOCH_2 = A{i+1}*S2;
if norm(N_A'*STOCH_1,1)>tol
INFO{i} = 0;
continue;
end
if norm(N_A'*STOCH_2,1)>tol
INFO{i} = 0;
continue;
end
SS1 = -A_bar\(M_A'*STOCH_1);
SS2 = -A_bar\(M_A'*STOCH_2);
STOCH_3 = B{i+1}*(R*Q1+S1)+W{i+1}+B{1}*RR*Q1+B{1}*SS1;
STOCH_4 = B{i+1}*(R*Q2+S2) +B{1}*RR*Q2+B{1}*SS2;
if norm(N_B'*STOCH_3,1)>tol
INFO{i} = 0;
continue;
end
if norm(N_B'*STOCH_4,1)>tol
INFO{i} = 0;
continue;
end
PHI_1 = -B_bar\(M_B'*STOCH_3);
PHI_2 = -B_bar\(M_B'*STOCH_4);
Q1dif{i} = PHI_1(1:size(R,2),:);
Q2dif{i} = PHI_2(1:size(R,2),:);
Z1 = PHI_1(size(R,2)+1:end,:);
Z2 = PHI_2(size(R,2)+1:end,:);
S1dif{i} = SS1+PSI*Z1;
S2dif{i} = SS2+PSI*Z2;
INFO{i} = 1;
if info~=0
INFO{i} = 2;
end
end |
|
