Skip to content
GitLab
Commit 4d51f38b authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Added the possibility to use the cyclic reduction algorithm without the block option.

parent b02c83a9
Hide whitespace changes
Inline Side-by-side
matlab/disp_dr.m
View file @ 4d51f38b
function disp_dr(dr,order,var_list)
% Display the decision rules
%
......@@ -102,11 +103,7 @@ end
%
for k=1:nx
flag = 0;
if options_.block
str1 = subst_auxvar(dr.state_var(k),-1);
else
str1 = subst_auxvar(k1(klag(k,1)),klag(k,2)-M_.maximum_lag-2);
end;
str1 = subst_auxvar(dr.state_var(k),-1);
str = sprintf('%-20s',str1);
for i=1:nvar
x = dr.ghx(ivar(i),k);
......
matlab/dyn_first_order_solver.m
View file @ 4d51f38b
function [dr,info] = dyn_first_order_solver(jacobia,M_,dr,options,task)
function [dr,info] = dyn_first_order_solver(jacobia,DynareModel,dr,DynareOptions,task)
%@info:
%! @deftypefn {Function File} {[@var{dr},@var{info}] =} dyn_first_order_solver (@var{jacobia},@var{M_},@var{dr},@var{options},@var{task})
%! @deftypefn {Function File} {[@var{dr},@var{info}] =} dyn_first_order_solver (@var{jacobia},@var{DynareModel},@var{dr},@var{DynareOptions},@var{task})
%! @anchor{dyn_first_order_solver}
%! @sp 1
%! Computes the first order reduced form of the DSGE model
......@@ -11,7 +11,7 @@ function [dr,info] = dyn_first_order_solver(jacobia,M_,dr,options,task)
%! @table @ @var
%! @item jacobia
%! Matrix containing the Jacobian of the model
%! @item M_
%! @item DynareModel
%! Matlab's structure describing the model (initialized by @code{dynare}).
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
......@@ -41,7 +41,7 @@ function [dr,info] = dyn_first_order_solver(jacobia,M_,dr,options,task)
%! @item info==5
%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
%! @item info==7
%! One of the generalized eigenvalues is close to 0/0
%! One of the generalized eigenvalues is close to 0/0
%! @end table
%! @end table
%! @end deftypefn
......@@ -63,75 +63,107 @@ function [dr,info] = dyn_first_order_solver(jacobia,M_,dr,options,task)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
info = 0;
dr.ghx = [];
dr.ghu = [];
klen = M_.maximum_endo_lag+M_.maximum_endo_lead+1;
kstate = dr.kstate;
kad = dr.kad;
kae = dr.kae;
nstatic = dr.nstatic;
nfwrd = dr.nfwrd;
npred = dr.npred;
nboth = dr.nboth;
persistent reorder_jacobian_columns innovations_idx index_s index_m index_c index_p row_indx index_0m index_0p k1 k2 j3 j4
persistent ndynamic nstatic nfwrd npred nboth nd nyf n
if ~nargin
reorder_jacobian_columns = [];
dr = [];
info = [];
return
end
if isempty(reorder_jacobian_columns)
kstate = dr.kstate;
nfwrd = dr.nfwrd;
nboth = dr.nboth;
npred = dr.npred-nboth;
nstatic = dr.nstatic;
ndynamic = npred+nfwrd+nboth;
nyf = nfwrd + nboth;
n = ndynamic+nstatic;
k1 = 1:(npred+nboth);
k2 = 1:(nfwrd+nboth);
order_var = dr.order_var;
nd = size(kstate,1);
lead_lag_incidence = M_.lead_lag_incidence;
lead_lag_incidence = DynareModel.lead_lag_incidence;
nz = nnz(lead_lag_incidence);
sdyn = M_.endo_nbr - nstatic;
%lead variables actually present in the model
j4 = nstatic+npred+1:nstatic+npred+nboth+nfwrd; % Index on the forward and both variables
j3 = nonzeros(lead_lag_incidence(2,j4)) - nstatic - 2 * npred - nboth; % Index on the non-zeros forward and both variables
j4 = find(lead_lag_incidence(2,j4));
[junk,cols_b,cols_j] = find(lead_lag_incidence(M_.maximum_endo_lag+1,...
order_var));
if nstatic > 0
[Q,R] = qr(jacobia(:,cols_j(1:nstatic)));
aa = Q'*jacobia;
else
aa = jacobia;
end
k1 = find([1:klen] ~= M_.maximum_endo_lag+1);
a = aa(:,nonzeros(lead_lag_incidence(k1,:)'));
b(:,cols_b) = aa(:,cols_j);
b10 = b(1:nstatic,1:nstatic);
b11 = b(1:nstatic,nstatic+1:end);
b2 = b(nstatic+1:end,nstatic+1:end);
if any(isinf(a(:)))
info = 1;
return
end
no_lead_id = find(lead_lag_incidence(3,:)==0);
no_lag_id = find(lead_lag_incidence(1,:)==0);
% buildind D and E
d = zeros(nd,nd) ;
e = d ;
k = find(kstate(:,2) >= M_.maximum_endo_lag+2 & kstate(:,3));
d(1:sdyn,k) = a(nstatic+1:end,kstate(k,3)) ;
k1 = find(kstate(:,2) == M_.maximum_endo_lag+2);
e(1:sdyn,k1) = -b2(:,kstate(k1,1)-nstatic);
k = find(kstate(:,2) <= M_.maximum_endo_lag+1 & kstate(:,4));
e(1:sdyn,k) = -a(nstatic+1:end,kstate(k,4)) ;
k2 = find(kstate(:,2) == M_.maximum_endo_lag+1);
k2 = k2(~ismember(kstate(k2,1),kstate(k1,1)));
d(1:sdyn,k2) = b2(:,kstate(k2,1)-nstatic);
if ~isempty(kad)
for j = 1:size(kad,1)
d(sdyn+j,kad(j)) = 1 ;
e(sdyn+j,kae(j)) = 1 ;
end
static_id = intersect(no_lead_id,no_lag_id);
lag_id = setdiff(no_lead_id,static_id);
lead_id = setdiff(no_lag_id,static_id);
both_id = intersect(setdiff(1:n,no_lead_id),setdiff(1:n,no_lag_id));
lead_idx = lead_lag_incidence(3,lead_id);
lag_idx = lead_lag_incidence(1,lag_id);
both_lagged_idx = lead_lag_incidence(1,both_id);
both_leaded_idx = lead_lag_incidence(3,both_id);
innovations_idx = (size(jacobia,2)-DynareModel.exo_nbr+1):size(jacobia,2);
dr.state_var = [lag_idx, both_lagged_idx];
indexi_0 = 0;
if DynareModel.maximum_endo_lag > 0 && (npred > 0 || nboth > 0)
indexi_0 = min(lead_lag_incidence(2,:));
end
% 1) if mjdgges.dll (or .mexw32 or ....) doesn't exit,
% matlab/qz is added to the path. There exists now qz/mjdgges.m that
% contains the calls to the old Sims code
% 2) In global_initialization.m, if mjdgges.m is visible exist(...)==2,
% this means that the DLL isn't avaiable and use_qzdiv is set to 1
[err,ss,tt,w,sdim,dr.eigval,info1] = mjdgges(e,d,options.qz_criterium);
index_c = lead_lag_incidence(2,:); % Index of all endogenous variables present at time=t
index_s = lead_lag_incidence(2,1:nstatic); % Index of all static endogenous variables present at time=t
index_0m = (nstatic+1:nstatic+npred)+indexi_0-1;
index_0p = (nstatic+npred+1:n)+indexi_0-1;
index_m = 1:(npred+nboth);
index_p = lead_lag_incidence(3,find(lead_lag_incidence(3,:)));
row_indx = nstatic+1:n;
reorder_jacobian_columns = [lag_idx, both_lagged_idx, npred+nboth+[static_id lag_id both_id lead_id], both_leaded_idx, lead_idx, innovations_idx ];
end
info = 0;
dr.ghx = [];
dr.ghu = [];
jacobia = jacobia(:,reorder_jacobian_columns);
if nstatic > 0
[Q, junk] = qr(jacobia(:,index_s));
aa = Q'*jacobia;
else
aa = jacobia;
end
A = aa(:,index_m); % Jacobain matrix for lagged endogeneous variables
B = aa(:,index_c); % Jacobian matrix for contemporaneous endogeneous variables
C = aa(:,index_p); % Jacobain matrix for led endogeneous variables
if task ~= 1 && DynareOptions.dr_cycle_reduction == 1
A1 = [aa(row_indx,index_m ) zeros(ndynamic,nfwrd)];
B1 = [aa(row_indx,index_0m) aa(row_indx,index_0p) ];
C1 = [zeros(ndynamic,npred) aa(row_indx,index_p)];
[ghx, info] = cycle_reduction(A1, B1, C1, DynareOptions.dr_cycle_reduction_tol);
ghx = ghx(:,index_m);
hx = ghx(1:npred+nboth,:);
gx = ghx(1+npred:end,:);
end
if (task ~= 1 && ((DynareOptions.dr_cycle_reduction == 1 && info ==1) || DynareOptions.dr_cycle_reduction == 0)) || task == 1
D = [[aa(row_indx,index_0m) zeros(ndynamic,nboth) aa(row_indx,index_p)] ; [zeros(nboth, npred) eye(nboth) zeros(nboth, nboth + nfwrd)]];
E = [-aa(row_indx,[index_m index_0p]) ; [zeros(nboth,nboth+npred) eye(nboth,nboth+nfwrd) ] ];
[err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E,D,DynareOptions.qz_criterium);
mexErrCheck('mjdgges', err);
if info1
......@@ -141,21 +173,19 @@ function [dr,info] = dyn_first_order_solver(jacobia,M_,dr,options,task)
else
info(1) = 2;
info(2) = info1;
info(3) = size(e,2);
info(3) = size(E,2);
end
return
end
nba = nd-sdim;
nyf = sum(kstate(:,2) > M_.maximum_endo_lag+1);
if task == 1
dr.rank = rank(w(1:nyf,nd-nyf+1:end));
% Under Octave, eig(A,B) doesn't exist, and
% lambda = qz(A,B) won't return infinite eigenvalues
if ~exist('OCTAVE_VERSION')
dr.eigval = eig(e,d);
dr.eigval = eig(E,D);
end
return
end
......@@ -163,74 +193,86 @@ function [dr,info] = dyn_first_order_solver(jacobia,M_,dr,options,task)
if nba ~= nyf
temp = sort(abs(dr.eigval));
if nba > nyf
temp = temp(nd-nba+1:nd-nyf)-1-options.qz_criterium;
temp = temp(nd-nba+1:nd-nyf)-1-DynareOptions.qz_criterium;
info(1) = 3;
elseif nba < nyf;
temp = temp(nd-nyf+1:nd-nba)-1-options.qz_criterium;
temp = temp(nd-nyf+1:nd-nba)-1-DynareOptions.qz_criterium;
info(1) = 4;
end
info(2) = temp'*temp;
return
end
np = nd - nyf;
n2 = np + 1;
n3 = nyf;
n4 = n3 + 1;
% derivatives with respect to dynamic state variables
% forward variables
w1 =w(1:n3,n2:nd);
if ~isscalar(w1) && (condest(w1) > 1e9);
% condest() fails on a scalar under Octave
info(1) = 5;
info(2) = condest(w1);
return;
else
gx = -w1'\w(n4:nd,n2:nd)';
end
% predetermined variables
hx = w(1:n3,1:np)'*gx+w(n4:nd,1:np)';
hx = (tt(1:np,1:np)*hx)\(ss(1:np,1:np)*hx);
%First order approximation
if task ~= 1
indx_stable_root = 1: (nd - nyf); %=> index of stable roots
indx_explosive_root = npred + nboth + 1:nd; %=> index of explosive roots
% derivatives with respect to dynamic state variables
% forward variables
Z = w';
Z11t = Z(indx_stable_root, indx_stable_root)';
Z21 = Z(indx_explosive_root, indx_stable_root);
Z22 = Z(indx_explosive_root, indx_explosive_root);
if ~isfloat(Z21) && (condest(Z21) > 1e9)
info(1) = 5;
info(2) = condest(Z21);
return;
else
gx = - Z22 \ Z21;
end
% predetermined variables
hx = Z11t * inv(tt(indx_stable_root, indx_stable_root)) * ss(indx_stable_root, indx_stable_root) * inv(Z11t);
ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];
end;
end
k1 = find(kstate(n4:nd,2) == M_.maximum_endo_lag+1);
k2 = find(kstate(1:n3,2) == M_.maximum_endo_lag+2);
dr.gx = gx;
dr.ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];
if task~= 1
%lead variables actually present in the model
j3 = nonzeros(kstate(:,3));
j4 = find(kstate(:,3));
% derivatives with respect to exogenous variables
if M_.exo_nbr
fu = aa(:,nz+(1:M_.exo_nbr));
a1 = b;
aa1 = [];
if nstatic > 0
aa1 = a1(:,1:nstatic);
end
dr.ghu = -[aa1 a(:,j3)*gx(j4,1:npred)+a1(:,nstatic+1:nstatic+ ...
npred) a1(:,nstatic+npred+1:end)]\fu;
if nstatic > 0
B_static = B(:,1:nstatic); % submatrix containing the derivatives w.r. to static variables
else
dr.ghu = [];
end
B_static = [];
end;
%static variables, backward variable, mixed variables and forward variables
B_pred = B(:,nstatic+1:nstatic+npred+nboth);
B_fyd = B(:,nstatic+npred+nboth+1:end);
% static variables
if nstatic > 0
temp = -a(1:nstatic,j3)*gx(j4,:)*hx;
j5 = find(kstate(n4:nd,4));
temp(:,j5) = temp(:,j5)-a(1:nstatic,nonzeros(kstate(:,4)));
temp = b10\(temp-b11*dr.ghx);
dr.ghx = [temp; dr.ghx];
temp = - C(1:nstatic,j3)*gx(j4,:)*hx;
b = aa(:,index_c);
b10 = b(1:nstatic, 1:nstatic);
b11 = b(1:nstatic, nstatic+1:n);
temp(:,index_m) = temp(:,index_m)-A(1:nstatic,:);
temp = b10\(temp-b11*ghx);
ghx = [temp; ghx];
temp = [];
end
if options.use_qzdiv
%% Necessary when using Sims' routines for QZ
gx = real(gx);
hx = real(hx);
dr.ghx = real(dr.ghx);
dr.ghu = real(dr.ghu);
end
A_ = real([B_static C(:,j3)*gx+B_pred B_fyd]); % The state_variable of the block are located at [B_pred B_both]
if DynareModel.exo_nbr
if nstatic > 0
fu = Q' * jacobia(:,innovations_idx);
else
fu = jacobia(:,innovations_idx);
end;
ghu = - A_ \ fu;
else
ghu = [];
end;
end
dr.ghx = ghx;
dr.ghu = ghu;
if DynareOptions.aim_solver ~= 1 && DynareOptions.use_qzdiv
% Necessary when using Sims' routines for QZ
dr.ghx = real(ghx);
dr.ghu = real(ghu);
hx = real(hx);
end
dr.Gy = hx;
\ No newline at end of file
dr.Gy = hx;
\ No newline at end of file
Supports Markdown
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