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

Use sparse representation for evaluating the static model 

We also take advantage of the fact that the non-block version is always
available next to the block one, so when we are only interested in the residual
as a whole, we simplify by using the non-block version.
estimate-initial-state
Sébastien Villemot 2023-01-10 16:05:33 +01:00
parent b343f0231d
commit 3c55aa57e1
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
5 changed files with 37 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
matlab/dyn_ramsey_static.m
View File

@ -17,7 +17,7 @@ function [steady_state, params, check] = dyn_ramsey_static(ys_init, M, options_,
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2003-2022 Dynare Team % Copyright © 2003-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -40,6 +40,7 @@ check = 0;
options_.steadystate.nocheck = 1; %locally disable checking because Lagrange multipliers are not accounted for in evaluate_steady_state_file options_.steadystate.nocheck = 1; %locally disable checking because Lagrange multipliers are not accounted for in evaluate_steady_state_file
% dyn_ramsey_static_1 is a subfunction % dyn_ramsey_static_1 is a subfunction
nl_func = @(x) dyn_ramsey_static_1(x,M,options_,oo); nl_func = @(x) dyn_ramsey_static_1(x,M,options_,oo);
exo_ss = [oo.exo_steady_state oo.exo_det_steady_state];
% check_static_model is a subfunction % check_static_model is a subfunction
if ~options_.steadystate_flag && check_static_model(ys_init,M,options_,oo) if ~options_.steadystate_flag && check_static_model(ys_init,M,options_,oo)
@ -72,7 +73,6 @@ elseif options_.steadystate_flag
end end
end end
ys_init(k_inst) = inst_val; ys_init(k_inst) = inst_val;
exo_ss = [oo.exo_steady_state oo.exo_det_steady_state];
[xx,params] = evaluate_steady_state_file(ys_init,exo_ss,M,options_,~options_.steadystate.nocheck); %run steady state file again to update parameters [xx,params] = evaluate_steady_state_file(ys_init,exo_ss,M,options_,~options_.steadystate.nocheck); %run steady state file again to update parameters
[~,~,steady_state] = nl_func(inst_val); %compute and return steady state [~,~,steady_state] = nl_func(inst_val); %compute and return steady state
else else
@ -107,6 +107,7 @@ inst_nbr = orig_endo_aux_nbr - orig_eq_nbr;
% indices of Lagrange multipliers % indices of Lagrange multipliers
fname = M.fname; fname = M.fname;
exo_ss = [oo.exo_steady_state oo.exo_det_steady_state];
if options_.steadystate_flag if options_.steadystate_flag
k_inst = []; k_inst = [];
@ -117,8 +118,7 @@ if options_.steadystate_flag
ys_init=zeros(size(oo.steady_state)); %create starting vector for steady state computation as only instrument value is handed over ys_init=zeros(size(oo.steady_state)); %create starting vector for steady state computation as only instrument value is handed over
ys_init(k_inst) = x; %set instrument, the only value required for steady state computation, to current value ys_init(k_inst) = x; %set instrument, the only value required for steady state computation, to current value
[x,params,check] = evaluate_steady_state_file(ys_init,... %returned x now has size endo_nbr as opposed to input size of n_instruments [x,params,check] = evaluate_steady_state_file(ys_init,... %returned x now has size endo_nbr as opposed to input size of n_instruments
[oo.exo_steady_state; ... exo_ss, ...
oo.exo_det_steady_state], ...
M,options_,~options_.steadystate.nocheck); M,options_,~options_.steadystate.nocheck);
if any(imag(x(1:M.orig_endo_nbr))) %return with penalty if any(imag(x(1:M.orig_endo_nbr))) %return with penalty
resids=ones(inst_nbr,1)+sum(abs(imag(x(1:M.orig_endo_nbr)))); %return with penalty resids=ones(inst_nbr,1)+sum(abs(imag(x(1:M.orig_endo_nbr)))); %return with penalty
@ -141,10 +141,7 @@ if any([M.aux_vars.type] ~= 6) %auxiliary variables other than multipliers
needs_set_auxiliary_variables = 1; needs_set_auxiliary_variables = 1;
if M.set_auxiliary_variables if M.set_auxiliary_variables
fh = str2func([M.fname '.set_auxiliary_variables']); fh = str2func([M.fname '.set_auxiliary_variables']);
s_a_v_func = @(z) fh(z,... s_a_v_func = @(z) fh(z, exo_ss, params);
[oo.exo_steady_state,...
oo.exo_det_steady_state],...
params);
else else
s_a_v_func = z; s_a_v_func = z;
end end
@ -156,12 +153,11 @@ end
% set multipliers and auxiliary variables that % set multipliers and auxiliary variables that
% depends on multipliers to 0 to compute residuals % depends on multipliers to 0 to compute residuals
if (options_.bytecode) if (options_.bytecode)
[res, junk] = bytecode('static',xx,[oo.exo_steady_state oo.exo_det_steady_state], ... [res, junk] = bytecode('static', xx, exo_ss, params, 'evaluate');
params, 'evaluate');
fJ = junk.g1; fJ = junk.g1;
else else
[res,fJ] = feval([fname '.static'],xx,[oo.exo_steady_state oo.exo_det_steady_state], ... [res, T_order, T] = feval([fname '.sparse.static_resid'], xx, exo_ss, params);
params); fJ = feval([fname '.sparse.static_g1'], xx, exo_ss, params, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, T_order, T);
end end
% index of multipliers and corresponding equations % index of multipliers and corresponding equations
% the auxiliary variables before the Lagrange multipliers are treated % the auxiliary variables before the Lagrange multipliers are treated
@ -192,12 +188,11 @@ end
function result = check_static_model(ys,M,options_,oo) function result = check_static_model(ys,M,options_,oo)
result = false; result = false;
exo_ss = [oo.exo_steady_state oo.exo_det_steady_state];
if (options_.bytecode) if (options_.bytecode)
[res, ~] = bytecode('static',ys,[oo.exo_steady_state oo.exo_det_steady_state], ... [res, ~] = bytecode('static', ys, exo_ss, M.params, 'evaluate');
M.params, 'evaluate');
else else
res = feval([M.fname '.static'],ys,[oo.exo_steady_state oo.exo_det_steady_state], ... res = feval([M.fname '.sparse.static_resid'], ys, exo_ss, M.params);
M.params);
end end
if norm(res) < options_.solve_tolf if norm(res) < options_.solve_tolf
result = true; result = true;

27
matlab/evaluate_static_model.m
View File

@ -20,7 +20,7 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2001-2021 Dynare Team % Copyright © 2001-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -48,27 +48,8 @@ if options.bytecode
jacob = junk.g1; jacob = junk.g1;
end end
else else
fh_static = str2func([M.fname '.static']); [residuals, T_order, T] = feval([M.fname '.sparse.static_resid'], ys, exo_ss, params);
if options.block if nargout >= 3
residuals = zeros(M.endo_nbr,1); jacob = feval([M.fname '.sparse.static_g1'], ys, exo_ss, params, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, T_order, T);
T = NaN(M.block_structure_stat.tmp_nbr, 1);
for b = 1:length(M.block_structure_stat.block)
[r, yy, T] = feval(fh_static,b,ys,exo_ss,params,T);
if M.block_structure_stat.block(b).Simulation_Type == 1 || ... % evaluateForward
M.block_structure_stat.block(b).Simulation_Type == 2 % evaluateBackward
vidx = M.block_structure_stat.block(b).variable;
r = yy(vidx) - ys(vidx);
end
residuals(M.block_structure_stat.block(b).equation) = r;
end
if nargout==3
jacob=NaN(length(ys));
end
else
if nargout<3
residuals = feval(fh_static,ys,exo_ss,params);
else
[residuals, jacob] = feval(fh_static,ys,exo_ss,params);
end
end end
end end

50
matlab/evaluate_steady_state.m
View File

@ -22,7 +22,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2001-2022 Dynare Team % Copyright © 2001-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -279,8 +279,9 @@ elseif steadystate_flag
return return
end end
elseif ~options.bytecode && ~options.block elseif ~options.bytecode && ~options.block
static_resid = str2func(sprintf('%s.sparse.static_resid', M.fname));
static_g1 = str2func(sprintf('%s.sparse.static_g1', M.fname));
if ~options.linear if ~options.linear
static_model = str2func(sprintf('%s.static', M.fname));
% non linear model % non linear model
if ismember(options.solve_algo,[10,11]) if ismember(options.solve_algo,[10,11])
[lb,ub,eq_index] = get_complementarity_conditions(M,options.ramsey_policy); [lb,ub,eq_index] = get_complementarity_conditions(M,options.ramsey_policy);
@ -295,16 +296,18 @@ elseif ~options.bytecode && ~options.block
ys_init,... ys_init,...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ... options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params,... options, exo_ss, params,...
M.endo_nbr,static_model,eq_index); M.endo_nbr, static_resid, static_g1, ...
M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, eq_index);
else else
[ys, check] = dynare_solve(@static_problem, ys_init, ... [ys, check] = dynare_solve(@static_problem, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ... options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params, M.endo_nbr, static_model); options, exo_ss, params, M.endo_nbr, static_resid, static_g1, ...
M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr);
end end
if check && options.debug if check && options.debug
[ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ... [ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ... options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params, M.endo_nbr, static_model); options, exo_ss, params, M.endo_nbr, static_resid, static_g1, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr);
dprintf('Nonlinear solver routine returned errorcode=%i.', errorcode) dprintf('Nonlinear solver routine returned errorcode=%i.', errorcode)
skipline() skipline()
[infrow,infcol]=find(isinf(fjac) | isnan(fjac)); [infrow,infcol]=find(isinf(fjac) | isnan(fjac));
@ -323,9 +326,8 @@ elseif ~options.bytecode && ~options.block
end end
else else
% linear model % linear model
fh_static = str2func([M.fname '.static']); [fvec, T_order, T] = static_resid(ys_init, exo_ss, params);
[fvec,jacob] = fh_static(ys_init,exo_ss, ... jacob = static_g1(ys_init, exo_ss, params, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, T_order, T);
params);
ii = find(~isfinite(fvec)); ii = find(~isfinite(fvec));
if ~isempty(ii) if ~isempty(ii)
@ -382,27 +384,12 @@ end
if M.static_and_dynamic_models_differ if M.static_and_dynamic_models_differ
% Evaluate residual of *dynamic* model using the steady state % Evaluate residual of *dynamic* model using the steady state
% computed on the *static* one % computed on the *static* one
z = repmat(ys,1,M.maximum_lead + M.maximum_lag + 1);
zx = repmat([exo_ss'], M.maximum_lead + M.maximum_lag + 1, 1);
if options.bytecode if options.bytecode
z = repmat(ys,1,M.maximum_lead + M.maximum_lag + 1);
zx = repmat([exo_ss'], M.maximum_lead + M.maximum_lag + 1, 1);
[r, ~]= bytecode('dynamic','evaluate', z, zx, params, ys, 1); [r, ~]= bytecode('dynamic','evaluate', z, zx, params, ys, 1);
elseif options.block
T=NaN(M.block_structure.dyn_tmp_nbr, 1);
for i = 1:length(M.block_structure.block)
[rr, yy, T, g] = feval([M.fname '.dynamic'], i, ...
dynvars_from_endo_simul(z, M.maximum_lag+1, M), ...
zx, params, ys, T, M.maximum_lag+1, false);
if M.block_structure.block(i).Simulation_Type == 1 || ... % evaluateForward
M.block_structure.block(i).Simulation_Type == 2 % evaluateBackward
vidx = M.block_structure.block(i).variable;
rr = yy(M.lead_lag_incidence(M.maximum_endo_lag+1, vidx)) - oo.steady_state(vidx);
end
idx = M.block_structure.block(i).equation;
r(idx) = rr;
end
else else
r = feval([M.fname '.dynamic'], dynvars_from_endo_simul(z, M.maximum_lag+1, M), ... r = feval([M.fname '.sparse.dynamic_resid'], repmat(ys, 3, 1), exo_ss, params, ys);
zx, params, ys, M.maximum_lag + 1);
end end
% Fail if residual greater than tolerance % Fail if residual greater than tolerance
if max(abs(r)) > options.solve_tolf if max(abs(r)) > options.solve_tolf
@ -424,13 +411,14 @@ if ~isempty(find(isnan(ys)))
return return
end end
function [resids,jac] = static_problem(y,x,params,nvar,fh_static_model) function [resids,jac] = static_problem(y, x, params, nvar, fh_static_resid, fh_static_g1, sparse_rowval, sparse_colval, sparse_colptr)
[r,j] = fh_static_model(y,x,params); [r, T_order, T] = fh_static_resid(y, x, params);
j = fh_static_g1(y, x, params, sparse_rowval, sparse_colval, sparse_colptr, T_order, T);
resids = r(1:nvar); resids = r(1:nvar);
jac = j(1:nvar,1:nvar); jac = j(1:nvar,1:nvar);
function [resids,jac] = static_mcp_problem(y, x, params, nvar, fh_static_resid, fh_static_g1, sparse_rowval, sparse_colval, sparse_colptr, eq_index)
function [resids,jac] = static_mcp_problem(y,x,params,nvar,fh_static_model,eq_index) [r, T_order, T] = fh_static_resid(y, x, params);
[r,j] = fh_static_model(y,x,params); j = fh_static_g1(y, x, params, sparse_rowval, sparse_colval, sparse_colptr, T_order, T);
resids = r(eq_index); resids = r(eq_index);
jac = j(eq_index,1:nvar); jac = j(eq_index,1:nvar);

3
matlab/model_diagnostics.m
View File

@ -157,7 +157,8 @@ for b=1:nb
int2str(b)]); int2str(b)]);
end end
else else
[res,jacob]=feval([M.fname '.static'],dr.ys,exo,M.params); [res, T_order, T] = feval([M.fname '.sparse.static_resid'], dr.ys, exo, M.params);
jacob = feval([M.fname '.sparse.static_g1'], dr.ys, exo, M.params, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, T_order, T);
end end
if any(any(isinf(jacob) | isnan(jacob))) if any(any(isinf(jacob) | isnan(jacob)))
problem_dummy=1; problem_dummy=1;

2
preprocessor

@ -1 +1 @@
Subproject commit a4f299c0b781676afe3d9d073a7e2cd0551af997 Subproject commit 35ac73fad8645302c10d0cd0068710ad3847366a