stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

optimal policy: bug correction in Lagrange multipliers elimination (mult_elimination.m)

time-shift
Michel Juillard 2010-05-08 10:10:19 +02:00
parent b903776e0c
commit f12c64cc84
1 changed files with 130 additions and 120 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

250
matlab/mult_elimination.m
View File

@ -1,120 +1,130 @@
function dr=mult_elimination(varlist,M_, options_, oo_) function dr=mult_elimination(varlist,M_, options_, oo_)
% function mult_elimination() % function mult_elimination()
% replaces Lagrange multipliers in Ramsey policy by lagged value of state % replaces Lagrange multipliers in Ramsey policy by lagged value of state
% and shock variables % and shock variables
% %
% INPUT % INPUT
% none % none
% %
% OUTPUT % OUTPUT
% dr: a structure with the new decision rule % dr: a structure with the new decision rule
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 2003-2008 Dynare Team % Copyright (C) 2003-2008 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
% Dynare is free software: you can redistribute it and/or modify % Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by % it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or % the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version. % (at your option) any later version.
% %
% Dynare is distributed in the hope that it will be useful, % Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of % but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details. % GNU General Public License for more details.
% %
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
dr = oo_.dr; dr = oo_.dr;
nstatic = dr.nstatic; nstatic = dr.nstatic;
npred = dr.npred; npred = dr.npred;
order_var = dr.order_var; order_var = dr.order_var;
nstates = M_.endo_names(order_var(nstatic+(1:npred)),:); nstates = M_.endo_names(order_var(nstatic+(1:npred)),:);
il = strmatch('mult_',nstates); il = strmatch('mult_',nstates);
nil = setdiff(1:dr.npred,il); nil = setdiff(1:dr.npred,il);
m_nbr = length(il); m_nbr = length(il);
nm_nbr = length(nil); nm_nbr = length(nil);
AA1 = dr.ghx(:,nil); AA1 = dr.ghx(:,nil);
AA2 = dr.ghx(:,il); AA2 = dr.ghx(:,il);
A1 = dr.ghx(nstatic+(1:npred),nil); A1 = dr.ghx(nstatic+(1:npred),nil);
A2 = dr.ghx(nstatic+(1:npred),il); A2 = dr.ghx(nstatic+(1:npred),il);
B = dr.ghu(nstatic+(1:npred),:); B = dr.ghu(nstatic+(1:npred),:);
A11 = A1(nil,:); A11 = A1(nil,:);
A21 = A1(il,:); A21 = A1(il,:);
A12 = A2(nil,:); A12 = A2(nil,:);
A22 = A2(il,:); A22 = A2(il,:);
B1 = B(nil,:);
[Q1,R1,E1] = qr(A2); B2 = B(il,:);
n1 = sum(abs(diag(R1)) > 1e-8);
[Q1,R1,E1] = qr([A12; A22]);
Q1_12 = Q1(1:nm_nbr,n1+1:end); n1 = sum(abs(diag(R1)) > 1e-8);
Q1_22 = Q1(nm_nbr+1:end,n1+1:end);
[Q2,R2,E2] = qr(Q1_22'); Q1_12 = Q1(1:nm_nbr,n1+1:end);
n2 = sum(abs(diag(R2)) > 1e-8); Q1_22 = Q1(nm_nbr+(1:m_nbr),n1+1:end);
[Q2,R2,E2] = qr(Q1_22');
R2_1 = inv(R2(1:n2,1:n2)); n2 = sum(abs(diag(R2)) > 1e-8);
M1(order_var,:) = AA1 - AA2*E2*[R2_1*Q2(:,1:n2)'*Q1_12'; zeros(m_nbr-n2,nm_nbr)]; R2_1 = inv(R2(1:n2,1:n2));
M2(order_var,:) = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*A1; zeros(m_nbr-n2,length(nil))];
M3(order_var,:) = dr.ghu; M1 = AA1 - AA2*E2*[R2_1*Q2(:,1:n2)'*Q1_12'; zeros(m_nbr-n2,nm_nbr)];
M4(order_var,:) = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*B; zeros(m_nbr-n2,size(B,2))]; M2 = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*[A11;A21]; zeros(m_nbr-n2,length(nil))];
M3 = dr.ghu;
endo_nbr = M_.orig_model.endo_nbr; M4 = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*[B1;B2]; zeros(m_nbr-n2,size(B,2))];
exo_nbr = M_.exo_nbr;
k1 = nstatic+(1:npred);
lead_lag_incidence = M_.lead_lag_incidence(:,1:endo_nbr+exo_nbr); k1 = k1(nil);
lead_lag_incidence1 = lead_lag_incidence(1,:) > 0;
maximum_lag = M_.maximum_lag; endo_nbr = M_.orig_model.endo_nbr;
for i=1:maximum_lag-1 exo_nbr = M_.exo_nbr;
lead_lag_incidence1 = [lead_lag_incidence1; lead_lag_incidence(i,:)| ...
lead_lag_incidence(i+1,:)]; lead_lag_incidence = M_.lead_lag_incidence(:,1:endo_nbr+exo_nbr);
end lead_lag_incidence1 = lead_lag_incidence(1,:) > 0;
lead_lag_incidence1 = [lead_lag_incidence1; ... maximum_lag = M_.maximum_lag;
lead_lag_incidence(M_.maximum_lag,:) > 0]; for i=1:maximum_lag-1
k = find(lead_lag_incidence1'); lead_lag_incidence1 = [lead_lag_incidence1; lead_lag_incidence(i,:)| ...
lead_lag_incidence1 = zeros(size(lead_lag_incidence1')); lead_lag_incidence(i+1,:)];
lead_lag_incidence1(k) = 1:length(k); end
lead_lag_incidence1 = lead_lag_incidence1'; lead_lag_incidence1 = [lead_lag_incidence1; ...
lead_lag_incidence(M_.maximum_lag,:) > 0];
kstate = zeros(0,2); k = find(lead_lag_incidence1');
for i=maximum_lag:-1:1 lead_lag_incidence1 = zeros(size(lead_lag_incidence1'));
k = find(lead_lag_incidence(i,:)); lead_lag_incidence1(k) = 1:length(k);
kstate = [kstate; [k' repmat(i+1,length(k),1)]]; lead_lag_incidence1 = lead_lag_incidence1';
end
kstate = zeros(0,2);
dr.M1 = M1; for i=maximum_lag:-1:1
dr.M2 = M2; k = find(lead_lag_incidence(i,:));
dr.M3 = M3; kstate = [kstate; [k' repmat(i+1,length(k),1)]];
dr.M4 = M4; end
nvar = length(varlist); dr.M1 = M1;
nspred = dr.nspred; dr.M2 = M2;
nspred = 6; dr.M3 = M3;
if nvar > 0 dr.M4 = M4;
res_table = zeros(2*(nspred+M_.exo_nbr),nvar);
headers = 'Variables'; nvar = length(varlist);
for i=1:length(varlist) nspred = dr.nspred;
k = strmatch(varlist{i},M_.endo_names(dr.order_var,:),'exact'); nspred = 6;
headers = strvcat(headers,varlist{i}); if nvar > 0
res_table = zeros(2*(nspred+M_.exo_nbr),nvar);
res_table(1:nspred,i) = M1(k,:)'; headers = 'Variables';
res_table(nspred+(1:nspred),i) = M2(k,:)'; for i=1:length(varlist)
res_table(2*nspred+(1:M_.exo_nbr),i) = M3(k,:)'; k = strmatch(varlist{i},M_.endo_names(dr.order_var,:),'exact');
res_table(2*nspred+M_.exo_nbr+(1:M_.exo_nbr),i) = M4(k,:)'; headers = strvcat(headers,varlist{i});
end
res_table(1:nspred,i) = M1(k,:)';
my_title='ELIMINATION OF THE MULTIPLIERS'; res_table(nspred+(1:nspred),i) = M2(k,:)';
lab1 = M_.endo_names(dr.order_var(dr.nstatic+[ 1 2 5:8]),:); res_table(2*nspred+(1:M_.exo_nbr),i) = M3(k,:)';
labels = strvcat(lab1,lab1,M_.exo_names,M_.exo_names); res_table(2*nspred+M_.exo_nbr+(1:M_.exo_nbr),i) = M4(k,:)';
lh = size(labels,2)+2; end
dyntable(my_title,headers,labels,res_table,lh,10,6);
disp(' ') my_title='ELIMINATION OF THE MULTIPLIERS';
end lab1 = M_.endo_names(dr.order_var(dr.nstatic+[ 1 2 5:8]),:);
labels = strvcat(lab1,lab1,M_.exo_names,M_.exo_names);
lh = size(labels,2)+2;
dyntable(my_title,headers,labels,res_table,lh,10,6);
disp(' ')
end
function x=varsol(a,b,c)
x = (kron(a,a)-kron(b,b))\c(:);
x = reshape(x,size(a,1),size(a,1));