Commit 717095ae authored by MichelJuillard's avatar MichelJuillard

fixed bug at order 2, when a variable is absent at the current period;

cleaned code that is useless since we transform leads and lags on
period > 1 (hand cherry-picked from 0303b1c0)
parent b4189182
......@@ -49,6 +49,8 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
lead_lag_incidence = M_.lead_lag_incidence;
info = 0;
if M_.maximum_endo_lag == 0 && options_.order > 1
......@@ -64,7 +66,7 @@ end
xlen = M_.maximum_endo_lead + M_.maximum_endo_lag + 1;
klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1;
iyv = M_.lead_lag_incidence';
iyv = lead_lag_incidence';
iyv = iyv(:);
iyr0 = find(iyv) ;
it_ = M_.maximum_lag + 1 ;
......@@ -207,11 +209,11 @@ npred = dr.npred;
nboth = dr.nboth;
order_var = dr.order_var;
nd = size(kstate,1);
nz = nnz(M_.lead_lag_incidence);
nz = nnz(lead_lag_incidence);
sdyn = M_.endo_nbr - nstatic;
[junk,cols_b,cols_j] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, ...
[junk,cols_b,cols_j] = find(lead_lag_incidence(M_.maximum_endo_lag+1, ...
order_var));
b = zeros(M_.endo_nbr,M_.endo_nbr);
b(:,cols_b) = jacobia_(:,cols_j);
......@@ -256,7 +258,7 @@ if M_.maximum_endo_lead == 0
info(2) = temp'*temp;
end
if options_.loglinear == 1
klags = find(M_.lead_lag_incidence(1,:));
klags = find(lead_lag_incidence(1,:));
dr.ghx = repmat(1./dr.ys,1,size(dr.ghx,2)).*dr.ghx.* ...
repmat(dr.ys(klags),size(dr.ghx,1),1);
dr.ghu = repmat(1./dr.ys,1,size(dr.ghu,2)).*dr.ghu;
......@@ -310,7 +312,7 @@ if (options_.aim_solver == 1) && (task == 0)
error('Problem with AIM solver - Try to remove the "aim_solver" option')
end
else % use original Dynare solver
k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
k1 = lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
a = aa(:,nonzeros(k1'));
b(:,cols_b) = aa(:,cols_j);
b10 = b(1:nstatic,1:nstatic);
......@@ -473,8 +475,8 @@ end
%exogenous deterministic variables
if M_.exo_det_nbr > 0
f1 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+2:end,order_var))));
f0 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var))));
f1 = sparse(jacobia_(:,nonzeros(lead_lag_incidence(M_.maximum_endo_lag+2:end,order_var))));
f0 = sparse(jacobia_(:,nonzeros(lead_lag_incidence(M_.maximum_endo_lag+1,order_var))));
fudet = sparse(jacobia_(:,nz+M_.exo_nbr+1:end));
M1 = inv(f0+[zeros(M_.endo_nbr,nstatic) f1*gx zeros(M_.endo_nbr,nyf-nboth)]);
M2 = M1*f1;
......@@ -490,56 +492,26 @@ if options_.order == 1
end
% Second order
%tempex = oo_.exo_simul ;
%hessian = real(hessext('ff1_',[z; oo_.exo_steady_state]))' ;
kk = flipud(cumsum(flipud(M_.lead_lag_incidence(M_.maximum_endo_lag+1:end,order_var)),1));
if M_.maximum_endo_lag > 0
kk = [cumsum(M_.lead_lag_incidence(1:M_.maximum_endo_lag,order_var),1); kk];
end
kk = kk';
kk = find(kk(:));
nk = size(kk,1) + M_.exo_nbr + M_.exo_det_nbr;
k1 = M_.lead_lag_incidence(:,order_var);
k1 = k1';
k1 = k1(:);
k1 = k1(kk);
k2 = find(k1);
kk1(k1(k2)) = k2;
kk1 = [kk1 length(k1)+1:length(k1)+M_.exo_nbr+M_.exo_det_nbr];
kk = reshape([1:nk^2],nk,nk);
kk1 = kk(kk1,kk1);
%[junk,junk,hessian] = feval([M_.fname '_dynamic'],z, oo_.exo_steady_state);
hessian(:,kk1(:)) = hessian1;
k1 = nonzeros(lead_lag_incidence(:,order_var)');
kk = [k1; length(k1)+(1:M_.exo_nbr+M_.exo_det_nbr)'];
nk = size(kk,1);
kk1 = reshape([1:nk^2],nk,nk);
kk1 = kk1(kk,kk);
hessian = hessian1(:,kk1(:));
clear hessian1
%oo_.exo_simul = tempex ;
%clear tempex
n1 = 0;
n2 = np;
zx = zeros(np,np);
zu=zeros(np,M_.exo_nbr);
for i=2:M_.maximum_endo_lag+1
k1 = sum(kstate(:,2) == i);
zx(n1+1:n1+k1,n2-k1+1:n2)=eye(k1);
n1 = n1+k1;
n2 = n2-k1;
end
kk = flipud(cumsum(flipud(M_.lead_lag_incidence(M_.maximum_endo_lag+1:end,order_var)),1));
zx(1:np,:)=eye(np);
k0 = [1:M_.endo_nbr];
gx1 = dr.ghx;
hu = dr.ghu(nstatic+[1:npred],:);
zx = [zx; gx1];
zu = [zu; dr.ghu];
for i=1:M_.maximum_endo_lead
k1 = find(kk(i+1,k0) > 0);
zu = [zu; gx1(k1,1:npred)*hu];
gx1 = gx1(k1,:)*hx;
zx = [zx; gx1];
kk = kk(:,k0);
k0 = k1;
end
k0 = find(lead_lag_incidence(M_.maximum_endo_lag+1,order_var)');
zx = [zx; gx1(k0,:)];
zu = [zu; dr.ghu(k0,:)];
k1 = find(lead_lag_incidence(M_.maximum_endo_lag+2,order_var)');
zu = [zu; gx1(k1,:)*hu];
zx = [zx; gx1(k1,:)*hx];
zx=[zx; zeros(M_.exo_nbr,np);zeros(M_.exo_det_nbr,np)];
zu=[zu; eye(M_.exo_nbr);zeros(M_.exo_det_nbr,M_.exo_nbr)];
[nrzx,nczx] = size(zx);
......@@ -550,46 +522,19 @@ rhs = -rhs;
%lhs
n = M_.endo_nbr+sum(kstate(:,2) > M_.maximum_endo_lag+1 & kstate(:,2) < M_.maximum_endo_lag+M_.maximum_endo_lead+1);
A = zeros(n,n);
B = zeros(n,n);
A(1:M_.endo_nbr,1:M_.endo_nbr) = jacobia_(:,M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var));
A = zeros(M_.endo_nbr,M_.endo_nbr);
B = zeros(M_.endo_nbr,M_.endo_nbr);
A(:,k0) = jacobia_(:,nonzeros(lead_lag_incidence(M_.maximum_endo_lag+1,order_var)));
% variables with the highest lead
k1 = find(kstate(:,2) == M_.maximum_endo_lag+M_.maximum_endo_lead+1);
if M_.maximum_endo_lead > 1
k2 = find(kstate(:,2) == M_.maximum_endo_lag+M_.maximum_endo_lead);
[junk,junk,k3] = intersect(kstate(k1,1),kstate(k2,1));
else
k2 = [1:M_.endo_nbr];
k3 = kstate(k1,1);
end
k1 = find(kstate(:,2) == M_.maximum_endo_lag+2);
% Jacobian with respect to the variables with the highest lead
B(1:M_.endo_nbr,end-length(k2)+k3) = jacobia_(:,kstate(k1,3)+M_.endo_nbr);
fyp = jacobia_(:,kstate(k1,3)+M_.endo_nbr);
B(:,nstatic+npred-dr.nboth+1:end) = fyp;
offset = M_.endo_nbr;
k0 = [1:M_.endo_nbr];
gx1 = dr.ghx;
for i=1:M_.maximum_endo_lead-1
k1 = find(kstate(:,2) == M_.maximum_endo_lag+i+1);
[k2,junk,k3] = find(kstate(k1,3));
A(1:M_.endo_nbr,offset+k2) = jacobia_(:,k3+M_.endo_nbr);
n1 = length(k1);
A(offset+[1:n1],nstatic+[1:npred]) = -gx1(kstate(k1,1),1:npred);
gx1 = gx1*hx;
A(offset+[1:n1],offset+[1:n1]) = eye(n1);
n0 = length(k0);
E = eye(n0);
if i == 1
[junk,junk,k4]=intersect(kstate(k1,1),[1:M_.endo_nbr]);
else
[junk,junk,k4]=intersect(kstate(k1,1),kstate(k0,1));
end
i1 = offset-n0+n1;
B(offset+[1:n1],offset-n0+[1:n0]) = -E(k4,:);
k0 = k1;
offset = offset + n1;
end
[junk,k1,k2] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+M_.maximum_endo_lead+1,order_var));
A(1:M_.endo_nbr,nstatic+1:nstatic+npred)=...
A(1:M_.endo_nbr,nstatic+[1:npred])+jacobia_(:,k2)*gx1(k1,1:npred);
A(1:M_.endo_nbr,nstatic+[1:npred])+fyp*gx1(k1,1:npred);
C = hx;
D = [rhs; zeros(n-M_.endo_nbr,size(rhs,2))];
......@@ -600,16 +545,10 @@ mexErrCheck('gensylv', err);
%ghxu
%rhs
hu = dr.ghu(nstatic+1:nstatic+npred,:);
%kk = reshape([1:np*np],np,np);
%kk = kk(1:npred,1:npred);
%rhs = -hessian*kron(zx,zu)-f1*dr.ghxx(end-nyf+1:end,kk(:))*kron(hx(1:npred,:),hu(1:npred,:));
[err, rhs] = sparse_hessian_times_B_kronecker_C(hessian,zx,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian,zx,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
nyf1 = sum(kstate(:,2) == M_.maximum_endo_lag+2);
hu1 = [hu;zeros(np-npred,M_.exo_nbr)];
%B1 = [B(1:M_.endo_nbr,:);zeros(size(A,1)-M_.endo_nbr,size(B,2))];
[nrhx,nchx] = size(hx);
[nrhu1,nchu1] = size(hu1);
......@@ -624,10 +563,7 @@ dr.ghxu = A\rhs;
%ghuu
%rhs
kk = reshape([1:np*np],np,np);
kk = kk(1:npred,1:npred);
[err, rhs] = sparse_hessian_times_B_kronecker_C(hessian,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
[err, B1] = A_times_B_kronecker_C(B*dr.ghxx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
......@@ -647,7 +583,8 @@ dr.ghuu = dr.ghuu(1:M_.endo_nbr,:);
% reordering predetermined variables in diminishing lag order
O1 = zeros(M_.endo_nbr,nstatic);
O2 = zeros(M_.endo_nbr,M_.endo_nbr-nstatic-npred);
LHS = jacobia_(:,M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var));
LHS = zeros(M_.endo_nbr,M_.endo_nbr);
LHS(:,k0) = jacobia_(:,nonzeros(lead_lag_incidence(M_.maximum_endo_lag+1,order_var)));
RHS = zeros(M_.endo_nbr,M_.exo_nbr^2);
kk = find(kstate(:,2) == M_.maximum_endo_lag+2);
gu = dr.ghu;
......@@ -655,51 +592,19 @@ guu = dr.ghuu;
Gu = [dr.ghu(nstatic+[1:npred],:); zeros(np-npred,M_.exo_nbr)];
Guu = [dr.ghuu(nstatic+[1:npred],:); zeros(np-npred,M_.exo_nbr*M_.exo_nbr)];
E = eye(M_.endo_nbr);
M_.lead_lag_incidenceordered = flipud(cumsum(flipud(M_.lead_lag_incidence(M_.maximum_endo_lag+1:end,order_var)),1));
if M_.maximum_endo_lag > 0
M_.lead_lag_incidenceordered = [cumsum(M_.lead_lag_incidence(1:M_.maximum_endo_lag,order_var),1); M_.lead_lag_incidenceordered];
end
M_.lead_lag_incidenceordered = M_.lead_lag_incidenceordered';
M_.lead_lag_incidenceordered = M_.lead_lag_incidenceordered(:);
k1 = find(M_.lead_lag_incidenceordered);
M_.lead_lag_incidenceordered(k1) = [1:length(k1)]';
M_.lead_lag_incidenceordered =reshape(M_.lead_lag_incidenceordered,M_.endo_nbr,M_.maximum_endo_lag+M_.maximum_endo_lead+1)';
kh = reshape([1:nk^2],nk,nk);
kp = sum(kstate(:,2) <= M_.maximum_endo_lag+1);
E1 = [eye(npred); zeros(kp-npred,npred)];
H = E1;
hxx = dr.ghxx(nstatic+[1:npred],:);
for i=1:M_.maximum_endo_lead
for j=i:M_.maximum_endo_lead
[junk,k2a,k2] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+j+1,order_var));
[junk,k3a,k3] = ...
find(M_.lead_lag_incidenceordered(M_.maximum_endo_lag+j+1,:));
nk3a = length(k3a);
[err, B1] = sparse_hessian_times_B_kronecker_C(hessian(:,kh(k3,k3)),gu(k3a,:),options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
RHS = RHS + jacobia_(:,k2)*guu(k2a,:)+B1;
end
% LHS
[junk,k2a,k2] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+i+1,order_var));
LHS = LHS + jacobia_(:,k2)*(E(k2a,:)+[O1(k2a,:) dr.ghx(k2a,:)*H O2(k2a,:)]);
if i == M_.maximum_endo_lead
break
end
[junk,k2a,k2] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+2,order_var));
[B1, err] = sparse_hessian_times_B_kronecker_C(hessian(:,kh(k2,k2)),gu(k2a,:),options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
RHS = RHS + jacobia_(:,k2)*guu(k2a,:)+B1;
% LHS
LHS = LHS + jacobia_(:,k2)*(E(k2a,:)+[O1(k2a,:) dr.ghx(k2a,:)*H O2(k2a,:)]);
kk = find(kstate(:,2) == M_.maximum_endo_lag+i+1);
gu = dr.ghx*Gu;
[nrGu,ncGu] = size(Gu);
[err, G1] = A_times_B_kronecker_C(dr.ghxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
mexErrCheck('A_times_B_kronecker_C', err);
[err, G2] = A_times_B_kronecker_C(hxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
mexErrCheck('A_times_B_kronecker_C', err);
guu = dr.ghx*Guu+G1;
Gu = hx*Gu;
Guu = hx*Guu;
Guu(end-npred+1:end,:) = Guu(end-npred+1:end,:) + G2;
H = E1 + hx*H;
end
RHS = RHS*M_.Sigma_e(:);
dr.fuu = RHS;
%RHS = -RHS-dr.fbias;
......
......@@ -38,9 +38,9 @@ endo_nbr = M_.endo_nbr;
lead_lag_incidence = M_.lead_lag_incidence;
klen = max_lag + max_lead + 1;
fwrd_var = find(any(lead_lag_incidence(max_lag+2:end,:),1))';
fwrd_var = find(lead_lag_incidence(max_lag+2:end,:))';
if max_lag > 0
pred_var = find(any(lead_lag_incidence(1,:),1))';
pred_var = find(lead_lag_incidence(1,:))';
both_var = intersect(pred_var,fwrd_var);
pred_var = setdiff(pred_var,both_var);
fwrd_var = setdiff(fwrd_var,both_var);
......@@ -61,11 +61,11 @@ inv_order_var(order_var) = (1:endo_nbr);
if max_lag > 0
kmask = [];
if max_lead > 0
kmask = [cumsum(flipud(lead_lag_incidence(max_lag+2:end,order_var)),1)] ;
kmask = lead_lag_incidence(max_lag+2,order_var) ;
end
kmask = [kmask; flipud(cumsum(lead_lag_incidence(1,order_var),1))] ;
kmask = [kmask; lead_lag_incidence(1,order_var)] ;
else
kmask = cumsum(flipud(lead_lag_incidence(max_lag+2:klen,order_var)),1) ;
kmask = lead_lag_incidence(max_lag+2,order_var) ;
end
kmask = kmask';
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment