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

Implement even more naming consistency 

Mostly removes M for M_
kalman-mex
Johannes Pfeifer 2023-10-25 17:29:55 +02:00
parent 0295e5ede8
commit 0561200f1c
34 changed files with 826 additions and 820 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
matlab/AIM_first_order_solver.m
View File

@ -1,7 +1,7 @@
function [dr,info]=AIM_first_order_solver(jacobia,M,dr,qz_criterium)
function [dr,info]=AIM_first_order_solver(jacobia,M_,dr,qz_criterium)
%@info:
%! @deftypefn {Function File} {[@var{dr},@var{info}] =} AIM_first_order_solver (@var{jacobia},@var{M},@var{dr},@var{qz_criterium})
%! @deftypefn {Function File} {[@var{dr},@var{info}] =} AIM_first_order_solver (@var{jacobia},@var{M_},@var{dr},@var{qz_criterium})
%! @anchor{AIM_first_order_solver}
%! @sp 1
%! Computes the first order reduced form of the DSGE model using AIM.
@ -11,7 +11,7 @@ function [dr,info]=AIM_first_order_solver(jacobia,M,dr,qz_criterium)
%! @table @ @var
%! @item jacobia
%! Matrix containing the Jacobian of the model
%! @item M
%! @item M_
%! Matlab's structure describing the model (initialized by @code{dynare}).
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
@ -70,20 +70,20 @@ function [dr,info]=AIM_first_order_solver(jacobia,M,dr,qz_criterium)
info = 0;
[dr,aimcode]=dynAIMsolver1(jacobia,M,dr);
[dr,aimcode]=dynAIMsolver1(jacobia,M_,dr);
if aimcode ~=1
info(1) = convertAimCodeToInfo(aimCode); %convert to be in the 100 range
info(2) = 1.0e+8;
return
end
A = kalman_transition_matrix(dr,M.nstatic+(1:M.nspred), 1:M.nspred);
A = kalman_transition_matrix(dr,M_.nstatic+(1:M_.nspred), 1:M_.nspred);
dr.eigval = eig(A);
disp(dr.eigval)
nd = size(dr.kstate,1);
nba = nd-sum( abs(dr.eigval) < qz_criterium );
nsfwrd = M.nsfwrd;
nsfwrd = M_.nsfwrd;
if nba ~= nsfwrd
temp = sort(abs(dr.eigval));

78
matlab/check_for_calibrated_covariances.m
View File

@ -1,12 +1,12 @@
function estim_params=check_for_calibrated_covariances(estim_params,M_)
% function check_for_calibrated_covariances(estim_params,M)
function estim_params_=check_for_calibrated_covariances(estim_params_,M_)
% function check_for_calibrated_covariances(estim_params_,M)
% find calibrated covariances to consider during estimation
% Inputs
% -estim_params [structure] describing parameters to be estimated
% -estim_params_ [structure] describing parameters to be estimated
% -M_ [structure] describing the model
%
% Outputs
% -estim_params [structure] describing parameters to be estimated
% -estim_params_ [structure] describing parameters to be estimated
%
% Notes: M is local to this function and not updated when calling
% set_all_parameters
@ -28,19 +28,19 @@ function estim_params=check_for_calibrated_covariances(estim_params,M_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if isfield(estim_params,'calibrated_covariances')
estim_params = rmfield(estim_params,'calibrated_covariances'); %remove if already present
if isfield(estim_params_,'calibrated_covariances')
estim_params_ = rmfield(estim_params_,'calibrated_covariances'); %remove if already present
end
if isfield(estim_params,'calibrated_covariances_ME')
estim_params = rmfield(estim_params,'calibrated_covariances_ME'); %remove if already present
if isfield(estim_params_,'calibrated_covariances_ME')
estim_params_ = rmfield(estim_params_,'calibrated_covariances_ME'); %remove if already present
end
[rows_calibrated, columns_calibrated]=ind2sub(size(M_.Sigma_e),find(tril(M_.Sigma_e,-1))); %find linear indices of preset lower triangular covariance entries
if estim_params.ncx %delete preset entries actually estimated
for i=1:estim_params.ncx
shock_1 = estim_params.corrx(i,1);
shock_2 = estim_params.corrx(i,2);
if estim_params_.ncx %delete preset entries actually estimated
for i=1:estim_params_.ncx
shock_1 = estim_params_.corrx(i,1);
shock_2 = estim_params_.corrx(i,2);
estimated_corr_pos=find(rows_calibrated==shock_1 & columns_calibrated==shock_2);
if ~isempty(estimated_corr_pos)
rows_calibrated(estimated_corr_pos)=[];
@ -53,17 +53,17 @@ if estim_params.ncx %delete preset entries actually estimated
end
end
if any(rows_calibrated)
estim_params.calibrated_covariances.position=[sub2ind(size(M_.Sigma_e),rows_calibrated,columns_calibrated);sub2ind(size(M_.Sigma_e),columns_calibrated,rows_calibrated)]; %get linear entries of upper triangular parts
estim_params.calibrated_covariances.cov_value=M_.Sigma_e(estim_params.calibrated_covariances.position);
estim_params_.calibrated_covariances.position=[sub2ind(size(M_.Sigma_e),rows_calibrated,columns_calibrated);sub2ind(size(M_.Sigma_e),columns_calibrated,rows_calibrated)]; %get linear entries of upper triangular parts
estim_params_.calibrated_covariances.cov_value=M_.Sigma_e(estim_params_.calibrated_covariances.position);
end
end
[rows_calibrated, columns_calibrated]=ind2sub(size(M_.H),find(tril(M_.H,-1))); %find linear indices of preset lower triangular covariance entries
if estim_params.ncn %delete preset entries actually estimated
for i=1:estim_params.ncn
shock_1 = estim_params.corrn(i,1);
shock_2 = estim_params.corrn(i,2);
if estim_params_.ncn %delete preset entries actually estimated
for i=1:estim_params_.ncn
shock_1 = estim_params_.corrn(i,1);
shock_2 = estim_params_.corrn(i,2);
estimated_corr_pos=find(rows_calibrated==shock_1 & columns_calibrated==shock_2);
if ~isempty(estimated_corr_pos)
rows_calibrated(estimated_corr_pos)=[];
@ -77,8 +77,8 @@ if estim_params.ncn %delete preset entries actually estimated
end
end
if any(rows_calibrated)
estim_params.calibrated_covariances_ME.position=[sub2ind(size(M_.H),rows_calibrated,columns_calibrated);sub2ind(size(M_.H),columns_calibrated,rows_calibrated)]; %get linear entries of upper triangular parts
estim_params.calibrated_covariances_ME.cov_value=M_.H(estim_params.calibrated_covariances_ME.position);
estim_params_.calibrated_covariances_ME.position=[sub2ind(size(M_.H),rows_calibrated,columns_calibrated);sub2ind(size(M_.H),columns_calibrated,rows_calibrated)]; %get linear entries of upper triangular parts
estim_params_.calibrated_covariances_ME.cov_value=M_.H(estim_params_.calibrated_covariances_ME.position);
end
return % --*-- Unit tests --*--
@ -89,14 +89,14 @@ M_.Sigma_e=[1 0; 0 1];
M_.H=[1 0; 0 1];
M_.Correlation_matrix= [1 -0.5; -0.5 1];
M_.Correlation_matrix_ME=[1 -0.5; -0.5 1];
estim_params.ncx=1;
estim_params.ncn=1;
estim_params_.ncx=1;
estim_params_.ncn=1;
estim_params.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params.corrn=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params_.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params_.corrn=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params=check_for_calibrated_covariances(estim_params,M_);
if isfield(estim_params,'calibrated_covariances_ME') || isfield(estim_params,'calibrated_covariances')
estim_params_=check_for_calibrated_covariances(estim_params_,M_);
if isfield(estim_params_,'calibrated_covariances_ME') || isfield(estim_params_,'calibrated_covariances')
t(1)=false;
else
t(1)=true;
@ -106,26 +106,26 @@ M_.Sigma_e=[1 -0.1; -0.1 1];
M_.H=[1 -0.1; -0.1 1];
M_.Correlation_matrix= [1 -0.5; -0.5 1];
M_.Correlation_matrix_ME=[1 0; 0 1];
estim_params.ncx=1;
estim_params.ncn=0;
estim_params_.ncx=1;
estim_params_.ncn=0;
estim_params.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params.corrn=[];
estim_params=check_for_calibrated_covariances(estim_params,M_);
t(2)=isequal(estim_params.calibrated_covariances_ME.position,[2;3]);
t(3)=isequal(estim_params.calibrated_covariances_ME.cov_value,[-0.1;-0.1]);
estim_params_.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params_.corrn=[];
estim_params_=check_for_calibrated_covariances(estim_params_,M_);
t(2)=isequal(estim_params_.calibrated_covariances_ME.position,[2;3]);
t(3)=isequal(estim_params_.calibrated_covariances_ME.cov_value,[-0.1;-0.1]);
M_.Sigma_e=[1 -0.1; -0.1 1];
M_.H=[1 -0.1; -0.1 1];
M_.Correlation_matrix= [1 -0.5; -0.5 1];
M_.Correlation_matrix_ME=[1 0; 0 1];
estim_params.ncx=1;
estim_params.ncn=1;
estim_params_.ncx=1;
estim_params_.ncn=1;
estim_params.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params.corrn=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params=check_for_calibrated_covariances(estim_params,M_);
if isfield(estim_params,'calibrated_covariances_ME') || isfield(estim_params,'calibrated_covariances')
estim_params_.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params_.corrn=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params_=check_for_calibrated_covariances(estim_params_,M_);
if isfield(estim_params_,'calibrated_covariances_ME') || isfield(estim_params_,'calibrated_covariances')
t(4)=false;
else
t(4)=true;

70
matlab/dyn_ramsey_static.m
View File

@ -1,5 +1,5 @@
function [steady_state, params, check] = dyn_ramsey_static(ys_init, exo_ss, M, options_)
function [steady_state, params, check] = dyn_ramsey_static(ys_init, exo_ss, M_, options_)
% [steady_state, params, check] = dyn_ramsey_static(ys_init, exo_ss, M_, options_)
% Computes the steady state for optimal policy
%
% When there is no steady state file, relies on the fact that Lagrange
@ -9,11 +9,11 @@ function [steady_state, params, check] = dyn_ramsey_static(ys_init, exo_ss, M, o
% the multipliers that minimizes the residuals, given the other variables
% (using a minimum norm solution, easy to compute because of the linearity
% property).
%
% INPUTS
% ys_init: vector of endogenous variables or instruments
% exo_ss vector of exogenous steady state (incl. deterministic exogenous)
% M: Dynare model structure
% M_: Dynare model structure
% options: Dynare options structure
%
% OUTPUTS
@ -43,19 +43,19 @@ function [steady_state, params, check] = dyn_ramsey_static(ys_init, exo_ss, M, o
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
params = M.params;
params = M_.params;
check = 0;
options_.steadystate.nocheck = 1; %locally disable checking because Lagrange multipliers are not accounted for in evaluate_steady_state_file
nl_func = @(x) dyn_ramsey_static_1(x,exo_ss,ys_init,M,options_);
nl_func = @(x) dyn_ramsey_static_1(x,exo_ss,ys_init,M_,options_);
if ~options_.steadystate_flag && check_static_model(ys_init,exo_ss,M,options_)
if ~options_.steadystate_flag && check_static_model(ys_init,exo_ss,M_,options_)
steady_state = ys_init;
return
elseif options_.steadystate_flag
k_inst = [];
inst_nbr = size(options_.instruments,1);
for i = 1:inst_nbr
k_inst = [k_inst; strmatch(options_.instruments{i}, M.endo_names, 'exact')];
k_inst = [k_inst; strmatch(options_.instruments{i}, M_.endo_names, 'exact')];
end
if inst_nbr == 1
%solve for instrument, using univariate solver, starting at initial value for instrument
@ -78,10 +78,10 @@ elseif options_.steadystate_flag
end
end
ys_init(k_inst) = inst_val;
[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
else
xx = ys_init(1:M.orig_endo_nbr);
xx = ys_init(1:M_.orig_endo_nbr);
o_jacobian_flag = options_.jacobian_flag;
options_.jacobian_flag = false;
[xx, errorflag] = dynare_solve(nl_func, xx, options_.ramsey.maxit, options_.solve_tolf, options_.solve_tolx, options_);
@ -93,42 +93,42 @@ else
end
function [resids,rJ,steady_state] = dyn_ramsey_static_1(x,exo_ss,ys_init,M,options_)
function [resids,rJ,steady_state] = dyn_ramsey_static_1(x,exo_ss,ys_init,M_,options_)
resids = [];
rJ = [];
mult = [];
inst_nbr = M.ramsey_orig_endo_nbr - M.ramsey_orig_eq_nbr;
inst_nbr = M_.ramsey_orig_endo_nbr - M_.ramsey_orig_eq_nbr;
if options_.steadystate_flag
k_inst = [];
for i = 1:size(options_.instruments,1)
k_inst = [k_inst; strmatch(options_.instruments{i}, M.endo_names, 'exact')];
k_inst = [k_inst; strmatch(options_.instruments{i}, M_.endo_names, 'exact')];
end
ys_init(k_inst) = x; %set instrument, the only value required for steady state computation, to current value
[x,M.params,check] = evaluate_steady_state_file(ys_init,... %returned x now has size endo_nbr as opposed to input size of n_instruments
[x,M_.params,check] = evaluate_steady_state_file(ys_init,... %returned x now has size endo_nbr as opposed to input size of n_instruments
exo_ss, ...
M,options_,~options_.steadystate.nocheck);
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
steady_state=NaN(M.endo_nbr,1);
M_,options_,~options_.steadystate.nocheck);
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
steady_state=NaN(M_.endo_nbr,1);
return
end
if check(1) %return
resids=ones(inst_nbr,1)+sum(abs(x(1:M.orig_endo_nbr))); %return with penalty
steady_state=NaN(M.endo_nbr,1);
resids=ones(inst_nbr,1)+sum(abs(x(1:M_.orig_endo_nbr))); %return with penalty
steady_state=NaN(M_.endo_nbr,1);
return
end
end
xx = zeros(M.endo_nbr,1); %initialize steady state vector
xx(1:M.orig_endo_nbr) = x(1:M.orig_endo_nbr); %set values of original endogenous variables based on steady state file or initial value
xx = zeros(M_.endo_nbr,1); %initialize steady state vector
xx(1:M_.orig_endo_nbr) = x(1:M_.orig_endo_nbr); %set values of original endogenous variables based on steady state file or initial value
% Determine whether other auxiliary variables will need to be updated
if any([M.aux_vars.type] ~= 6) %auxiliary variables other than multipliers
if any([M_.aux_vars.type] ~= 6) %auxiliary variables other than multipliers
needs_set_auxiliary_variables = true;
fh = str2func([M.fname '.set_auxiliary_variables']);
s_a_v_func = @(z) fh(z, exo_ss, M.params);
fh = str2func([M_.fname '.set_auxiliary_variables']);
s_a_v_func = @(z) fh(z, exo_ss, M_.params);
xx = s_a_v_func(xx);
else
needs_set_auxiliary_variables = false;
@ -137,12 +137,12 @@ end
% Compute the value of the Lagrange multipliers that minimizes the norm of the
% residuals, given the other endogenous
if options_.bytecode
res = bytecode('static', M, options, xx, exo_ss, M.params, 'evaluate');
res = bytecode('static', M_, options, xx, exo_ss, M_.params, 'evaluate');
else
res = feval([M.fname '.sparse.static_resid'], xx, exo_ss, M.params);
res = feval([M_.fname '.sparse.static_resid'], xx, exo_ss, M_.params);
end
A = feval([M.fname '.ramsey_multipliers_static_g1'], xx, exo_ss, M.params, M.ramsey_multipliers_static_g1_sparse_rowval, M.ramsey_multipliers_static_g1_sparse_colval, M.ramsey_multipliers_static_g1_sparse_colptr);
y = res(1:M.ramsey_orig_endo_nbr);
A = feval([M_.fname '.ramsey_multipliers_static_g1'], xx, exo_ss, M_.params, M_.ramsey_multipliers_static_g1_sparse_rowval, M_.ramsey_multipliers_static_g1_sparse_colval, M_.ramsey_multipliers_static_g1_sparse_colptr);
y = res(1:M_.ramsey_orig_endo_nbr);
mult = -A\y;
resids1 = y+A*mult;
@ -156,20 +156,20 @@ end
if options_.steadystate_flag
resids = r1;
else
resids = [res(M.ramsey_orig_endo_nbr+(1:M.orig_endo_nbr-inst_nbr)); r1];
resids = [res(M_.ramsey_orig_endo_nbr+(1:M_.orig_endo_nbr-inst_nbr)); r1];
end
if needs_set_auxiliary_variables
steady_state = s_a_v_func([xx(1:M.ramsey_orig_endo_nbr); mult]);
steady_state = s_a_v_func([xx(1:M_.ramsey_orig_endo_nbr); mult]);
else
steady_state = [xx(1:M.ramsey_orig_endo_nbr); mult];
steady_state = [xx(1:M_.ramsey_orig_endo_nbr); mult];
end
function result = check_static_model(ys,exo_ss,M,options_)
function result = check_static_model(ys,exo_ss,M_,options_)
result = false;
if (options_.bytecode)
[res, ~] = bytecode('static', M, options, ys, exo_ss, M.params, 'evaluate');
[res, ~] = bytecode('static', M_, options, ys, exo_ss, M_.params, 'evaluate');
else
res = feval([M.fname '.sparse.static_resid'], ys, exo_ss, M.params);
res = feval([M_.fname '.sparse.static_resid'], ys, exo_ss, M_.params);
end
if norm(res) < options_.solve_tolf
result = true;

42
matlab/dyn_second_order_solver.m
View File

@ -1,4 +1,4 @@
function dr = dyn_second_order_solver(jacobia,hessian_mat,dr,M,threads_BC)
function dr = dyn_second_order_solver(jacobia,hessian_mat,dr,M_,threads_BC)
%@info:
%! @deftypefn {Function File} {@var{dr} =} dyn_second_order_solver (@var{jacobia},@var{hessian_mat},@var{dr},@var{M_},@var{threads_BC})
@ -58,33 +58,33 @@ dr.ghuu = [];
dr.ghxu = [];
dr.ghs2 = [];
k1 = nonzeros(M.lead_lag_incidence(:,dr.order_var)');
kk1 = [k1; length(k1)+(1:M.exo_nbr+M.exo_det_nbr)'];
k1 = nonzeros(M_.lead_lag_incidence(:,dr.order_var)');
kk1 = [k1; length(k1)+(1:M_.exo_nbr+M_.exo_det_nbr)'];
nk = size(kk1,1);
kk2 = reshape(1:nk^2,nk,nk);
ic = [ M.nstatic+(1:M.nspred) ]';
ic = [ M_.nstatic+(1:M_.nspred) ]';
klag = M.lead_lag_incidence(1,dr.order_var); %columns are in DR order
kcurr = M.lead_lag_incidence(2,dr.order_var); %columns are in DR order
klead = M.lead_lag_incidence(3,dr.order_var); %columns are in DR order
klag = M_.lead_lag_incidence(1,dr.order_var); %columns are in DR order
kcurr = M_.lead_lag_incidence(2,dr.order_var); %columns are in DR order
klead = M_.lead_lag_incidence(3,dr.order_var); %columns are in DR order
%% ghxx
A = zeros(M.endo_nbr,M.endo_nbr);
A = zeros(M_.endo_nbr,M_.endo_nbr);
A(:,kcurr~=0) = jacobia(:,nonzeros(kcurr));
A(:,ic) = A(:,ic) + jacobia(:,nonzeros(klead))*dr.ghx(klead~=0,:);
B = zeros(M.endo_nbr,M.endo_nbr);
B(:,M.nstatic+M.npred+1:end) = jacobia(:,nonzeros(klead));
B = zeros(M_.endo_nbr,M_.endo_nbr);
B(:,M_.nstatic+M_.npred+1:end) = jacobia(:,nonzeros(klead));
C = dr.ghx(ic,:);
zx = [eye(length(ic));
dr.ghx(kcurr~=0,:);
dr.ghx(klead~=0,:)*dr.ghx(ic,:);
zeros(M.exo_nbr,length(ic));
zeros(M.exo_det_nbr,length(ic))];
zu = [zeros(length(ic),M.exo_nbr);
zeros(M_.exo_nbr,length(ic));
zeros(M_.exo_det_nbr,length(ic))];
zu = [zeros(length(ic),M_.exo_nbr);
dr.ghu(kcurr~=0,:);
dr.ghx(klead~=0,:)*dr.ghu(ic,:);
eye(M.exo_nbr);
zeros(M.exo_det_nbr,M.exo_nbr)];
eye(M_.exo_nbr);
zeros(M_.exo_det_nbr,M_.exo_nbr)];
rhs = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,threads_BC); %hessian_mat: reordering to DR order
rhs = -rhs;
dr.ghxx = gensylv(2,A,B,C,rhs);
@ -108,17 +108,17 @@ dr.ghuu = A\rhs;
%% ghs2
% derivatives of F with respect to forward variables
O1 = zeros(M.endo_nbr,M.nstatic);
O2 = zeros(M.endo_nbr,M.nfwrd);
LHS = zeros(M.endo_nbr,M.endo_nbr);
O1 = zeros(M_.endo_nbr,M_.nstatic);
O2 = zeros(M_.endo_nbr,M_.nfwrd);
LHS = zeros(M_.endo_nbr,M_.endo_nbr);
LHS(:,kcurr~=0) = jacobia(:,nonzeros(kcurr));
RHS = zeros(M.endo_nbr,M.exo_nbr^2);
E = eye(M.endo_nbr);
RHS = zeros(M_.endo_nbr,M_.exo_nbr^2);
E = eye(M_.endo_nbr);
B1 = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(nonzeros(klead),nonzeros(klead))), dr.ghu(klead~=0,:),threads_BC); %hessian_mat:focus only on forward variables and reorder to DR order
RHS = RHS + jacobia(:,nonzeros(klead))*dr.ghuu(klead~=0,:)+B1;
% LHS
LHS = LHS + jacobia(:,nonzeros(klead))*(E(klead~=0,:)+[O1(klead~=0,:) dr.ghx(klead~=0,:) O2(klead~=0,:)]);
RHS = RHS*M.Sigma_e(:);
RHS = RHS*M_.Sigma_e(:);
dr.fuu = RHS;
RHS = -RHS;
dr.ghs2 = LHS\RHS;

34
matlab/dynare_estimation_init.m
View File

@ -335,29 +335,29 @@ if options_.analytic_derivation
options_.analytic_derivation = 1;
if estim_params_.np || isfield(options_,'identification_check_endogenous_params_with_no_prior')
% check if steady state changes param values
M=M_;
M_local=M_;
if isfield(options_,'identification_check_endogenous_params_with_no_prior')
M.params = M.params*1.01; %vary parameters
M_local.params = M_local.params*1.01; %vary parameters
else
M.params(estim_params_.param_vals(:,1)) = xparam1(estim_params_.nvx+estim_params_.ncx+estim_params_.nvn+estim_params_.ncn+1:end); %set parameters
M.params(estim_params_.param_vals(:,1)) = M.params(estim_params_.param_vals(:,1))*1.01; %vary parameters
M_local.params(estim_params_.param_vals(:,1)) = xparam1(estim_params_.nvx+estim_params_.ncx+estim_params_.nvn+estim_params_.ncn+1:end); %set parameters
M_local.params(estim_params_.param_vals(:,1)) = M_local.params(estim_params_.param_vals(:,1))*1.01; %vary parameters
end
if options_.diffuse_filter || options_.steadystate.nocheck
steadystate_check_flag = 0;
else
steadystate_check_flag = 1;
end
[tmp1, params] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M,options_,steadystate_check_flag);
change_flag=any(find(params-M.params));
[tmp1, params] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_local,options_,steadystate_check_flag);
change_flag=any(find(params-M_local.params));
if change_flag
skipline()
if any(isnan(params))
disp('After computing the steadystate, the following parameters are still NaN: '),
disp(char(M.param_names(isnan(params))))
disp(char(M_local.param_names(isnan(params))))
end
if any(find(params(~isnan(params))-M.params(~isnan(params))))
if any(find(params(~isnan(params))-M_local.params(~isnan(params))))
disp('The steadystate file changed the values for the following parameters: '),
disp(char(M.param_names(find(params(~isnan(params))-M.params(~isnan(params))))))
disp(char(M_local.param_names(find(params(~isnan(params))-M_local.params(~isnan(params))))))
end
disp('The derivatives of jacobian and steady-state will be computed numerically'),
disp('(re-set options_.analytic_derivation_mode= -2)'),
@ -414,21 +414,21 @@ else
end
%check steady state at initial parameters
M = M_;
M_local = M_;
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
if estim_params_.np
M.params(estim_params_.param_vals(:,1)) = xparam1(nvx+ncx+nvn+ncn+1:end);
M_local.params(estim_params_.param_vals(:,1)) = xparam1(nvx+ncx+nvn+ncn+1:end);
end
[oo_.steady_state, params,info] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M,options_,steadystate_check_flag);
[oo_.steady_state, params,info] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_local,options_,steadystate_check_flag);
if info(1)
fprintf('\ndynare_estimation_init:: The steady state at the initial parameters cannot be computed.\n')
if options_.debug
M.params=params;
plist = list_of_parameters_calibrated_as_NaN(M);
M_local.params=params;
plist = list_of_parameters_calibrated_as_NaN(M_local);
if ~isempty(plist)
message = ['dynare_estimation_init:: Some of the parameters are NaN (' ];
for i=1:length(plist)
@ -440,7 +440,7 @@ if info(1)
end
end
fprintf('%s\n',message)
plist = list_of_parameters_calibrated_as_Inf(M);
plist = list_of_parameters_calibrated_as_Inf(M_local);
if ~isempty(plist)
message = ['dynare_estimation_init:: Some of the parameters are Inf (' ];
for i=1:length(plist)
@ -526,9 +526,9 @@ end
if (options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter) || (options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter)
if isempty(options_.occbin.likelihood.IVF_shock_observable_mapping)
options_.occbin.likelihood.IVF_shock_observable_mapping=find(diag(M.Sigma_e)~=0);
options_.occbin.likelihood.IVF_shock_observable_mapping=find(diag(M_local.Sigma_e)~=0);
else
zero_var_shocks=find(diag(M.Sigma_e)==0);
zero_var_shocks=find(diag(M_local.Sigma_e)==0);
if any(ismember(options_.occbin.likelihood.IVF_shock_observable_mapping,zero_var_shocks))
error('IVF-filter: an observable is mapped to a zero variance shock.')
end

18
matlab/evaluate_static_model.m
View File

@ -1,6 +1,6 @@
function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,options)
function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M_,options_)
% function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,options)
% function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M_,options_)
% Evaluates the static model
%
% INPUTS
@ -8,8 +8,8 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt
% state
% exo_ss vector exogenous steady state
% params vector parameters
% M struct model structure
% options struct options
% M_ struct model structure
% options_ struct options
%
% OUTPUTS
% residuals vector residuals when ys is not
@ -38,18 +38,18 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
check1 = 0;
if options.bytecode
if options_.bytecode
if nargout<3
[residuals]= bytecode('evaluate', 'static', M, options, ys, ...
[residuals]= bytecode('evaluate', 'static', M_, options_, ys, ...
exo_ss, params, ys, 1);
else
[residuals, junk]= bytecode('evaluate', 'static', M, options, ys, ...
[residuals, junk]= bytecode('evaluate', 'static', M_, options_, ys, ...
exo_ss, params, ys, 1);
jacob = junk.g1;
end
else
[residuals, T_order, T] = feval([M.fname '.sparse.static_resid'], ys, exo_ss, params);
[residuals, T_order, T] = feval([M_.fname '.sparse.static_resid'], ys, exo_ss, params);
if nargout >= 3
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);
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);
end
end

296
matlab/evaluate_steady_state.m
View File

@ -1,13 +1,13 @@
function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M,options,steadystate_check_flag)
% function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M,options,steadystate_check_flag)
function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M_,options_,steadystate_check_flag)
% function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M_,options_,steadystate_check_flag)
% Computes the steady state
%
% INPUTS
% ys_init vector initial values used to compute the steady
% state
% exo_ss vector exogenous steady state (incl. deterministic exogenous)
% M struct model structure
% options struct options
% M_ struct model structure
% options_ struct options
% steadystate_check_flag boolean if true, check that the
% steadystate verifies the
% static model
@ -39,21 +39,21 @@ function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M,options,stead
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if options.solve_algo < 0 || options.solve_algo > 14
if options_.solve_algo < 0 || options_.solve_algo > 14
error('STEADY: solve_algo must be between 0 and 14')
end
if ~options.bytecode && ~options.block && options.solve_algo > 4 && ...
options.solve_algo < 9
error('STEADY: you can''t use solve_algo = {5,6,7,8} without block nor bytecode options')
if ~options_.bytecode && ~options_.block && options_.solve_algo > 4 && ...
options_.solve_algo < 9
error('STEADY: you can''t use solve_algo = {5,6,7,8} without block nor bytecode options_')
end
if ~options.bytecode && options.block && options.solve_algo == 5
if ~options_.bytecode && options_.block && options_.solve_algo == 5
error('STEADY: you can''t use solve_algo = 5 without bytecode option')
end
if isoctave && options.solve_algo == 11
error(['STEADY: you can''t use solve_algo = %u under Octave'],options.solve_algo)
if isoctave && options_.solve_algo == 11
error(['STEADY: you can''t use solve_algo = %u under Octave'],options_.solve_algo)
end
% To ensure that the z and zx matrices constructed by repmat and passed to bytecode
@ -68,36 +68,36 @@ end
info = 0;
check = 0;
steadystate_flag = options.steadystate_flag;
params = M.params;
steadystate_flag = options_.steadystate_flag;
params = M_.params;
if length(M.aux_vars) > 0 && ~steadystate_flag && M.set_auxiliary_variables
h_set_auxiliary_variables = str2func([M.fname '.set_auxiliary_variables']);
if length(M_.aux_vars) > 0 && ~steadystate_flag && M_.set_auxiliary_variables
h_set_auxiliary_variables = str2func([M_.fname '.set_auxiliary_variables']);
ys_init = h_set_auxiliary_variables(ys_init,exo_ss,params);
end
if options.ramsey_policy
if ~isfinite(M.params(strmatch('optimal_policy_discount_factor',M.param_names,'exact')))
if options_.ramsey_policy
if ~isfinite(M_.params(strmatch('optimal_policy_discount_factor',M_.param_names,'exact')))
fprintf('\nevaluate_steady_state: the planner_discount is NaN/Inf. That will cause problems.\n')
end
if steadystate_flag
% explicit steady state file
[ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M, ...
options,steadystate_check_flag);
[ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M_, ...
options_,steadystate_check_flag);
%test whether it solves model conditional on the instruments
if ~options.debug
resids = evaluate_static_model(ys,exo_ss,params,M,options);
if ~options_.debug
resids = evaluate_static_model(ys,exo_ss,params,M_,options_);
else
[resids, ~ , jacob]= evaluate_static_model(ys,exo_ss,params,M,options);
[resids, ~ , jacob]= evaluate_static_model(ys,exo_ss,params,M_,options_);
end
nan_indices=find(isnan(resids(M.ramsey_orig_endo_nbr+(1:M.ramsey_orig_eq_nbr))));
nan_indices=find(isnan(resids(M_.ramsey_orig_endo_nbr+(1:M_.ramsey_orig_eq_nbr))));
if ~isempty(nan_indices)
if options.debug
if options_.debug
fprintf('\nevaluate_steady_state: The steady state file computation for the Ramsey problem resulted in NaNs.\n')
fprintf('evaluate_steady_state: The steady state was computed conditional on the following initial instrument values: \n')
for ii = 1:size(options.instruments,1)
fprintf('\t %s \t %f \n',options.instruments{ii},ys_init(strmatch(options.instruments{ii},M.endo_names,'exact')))
for ii = 1:size(options_.instruments,1)
fprintf('\t %s \t %f \n',options_.instruments{ii},ys_init(strmatch(options_.instruments{ii},M_.endo_names,'exact')))
end
fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
fprintf('\t Equation(s): ')
@ -113,12 +113,12 @@ if options.ramsey_policy
return
end
if any(imag(ys(M.ramsey_orig_endo_nbr+(1:M.ramsey_orig_eq_nbr))))
if options.debug
if any(imag(ys(M_.ramsey_orig_endo_nbr+(1:M_.ramsey_orig_eq_nbr))))
if options_.debug
fprintf('\nevaluate_steady_state: The steady state file computation for the Ramsey problem resulted in complex numbers.\n')
fprintf('evaluate_steady_state: The steady state was computed conditional on the following initial instrument values: \n')
for ii = 1:size(options.instruments,1)
fprintf('\t %s \t %f \n',options.instruments{ii},ys_init(strmatch(options.instruments{ii},M.endo_names,'exact')))
for ii = 1:size(options_.instruments,1)
fprintf('\t %s \t %f \n',options_.instruments{ii},ys_init(strmatch(options_.instruments{ii},M_.endo_names,'exact')))
end
fprintf('evaluate_steady_state: If those initial values are not admissable, change them using an initval-block.\n')
skipline(2)
@ -128,17 +128,17 @@ if options.ramsey_policy
return
end
if max(abs(resids(M.ramsey_orig_endo_nbr+(1:M.ramsey_orig_eq_nbr)))) > options.solve_tolf %does it solve for all variables except for the Lagrange multipliers
if options.debug
if max(abs(resids(M_.ramsey_orig_endo_nbr+(1:M_.ramsey_orig_eq_nbr)))) > options_.solve_tolf %does it solve for all variables except for the Lagrange multipliers
if options_.debug
fprintf('\nevaluate_steady_state: The steady state file does not solve the steady state for the Ramsey problem.\n')
fprintf('evaluate_steady_state: Conditional on the following instrument values: \n')
for ii = 1:size(options.instruments,1)
fprintf('\t %s \t %f \n',options.instruments{ii},ys_init(strmatch(options.instruments{ii},M.endo_names,'exact')))
for ii = 1:size(options_.instruments,1)
fprintf('\t %s \t %f \n',options_.instruments{ii},ys_init(strmatch(options_.instruments{ii},M_.endo_names,'exact')))
end
fprintf('evaluate_steady_state: the following equations have non-zero residuals: \n')
for ii=M.ramsey_orig_endo_nbr+1:M.endo_nbr
if abs(resids(ii)) > options.solve_tolf
fprintf('\t Equation number %d: %f\n',ii-M.ramsey_orig_endo_nbr, resids(ii))
for ii=M_.ramsey_orig_endo_nbr+1:M_.endo_nbr
if abs(resids(ii)) > options_.solve_tolf
fprintf('\t Equation number %d: %f\n',ii-M_.ramsey_orig_endo_nbr, resids(ii))
end
end
skipline(2)
@ -148,31 +148,31 @@ if options.ramsey_policy
return
end
end
if options.debug
if options_.debug
if steadystate_flag
infrow=find(isinf(ys_init(1:M.orig_endo_nbr)));
infrow=find(isinf(ys_init(1:M_.orig_endo_nbr)));
else
infrow=find(isinf(ys_init));
end
if ~isempty(infrow)
fprintf('\nevaluate_steady_state: The initial values for the steady state of the following variables are Inf:\n');
for iter=1:length(infrow)
fprintf('%s\n',M.endo_names{infrow(iter)});
fprintf('%s\n',M_.endo_names{infrow(iter)});
end
end
if steadystate_flag
nanrow=find(isnan(ys_init(1:M.orig_endo_nbr)));
nanrow=find(isnan(ys_init(1:M_.orig_endo_nbr)));
else
nanrow=find(isnan(ys_init));
end
if ~isempty(nanrow)
fprintf('\nevaluate_steady_state: The initial values for the steady state of the following variables are NaN:\n');
for iter=1:length(nanrow)
fprintf('%s\n',M.endo_names{nanrow(iter)});
fprintf('%s\n',M_.endo_names{nanrow(iter)});
end
end
if steadystate_flag
nan_indices_mult=find(isnan(resids(1:M.ramsey_orig_endo_nbr)));
nan_indices_mult=find(isnan(resids(1:M_.ramsey_orig_endo_nbr)));
if any(nan_indices_mult)
fprintf('evaluate_steady_state: The steady state results NaN for auxiliary equation %u.\n',nan_indices_mult);
fprintf('evaluate_steady_state: This is often a sign of problems.\n');
@ -181,41 +181,41 @@ if options.ramsey_policy
if ~isempty(infrow)
fprintf('\nevaluate_steady_state: The Jacobian of the dynamic model contains Inf. The problem is associated with:\n\n')
display_problematic_vars_Jacobian(infrow,infcol,M,ys,'static','evaluate_steady_state: ')
display_problematic_vars_Jacobian(infrow,infcol,M_,ys,'static','evaluate_steady_state: ')
end
if ~isreal(jacob)
[imagrow,imagcol]=find(abs(imag(jacob))>1e-15);
fprintf('\nevaluate_steady_state: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n')
display_problematic_vars_Jacobian(imagrow,imagcol,M,ys,'static','evaluate_steady_state: ')
display_problematic_vars_Jacobian(imagrow,imagcol,M_,ys,'static','evaluate_steady_state: ')
end
[nanrow,nancol]=find(isnan(jacob));
if ~isempty(nanrow)
fprintf('\nevaluate_steady_state: The Jacobian of the dynamic model contains NaN. The problem is associated with:\n\n')
display_problematic_vars_Jacobian(nanrow,nancol,M,ys,'static','evaluate_steady_state: ')
display_problematic_vars_Jacobian(nanrow,nancol,M_,ys,'static','evaluate_steady_state: ')
end
end
end
%either if no steady state file or steady state file without problems
[ys,params,info] = dyn_ramsey_static(ys_init,exo_ss,M,options);
[ys,params,info] = dyn_ramsey_static(ys_init,exo_ss,M_,options_);
if info
return
end
%check whether steady state really solves the model
resids = evaluate_static_model(ys,exo_ss,params,M,options);
resids = evaluate_static_model(ys,exo_ss,params,M_,options_);
nan_indices_multiplier=find(isnan(resids(1:M.ramsey_orig_endo_nbr)));
nan_indices=find(isnan(resids(M.ramsey_orig_endo_nbr+1:end)));
nan_indices_multiplier=find(isnan(resids(1:M_.ramsey_orig_endo_nbr)));
nan_indices=find(isnan(resids(M_.ramsey_orig_endo_nbr+1:end)));
if ~isempty(nan_indices)
if options.debug
if options_.debug
fprintf('\nevaluate_steady_state: The steady state computation for the Ramsey problem resulted in NaNs.\n')
fprintf('evaluate_steady_state: The steady state computation resulted in the following instrument values: \n')
for i = 1:size(options.instruments,1)
fprintf('\t %s \t %f \n',options.instruments{i},ys(strmatch(options.instruments{i},M.endo_names,'exact')))
for i = 1:size(options_.instruments,1)
fprintf('\t %s \t %f \n',options_.instruments{i},ys(strmatch(options_.instruments{i},M_.endo_names,'exact')))
end
fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
fprintf('\t Equation(s): ')
@ -229,11 +229,11 @@ if options.ramsey_policy
end
if ~isempty(nan_indices_multiplier)
if options.debug
if options_.debug
fprintf('\nevaluate_steady_state: The steady state computation for the Ramsey problem resulted in NaNs in the auxiliary equations.\n')
fprintf('evaluate_steady_state: The steady state computation resulted in the following instrument values: \n')
for i = 1:size(options.instruments,1)
fprintf('\t %s \t %f \n',options.instruments{i},ys(strmatch(options.instruments{i},M.endo_names,'exact')))
for i = 1:size(options_.instruments,1)
fprintf('\t %s \t %f \n',options_.instruments{i},ys(strmatch(options_.instruments{i},M_.endo_names,'exact')))
end
fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
fprintf('\t Auxiliary equation(s): ')
@ -246,22 +246,22 @@ if options.ramsey_policy
return
end
if max(abs(resids)) > options.solve_tolf %does it solve for all variables including the auxiliary ones
if options.debug
if max(abs(resids)) > options_.solve_tolf %does it solve for all variables including the auxiliary ones
if options_.debug
fprintf('\nevaluate_steady_state: The steady state for the Ramsey problem could not be computed.\n')
fprintf('evaluate_steady_state: The steady state computation stopped with the following instrument values:: \n')
for i = 1:size(options.instruments,1)
fprintf('\t %s \t %f \n',options.instruments{i},ys(strmatch(options.instruments{i},M.endo_names,'exact')))
for i = 1:size(options_.instruments,1)
fprintf('\t %s \t %f \n',options_.instruments{i},ys(strmatch(options_.instruments{i},M_.endo_names,'exact')))
end
fprintf('evaluate_steady_state: The following equations have non-zero residuals: \n')
for ii=1:M.ramsey_orig_endo_nbr
if abs(resids(ii)) > options.solve_tolf/100
for ii=1:M_.ramsey_orig_endo_nbr
if abs(resids(ii)) > options_.solve_tolf/100
fprintf('\t Auxiliary Ramsey equation number %d: %f\n',ii, resids(ii))
end
end
for ii=M.ramsey_orig_endo_nbr+1:M.endo_nbr
if abs(resids(ii)) > options.solve_tolf/100
fprintf('\t Equation number %d: %f\n',ii-M.ramsey_orig_endo_nbr, resids(ii))
for ii=M_.ramsey_orig_endo_nbr+1:M_.endo_nbr
if abs(resids(ii)) > options_.solve_tolf/100
fprintf('\t Equation number %d: %f\n',ii-M_.ramsey_orig_endo_nbr, resids(ii))
end
end
skipline(2)
@ -272,46 +272,46 @@ if options.ramsey_policy
end
elseif steadystate_flag
% explicit steady state file
[ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M, options,steadystate_check_flag);
[ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M_, options_,steadystate_check_flag);
if size(ys,2)>size(ys,1)
error('STEADY: steady_state-file must return a column vector, not a row vector.')
end
if info(1)
return
end
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
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
% non linear model
if ismember(options.solve_algo,[10,11])
[lb,ub,eq_index] = get_complementarity_conditions(M,options.ramsey_policy);
if options.solve_algo == 10
options.lmmcp.lb = lb;
options.lmmcp.ub = ub;
elseif options.solve_algo == 11
options.mcppath.lb = lb;
options.mcppath.ub = ub;
if ismember(options_.solve_algo,[10,11])
[lb,ub,eq_index] = get_complementarity_conditions(M_,options_.ramsey_policy);
if options_.solve_algo == 10
options_.lmmcp.lb = lb;
options_.lmmcp.ub = ub;
elseif options_.solve_algo == 11
options_.mcppath.lb = lb;
options_.mcppath.ub = ub;
end
[ys,check,fvec, fjac, errorcode] = dynare_solve(@static_mcp_problem,...
ys_init,...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
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, eq_index);
options_.steady.maxit, options_.solve_tolf, options_.solve_tolx, ...
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, eq_index);
else
[ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
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);
options_.steady.maxit, options_.solve_tolf, options_.solve_tolx, ...
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
if check && options.debug
if check && options_.debug
dprintf('Nonlinear solver routine returned errorcode=%i.', errorcode)
skipline()
[infrow,infcol]=find(isinf(fjac) | isnan(fjac));
if ~isempty(infrow)
fprintf('\nSTEADY: The Jacobian at the initial values contains Inf or NaN. The problem arises from: \n')
display_problematic_vars_Jacobian(infrow,infcol,M,ys_init,'static','STEADY: ')
display_problematic_vars_Jacobian(infrow,infcol,M_,ys_init,'static','STEADY: ')
end
problematic_equation = find(~isfinite(fvec));
if ~isempty(problematic_equation)
@ -325,7 +325,7 @@ elseif ~options.bytecode && ~options.block
else
% linear model
[fvec, T_order, T] = static_resid(ys_init, exo_ss, params);
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);
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);
ii = find(~isfinite(fvec));
if ~isempty(ii)
@ -337,7 +337,7 @@ elseif ~options.bytecode && ~options.block
disp('Check whether your model is truly linear. Put "resid(1);" before "steady;" to see the problematic equations.')
elseif isempty(ii) && max(abs(fvec)) > 1e-12
ys = ys_init-jacob\fvec;
resid = evaluate_static_model(ys,exo_ss,params,M,options);
resid = evaluate_static_model(ys,exo_ss,params,M_,options_);
if max(abs(resid)) > 1e-6
check=1;
fprintf('STEADY: No steady state for your model could be found\n')
@ -346,39 +346,39 @@ elseif ~options.bytecode && ~options.block
else
ys = ys_init;
end
if options.debug
if options_.debug
if any(any(isinf(jacob) | isnan(jacob)))
[infrow,infcol]=find(isinf(jacob) | isnan(jacob));
fprintf('\nSTEADY: The Jacobian contains Inf or NaN. The problem arises from: \n\n')
for ii=1:length(infrow)
fprintf('STEADY: Derivative of Equation %d with respect to Variable %s (initial value of %s: %g) \n',infrow(ii),M.endo_names{infcol(ii),:},M.endo_names{infcol(ii),:},ys_init(infcol(ii)))
fprintf('STEADY: Derivative of Equation %d with respect to Variable %s (initial value of %s: %g) \n',infrow(ii),M_.endo_names{infcol(ii),:},M_.endo_names{infcol(ii),:},ys_init(infcol(ii)))
end
fprintf('Check whether your model is truly linear. Put "resid(1);" before "steady;" to see the problematic equations.\n')
end
end
end
elseif ~options.bytecode && options.block
elseif ~options_.bytecode && options_.block
ys = ys_init;
T = NaN(M.block_structure_stat.tmp_nbr, 1);
for b = 1:length(M.block_structure_stat.block)
fh_static = str2func(sprintf('%s.sparse.block.static_%d', M.fname, b));
if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
M.block_structure_stat.block(b).Simulation_Type ~= 2
mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end);
if options.solve_algo <= 4 || options.solve_algo >= 9
T = NaN(M_.block_structure_stat.tmp_nbr, 1);
for b = 1:length(M_.block_structure_stat.block)
fh_static = str2func(sprintf('%s.sparse.block.static_%d', M_.fname, b));
if M_.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
M_.block_structure_stat.block(b).Simulation_Type ~= 2
mfs_idx = M_.block_structure_stat.block(b).variable(end-M_.block_structure_stat.block(b).mfs+1:end);
if options_.solve_algo <= 4 || options_.solve_algo >= 9
[ys(mfs_idx), errorflag] = dynare_solve(@block_mfs_steadystate, ys(mfs_idx), ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, fh_static, b, ys, exo_ss, params, T, M);
options_.steady.maxit, options_.solve_tolf, options_.solve_tolx, ...
options_, fh_static, b, ys, exo_ss, params, T, M_);
if errorflag
check = 1;
break
end
else
nze = length(M.block_structure_stat.block(b).g1_sparse_rowval);
nze = length(M_.block_structure_stat.block(b).g1_sparse_rowval);
[ys, T, success] = solve_one_boundary(fh_static, ys, exo_ss, ...
params, [], T, mfs_idx, nze, 1, false, b, 0, options.steady.maxit, ...
options.solve_tolf, ...
0, options.solve_algo, true, false, false, M, options);
params, [], T, mfs_idx, nze, 1, false, b, 0, options_.steady.maxit, ...
options_.solve_tolf, ...
0, options_.solve_algo, true, false, false, M_, options_);
if ~success
check = 1;
break
@ -387,36 +387,36 @@ elseif ~options.bytecode && options.block
end
% Compute endogenous if the block is of type evaluate forward/backward or if there are recursive variables in a solve block.
% Also update the temporary terms vector (needed for the dynare_solve case)
[ys, T] = fh_static(ys, exo_ss, params, M.block_structure_stat.block(b).g1_sparse_rowval, ...
M.block_structure_stat.block(b).g1_sparse_colval, ...
M.block_structure_stat.block(b).g1_sparse_colptr, T);
[ys, T] = fh_static(ys, exo_ss, params, M_.block_structure_stat.block(b).g1_sparse_rowval, ...
M_.block_structure_stat.block(b).g1_sparse_colval, ...
M_.block_structure_stat.block(b).g1_sparse_colptr, T);
end
elseif options.bytecode
if options.solve_algo >= 5 && options.solve_algo <= 8
elseif options_.bytecode
if options_.solve_algo >= 5 && options_.solve_algo <= 8
try
if options.block
ys = bytecode('static', 'block_decomposed', M, options, ys_init, exo_ss, params);
if options_.block
ys = bytecode('static', 'block_decomposed', M_, options_, ys_init, exo_ss, params);
else
ys = bytecode('static', M, options, ys_init, exo_ss, params);
ys = bytecode('static', M_, options_, ys_init, exo_ss, params);
end
catch ME
if options.verbosity >= 1
if options_.verbosity >= 1
disp(ME.message);
end
ys = ys_init;
check = 1;
end
elseif options.block
elseif options_.block
ys = ys_init;
T = NaN(M.block_structure_stat.tmp_nbr, 1);
for b = 1:length(M.block_structure_stat.block)
if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
M.block_structure_stat.block(b).Simulation_Type ~= 2
mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end);
T = NaN(M_.block_structure_stat.tmp_nbr, 1);
for b = 1:length(M_.block_structure_stat.block)
if M_.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
M_.block_structure_stat.block(b).Simulation_Type ~= 2
mfs_idx = M_.block_structure_stat.block(b).variable(end-M_.block_structure_stat.block(b).mfs+1:end);
[ys(mfs_idx), errorflag] = dynare_solve(@block_bytecode_mfs_steadystate, ...
ys(mfs_idx), options.steady.maxit, ...
options.solve_tolf, options.solve_tolx, ...
options, b, ys, exo_ss, params, T, M, options);
ys(mfs_idx), options_.steady.maxit, ...
options_.solve_tolf, options_.solve_tolx, ...
options_, b, ys, exo_ss, params, T, M_, options_);
if errorflag
check = 1;
break
@ -425,10 +425,10 @@ elseif options.bytecode
% Compute endogenous if the block is of type evaluate forward/backward or if there are recursive variables in a solve block.
% Also update the temporary terms vector (needed for the dynare_solve case)
try
[~, ~, ys, T] = bytecode(M, options, ys, exo_ss, params, ys, 1, ys, T, 'evaluate', 'static', ...
[~, ~, ys, T] = bytecode(M_, options_, ys, exo_ss, params, ys, 1, ys, T, 'evaluate', 'static', ...
'block_decomposed', ['block = ' int2str(b)]);
catch ME
if options.verbosity >= 1
if options_.verbosity >= 1
disp(ME.message);
end
check = 1;
@ -437,20 +437,20 @@ elseif options.bytecode
end
else
[ys, check] = dynare_solve(@bytecode_steadystate, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params, M, options);
options_.steady.maxit, options_.solve_tolf, options_.solve_tolx, ...
options_, exo_ss, params, M_, options_);
end
end
if check
info(1)= 20;
%make sure ys contains auxiliary variables in case of problem with dynare_solve
if length(M.aux_vars) > 0 && ~steadystate_flag
if M.set_auxiliary_variables
if length(M_.aux_vars) > 0 && ~steadystate_flag
if M_.set_auxiliary_variables
ys = h_set_auxiliary_variables(ys,exo_ss,params);
end
end
resid = evaluate_static_model(ys,exo_ss,params,M,options);
resid = evaluate_static_model(ys,exo_ss,params,M_,options_);
info(2) = resid'*resid ;
if isnan(info(2))
info(1)=22;
@ -459,18 +459,18 @@ if check
end
% If some equations are tagged [static] or [dynamic], verify consistency
if M.static_and_dynamic_models_differ
if M_.static_and_dynamic_models_differ
% Evaluate residual of *dynamic* model using the steady state
% computed on the *static* one
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', M, options, z, zx, params, ys, 1);
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', M_, options_, z, zx, params, ys, 1);
else
r = feval([M.fname '.sparse.dynamic_resid'], repmat(ys, 3, 1), exo_ss, params, ys);
r = feval([M_.fname '.sparse.dynamic_resid'], repmat(ys, 3, 1), exo_ss, params, ys);
end
% Fail if residual greater than tolerance
if max(abs(r)) > options.solve_tolf
if max(abs(r)) > options_.solve_tolf
info(1) = 25;
return
end
@ -505,21 +505,21 @@ j = fh_static_g1(y, x, params, sparse_rowval, sparse_colval, sparse_colptr, T_or
resids = r(eq_index);
jac = j(eq_index,1:nvar);
function [r, g1] = block_mfs_steadystate(y, fh_static, b, y_all, exo, params, T, M)
function [r, g1] = block_mfs_steadystate(y, fh_static, b, y_all, exo, params, T, M_)
% Wrapper around the static files, for block without bytecode
mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end);
mfs_idx = M_.block_structure_stat.block(b).variable(end-M_.block_structure_stat.block(b).mfs+1:end);
y_all(mfs_idx) = y;
[~,~,r,g1] = fh_static(y_all, exo, params, M.block_structure_stat.block(b).g1_sparse_rowval, ...
M.block_structure_stat.block(b).g1_sparse_colval, ...
M.block_structure_stat.block(b).g1_sparse_colptr, T);
[~,~,r,g1] = fh_static(y_all, exo, params, M_.block_structure_stat.block(b).g1_sparse_rowval, ...
M_.block_structure_stat.block(b).g1_sparse_colval, ...
M_.block_structure_stat.block(b).g1_sparse_colptr, T);
function [r, g1] = bytecode_steadystate(y, exo, params, M, options)
function [r, g1] = bytecode_steadystate(y, exo, params, M_, options_)
% Wrapper around the static file, for bytecode (without block)
[r, g1] = bytecode(M, options, y, exo, params, y, 1, exo, 'evaluate', 'static');
[r, g1] = bytecode(M_, options_, y, exo, params, y, 1, exo, 'evaluate', 'static');
function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, T, M, options)
function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, T, M_, options_)
% Wrapper around the static files, for bytecode with block
mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end);
mfs_idx = M_.block_structure_stat.block(b).variable(end-M_.block_structure_stat.block(b).mfs+1:end);
y_all(mfs_idx) = y;
[r, g1] = bytecode(M, options, y_all, exo, params, y_all, 1, y_all, T, 'evaluate', 'static', 'block_decomposed', ['block = ' int2str(b) ]);
g1 = g1(:,end-M.block_structure_stat.block(b).mfs+1:end); % Make Jacobian square if mfs>0
[r, g1] = bytecode(M_, options_, y_all, exo, params, y_all, 1, y_all, T, 'evaluate', 'static', 'block_decomposed', ['block = ' int2str(b) ]);
g1 = g1(:,end-M_.block_structure_stat.block(b).mfs+1:end); % Make Jacobian square if mfs>0

32
matlab/evaluate_steady_state_file.m
View File

@ -1,13 +1,13 @@
function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options,steady_state_checkflag)
% function [ys,params1,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options,steady_state_checkflag)
function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M_,options_,steady_state_checkflag)
% function [ys,params1,info] = evaluate_steady_state_file(ys_init,exo_ss,M_,options_,steady_state_checkflag)
% Evaluates steady state files
%
% INPUTS
% ys_init vector initial values used to compute the steady
% state
% exo_ss vector exogenous steady state (incl. deterministic exogenous)
% M struct model parameters
% options struct options
% M_ struct model parameters
% options_ struct options
% steady_state_checkflag boolean indicator whether to check steady state returned
% OUTPUTS
% ys vector steady state
@ -36,17 +36,17 @@ function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
params = M.params;
params = M_.params;
info = 0;
fname = M.fname;
fname = M_.fname;
if options.steadystate_flag == 1
if options_.steadystate_flag == 1
% old format
assignin('base','tmp_00_',params);
evalin('base','M_.params=tmp_00_; clear(''tmp_00_'')');
h_steadystate = str2func([fname '_steadystate']);
[ys,params1,check] = h_steadystate(ys_init, exo_ss,M,options);
[ys,params1,check] = h_steadystate(ys_init, exo_ss,M_,options_);
else % steadystate_flag == 2
% new format
h_steadystate = str2func([fname '.steadystate']);
@ -59,19 +59,19 @@ if check
return
end
if M.param_nbr > 0
if M_.param_nbr > 0
updated_params_flag = max(abs(params1-params)) > 1e-12 ...
|| ~isequal(isnan(params1),isnan(params)); %checks whether numbers or NaN changed
else
updated_params_flag = 0;
end
h_set_auxiliary_variables = str2func([M.fname '.set_auxiliary_variables']);
h_set_auxiliary_variables = str2func([M_.fname '.set_auxiliary_variables']);
if isnan(updated_params_flag) || (updated_params_flag && any(isnan(params(~isnan(params))-params1(~isnan(params))))) %checks if new NaNs were added
info(1) = 24;
info(2) = NaN;
if M.set_auxiliary_variables
if M_.set_auxiliary_variables
ys = h_set_auxiliary_variables(ys,exo_ss,params);
end
return
@ -80,7 +80,7 @@ end
if updated_params_flag && ~isreal(params1)
info(1) = 23;
info(2) = sum(imag(params).^2,'omitnan');
if M.set_auxiliary_variables
if M_.set_auxiliary_variables
ys = h_set_auxiliary_variables(ys,exo_ss,params);
end
return
@ -91,21 +91,21 @@ if updated_params_flag
end
% adding values for auxiliary variables
if ~isempty(M.aux_vars) && ~options.ramsey_policy
if M.set_auxiliary_variables
if ~isempty(M_.aux_vars) && ~options_.ramsey_policy
if M_.set_auxiliary_variables
ys = h_set_auxiliary_variables(ys,exo_ss,params);
end
end
if steady_state_checkflag
% Check whether the steady state obtained from the _steadystate file is a steady state.
[residuals, check] = evaluate_static_model(ys, exo_ss, params, M, options);
[residuals, check] = evaluate_static_model(ys, exo_ss, params, M_, options_);
if check
info(1) = 19;
info(2) = check; % to be improved
return
end
if max(abs(residuals)) > options.solve_tolf
if max(abs(residuals)) > options_.solve_tolf
info(1) = 19;
info(2) = residuals'*residuals;
return

46
matlab/expand_group.m
View File

@ -32,16 +32,16 @@ if nargin<4
end
filename = get(gcf,'filename');
[filepath, name, ext]=fileparts(filename);
M = evalin('base','M_');
oo = evalin('base','oo_');
options = evalin('base','options_');
M_ = evalin('base','M_');
oo_ = evalin('base','oo_');
options_ = evalin('base','options_');
mydata=get(findobj(gcf,'tag',['group' int2str(ic)]),'userdata');
if isfield(mydata,'shock_decomp')
options.shock_decomp=mydata.shock_decomp;
options_.shock_decomp=mydata.shock_decomp;
end
options.plot_shock_decomp=mydata.plot_shock_decomp;
options.first_obs=mydata.first_obs;
options.nobs=mydata.nobs;
options_.plot_shock_decomp=mydata.plot_shock_decomp;
options_.first_obs=mydata.first_obs;
options_.nobs=mydata.nobs;
% define expanded group
label = mydata.shock_group.label;
label = strrep(label,' ','_');
@ -50,32 +50,32 @@ label = strrep(label,'(','');
label = strrep(label,')','');
label = strrep(label,'.','');
shocks = mydata.shock_group.shocks;
options.plot_shock_decomp.fig_name = [mydata.fig_name '. Expand'];
options.plot_shock_decomp.use_shock_groups = label; %[use_shock_groups_old int2str(ic)];
options.plot_shock_decomp.filepath = filepath;
options_.plot_shock_decomp.fig_name = [mydata.fig_name '. Expand'];
options_.plot_shock_decomp.use_shock_groups = label; %[use_shock_groups_old int2str(ic)];
options_.plot_shock_decomp.filepath = filepath;
for j=1:length(shocks)
M.shock_groups.(options.plot_shock_decomp.use_shock_groups).(['group' int2str(j)]).label=shocks{j};
M.shock_groups.(options.plot_shock_decomp.use_shock_groups).(['group' int2str(j)]).shocks=shocks(j);
M_.shock_groups.(options_.plot_shock_decomp.use_shock_groups).(['group' int2str(j)]).label=shocks{j};
M_.shock_groups.(options_.plot_shock_decomp.use_shock_groups).(['group' int2str(j)]).shocks=shocks(j);
end
M.exo_names = mydata.exo_names;
M_.exo_names = mydata.exo_names;
options.plot_shock_decomp.interactive=0;
options.plot_shock_decomp.expand=1;
options.plot_shock_decomp.nodisplay=0;
options_.plot_shock_decomp.interactive=0;
options_.plot_shock_decomp.expand=1;
options_.plot_shock_decomp.nodisplay=0;
if no_graph
options.no_graph.plot_shock_decomposition=1;
options.plot_shock_decomp.write_xls=1;
options_.no_graph.plot_shock_decomposition=1;
options_.plot_shock_decomp.write_xls=1;
else
options.plot_shock_decomp.write_xls=0;
options_.plot_shock_decomp.write_xls=0;
end
%% set optimal colormap
func = @(x) colorspace('RGB->Lab',x);
MAP = distinguishable_colors(length(shocks)+1,'w',func);
MAP(end,:) = [0.7 0.7 0.7];
options.plot_shock_decomp.colormap = MAP;
options_.plot_shock_decomp.colormap = MAP;
options.plot_shock_decomp.endo_names = mydata.endo_names;
options.plot_shock_decomp.endo_names_tex = mydata.endo_names_tex;
options_.plot_shock_decomp.endo_names = mydata.endo_names;
options_.plot_shock_decomp.endo_names_tex = mydata.endo_names_tex;
plot_shock_decomposition(M,oo,options,{var_list_});
plot_shock_decomposition(M_,oo_,options_,{var_list_});

36
matlab/get_identification_jacobians.m
View File

@ -1,12 +1,12 @@
function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% previously getJJ.m in Dynare 4.5
% Sets up the Jacobians needed for identification analysis
% =========================================================================
% INPUTS
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% options: [structure] storing the options
% M_: [structure] storing the model information
% options_: [structure] storing the options
% options_ident: [structure] storing the options for identification
% indpmodel: [modparam_nbr by 1] index of estimated parameters in M_.params;
% corresponds to model parameters (no stderr and no corr)
@ -126,19 +126,19 @@ useautocorr = options_ident.useautocorr;
grid_nbr = options_ident.grid_nbr;
kronflag = options_ident.analytic_derivation_mode;
% get fields from M
endo_nbr = M.endo_nbr;
exo_nbr = M.exo_nbr;
fname = M.fname;
lead_lag_incidence = M.lead_lag_incidence;
nspred = M.nspred;
nstatic = M.nstatic;
params = M.params;
Sigma_e = M.Sigma_e;
% get fields from M_
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
fname = M_.fname;
lead_lag_incidence = M_.lead_lag_incidence;
nspred = M_.nspred;
nstatic = M_.nstatic;
params = M_.params;
Sigma_e = M_.Sigma_e;
stderr_e = sqrt(diag(Sigma_e));
% set all selected values: stderr and corr come first, then model parameters
xparam1 = get_all_parameters(estim_params, M); %try using estimated_params block
xparam1 = get_all_parameters(estim_params, M_); %try using estimated_params block
if isempty(xparam1)
%if there is no estimated_params block, consider all stderr and all model parameters, but no corr parameters
xparam1 = [stderr_e', params'];
@ -153,7 +153,7 @@ obs_nbr = length(indvobs);
d2flag = 0; % do not compute second parameter derivatives
% Get Jacobians (wrt selected params) of steady state, dynamic model derivatives and perturbation solution matrices for all endogenous variables
dr.derivs = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag);
dr.derivs = get_perturbation_params_derivs(M_, options_, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag);
[I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files
yy0 = dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order
@ -230,7 +230,7 @@ elseif order == 3
end
% Get (pruned) state space representation:
pruned = pruned_state_space_system(M, options, dr, indvobs, nlags, useautocorr, 1);
pruned = pruned_state_space_system(M_, options_, dr, indvobs, nlags, useautocorr, 1);
MEAN = pruned.E_y;
dMEAN = pruned.dE_y;
%storage for Jacobians used in dsge_likelihood.m for analytical Gradient and Hession of likelihood (only at order=1)
@ -258,7 +258,7 @@ if ~no_identification_moments
if kronflag == -1
%numerical derivative of autocovariogram
dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M, options, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1]
dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1]
dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables
else
dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr);
@ -315,7 +315,7 @@ if ~no_identification_spectrum
IA = eye(size(pruned.A,1));
if kronflag == -1
%numerical derivative of spectral density
dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M, options, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2]
dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2]
kk = 0;
for ig = 1:length(freqs)
kk = kk+1;

124
matlab/get_perturbation_params_derivs.m
View File

@ -1,5 +1,5 @@
function DERIVS = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag)
% DERIVS = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag)
function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag)
% DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag)
% previously getH.m in dynare 4.5
% -------------------------------------------------------------------------
% Computes derivatives (with respect to parameters) of
@ -16,14 +16,14 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, dr, e
%
% =========================================================================
% INPUTS
% M: [structure] storing the model information
% options: [structure] storing the options
% estim_params: [structure] storing the estimation information
% M_: [structure] storing the model information
% options_: [structure] storing the options
% estim_params_: [structure] storing the estimation information
% dr [structure] Reduced form model.
% endo_steady_state [vector] steady state value for endogenous variables
% exo_steady_state [vector] steady state value for exogenous variables
% exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M.params;
% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M_.params;
% corresponds to model parameters (no stderr and no corr)
% in estimated_params block
% indpstderr: [stderrparam_nbr by 1] index of selected (estimated) standard errors,
@ -126,40 +126,40 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, dr, e
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
% Get fields from M
Correlation_matrix = M.Correlation_matrix;
dname = M.dname;
dynamic_tmp_nbr = M.dynamic_tmp_nbr;
endo_nbr = M.endo_nbr;
exo_nbr = M.exo_nbr;
exo_det_nbr = M.exo_det_nbr;
fname = M.fname;
lead_lag_incidence = M.lead_lag_incidence;
nfwrd = M.nfwrd;
npred = M.npred;
nspred = M.nspred;
nstatic = M.nstatic;
params = M.params;
param_nbr = M.param_nbr;
Sigma_e = M.Sigma_e;
% Get fields from M_
Correlation_matrix = M_.Correlation_matrix;
dname = M_.dname;
dynamic_tmp_nbr = M_.dynamic_tmp_nbr;
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
exo_det_nbr = M_.exo_det_nbr;
fname = M_.fname;
lead_lag_incidence = M_.lead_lag_incidence;
nfwrd = M_.nfwrd;
npred = M_.npred;
nspred = M_.nspred;
nstatic = M_.nstatic;
params = M_.params;
param_nbr = M_.param_nbr;
Sigma_e = M_.Sigma_e;
% Get fields from options
analytic_derivation_mode = options.analytic_derivation_mode;
% Get fields from options_
analytic_derivation_mode = options_.analytic_derivation_mode;
% analytic_derivation_mode: select method to compute Jacobians, default is 0
% * 0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012)
% * 1: kronecker products method to compute analytical derivatives as in Iskrev (2010), only for order=1
% * -1: numerical two-sided finite difference method to compute numerical derivatives of all output arguments using function get_perturbation_params_derivs_numerical_objective.m
% * -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the other output arguments are computed analytically as in kronflag=0
gstep = options.gstep;
order = options.order;
if isempty(options.qz_criterium)
gstep = options_.gstep;
order = options_.order;
if isempty(options_.qz_criterium)
% set default value for qz_criterium: if there are no unit roots one can use 1.0
% If they are possible, you may have have multiple unit roots and the accuracy
% decreases when computing the eigenvalues in lyapunov_symm. Hence, we normally use 1+1e-6
options = select_qz_criterium_value(options);
options_ = select_qz_criterium_value(options_);
end
qz_criterium = options.qz_criterium;
threads_BC = options.threads.kronecker.sparse_hessian_times_B_kronecker_C;
qz_criterium = options_.qz_criterium;
threads_BC = options_.threads.kronecker.sparse_hessian_times_B_kronecker_C;
% Get fields from dr
ghx = dr.ghx;
@ -187,7 +187,7 @@ if exo_det_nbr > 0
end
if order > 1 && analytic_derivation_mode == 1
%analytic derivatives using Kronecker products is implemented only at first-order, at higher order we reset to analytic derivatives with sylvester equations
%options.analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n');
%options_.analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n');
analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n');
end
@ -227,35 +227,35 @@ else %normal models
end
if analytic_derivation_mode < 0
%Create auxiliary estim_params blocks if not available for numerical derivatives, estim_params_model contains only model parameters
%Create auxiliary estim_params_ blocks if not available for numerical derivatives, estim_params_model contains only model parameters
estim_params_model.np = length(indpmodel);
estim_params_model.param_vals(:,1) = indpmodel;
estim_params_model.nvx = 0; estim_params_model.ncx = 0; estim_params_model.nvn = 0; estim_params_model.ncn = 0;
modparam1 = get_all_parameters(estim_params_model, M); %get all selected model parameters
if ~isempty(indpstderr) && isempty(estim_params.var_exo) %if there are stderr parameters but no estimated_params_block
modparam1 = get_all_parameters(estim_params_model, M_); %get all selected model parameters
if ~isempty(indpstderr) && isempty(estim_params_.var_exo) %if there are stderr parameters but no estimated_params_block
%provide temporary necessary information for stderr parameters
estim_params.nvx = length(indpstderr);
estim_params.var_exo(:,1) = indpstderr;
estim_params_.nvx = length(indpstderr);
estim_params_.var_exo(:,1) = indpstderr;
end
if ~isempty(indpcorr) && isempty(estim_params.corrx) %if there are corr parameters but no estimated_params_block
if ~isempty(indpcorr) && isempty(estim_params_.corrx) %if there are corr parameters but no estimated_params_block
%provide temporary necessary information for stderr parameters
estim_params.ncx = size(indpcorr,1);
estim_params.corrx(:,1:2) = indpcorr;
estim_params_.ncx = size(indpcorr,1);
estim_params_.corrx(:,1:2) = indpcorr;
end
if ~isfield(estim_params,'nvn') %stderr of measurement errors not yet
estim_params.nvn = 0;
estim_params.var_endo = [];
if ~isfield(estim_params_,'nvn') %stderr of measurement errors not yet
estim_params_.nvn = 0;
estim_params_.var_endo = [];
end
if ~isfield(estim_params,'ncn') %corr of measurement errors not yet
estim_params.ncn = 0;
estim_params.corrn = [];
if ~isfield(estim_params_,'ncn') %corr of measurement errors not yet
estim_params_.ncn = 0;
estim_params_.corrn = [];
end
if ~isempty(indpmodel) && isempty(estim_params.param_vals) %if there are model parameters but no estimated_params_block
if ~isempty(indpmodel) && isempty(estim_params_.param_vals) %if there are model parameters but no estimated_params_block
%provide temporary necessary information for model parameters
estim_params.np = length(indpmodel);
estim_params.param_vals(:,1) = indpmodel;
estim_params_.np = length(indpmodel);
estim_params_.param_vals(:,1) = indpmodel;
end
xparam1 = get_all_parameters(estim_params, M); %get all selected stderr, corr, and model parameters in estimated_params block, stderr and corr come first, then model parameters
xparam1 = get_all_parameters(estim_params_, M_); %get all selected stderr, corr, and model parameters in estimated_params block, stderr and corr come first, then model parameters
end
if d2flag
modparam_nbr2 = modparam_nbr*(modparam_nbr+1)/2; %number of unique entries of selected model parameters only in second-order derivative matrix
@ -295,7 +295,7 @@ if analytic_derivation_mode == -1
% - perturbation solution matrices: dghx, dghu, dghxx, dghxu, dghuu, dghs2, dghxxx, dghxxu, dghxuu, dghuuu, dghxss, dghuss
%Parameter Jacobian of covariance matrix and solution matrices (wrt selected stderr, corr and model paramters)
dSig_gh = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
dSig_gh = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Sigma_e = (1:exo_nbr^2);
ind_ghx = ind_Sigma_e(end) + (1:endo_nbr*nspred);
ind_ghu = ind_ghx(end) + (1:endo_nbr*exo_nbr);
@ -348,16 +348,16 @@ if analytic_derivation_mode == -1
end
%Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only)
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Yss = 1:endo_nbr;
if options.discretionary_policy || options.ramsey_policy
ind_g1 = ind_Yss(end) + (1:M.eq_nbr*yy0ex0_nbr);
if options_.discretionary_policy || options_.ramsey_policy
ind_g1 = ind_Yss(end) + (1:M_.eq_nbr*yy0ex0_nbr);
else
ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr);
end
DERIVS.dYss = dYss_g(ind_Yss, :); %in tensor notation, wrt selected model parameters only
if options.discretionary_policy || options.ramsey_policy
DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[M.eq_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only
if options_.discretionary_policy || options_.ramsey_policy
DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[M_.eq_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only
else
DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[endo_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only
end
@ -373,9 +373,9 @@ if analytic_derivation_mode == -1
if d2flag
% Hessian (wrt paramters) of steady state and first-order solution matrices ghx and Om
% note that hessian_sparse.m (contrary to hessian.m) does not take symmetry into account, but focuses already on unique values
options.order = 1; %make sure only first order
d2Yss_KalmanA_Om = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
options.order = order; %make sure to set back
options_.order = 1; %make sure only first order
d2Yss_KalmanA_Om = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
options_.order = order; %make sure to set back
ind_KalmanA = ind_Yss(end) + (1:endo_nbr^2);
DERIVS.d2KalmanA = d2Yss_KalmanA_Om(ind_KalmanA, indp2tottot2); %only unique elements
DERIVS.d2Om = d2Yss_KalmanA_Om(ind_KalmanA(end)+1:end , indp2tottot2); %only unique elements
@ -404,9 +404,9 @@ if analytic_derivation_mode == -2
if d2flag
% computation of d2Yss and d2g1
% note that hessian_sparse does not take symmetry into account, i.e. compare hessian_sparse.m to hessian.m, but focuses already on unique values, which are duplicated below
options.order = 1; %d2flag requires only first order
d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % d2flag requires only first-order
options.order = order; %make sure to set back the order
options_.order = 1; %d2flag requires only first order
d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % d2flag requires only first-order
options_.order = order; %make sure to set back the order
d2Yss = reshape(full(d2Yss_g1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation
for j=1:endo_nbr
d2Yss(j,:,:) = dyn_unvech(dyn_vech(d2Yss(j,:,:))); %add duplicate values to full hessian
@ -431,7 +431,7 @@ if analytic_derivation_mode == -2
end
%Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only)
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Yss = 1:endo_nbr;
ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr);
dYss = dYss_g(ind_Yss,:); %in tensor notation, wrt selected model parameters only

46
matlab/get_perturbation_params_derivs_numerical_objective.m
View File

@ -1,5 +1,5 @@
function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, options, dr, steady_state, exo_steady_state, exo_det_steady_state)
function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M_, options_, dr, steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Objective function to compute numerically the Jacobians used for get_perturbation_params_derivs
% =========================================================================
@ -8,8 +8,8 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% stderr parameters come first, corr parameters second, model parameters third
% outputflag: [string] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% options: [structure] storing the options
% M_: [structure] storing the model information
% options_: [structure] storing the options
% dr [structure] Reduced form model.
% endo_steady_state [vector] steady state value for endogenous variables
% exo_steady_state [vector] steady state value for exogenous variables
@ -30,7 +30,7 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% * get_perturbation_params_derivs.m (previously getH.m)
% -------------------------------------------------------------------------
% This function calls
% * [M.fname,'.dynamic']
% * [M_.fname,'.dynamic']
% * resol
% * dyn_vech
% =========================================================================
@ -53,9 +53,9 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% =========================================================================
%% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters
M = set_all_parameters(params,estim_params,M);
[dr,info,M.params] = compute_decision_rules(M,options,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
Sigma_e = M.Sigma_e;
M_ = set_all_parameters(params,estim_params,M_);
[dr,info,M_.params] = compute_decision_rules(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
Sigma_e = M_.Sigma_e;
if info(1) > 0
% there are errors in the solution algorithm
@ -64,10 +64,10 @@ if info(1) > 0
else
ys = dr.ys; %steady state of model variables in declaration order
ghx = dr.ghx; ghu = dr.ghu;
if options.order > 1
if options_.order > 1
ghxx = dr.ghxx; ghxu = dr.ghxu; ghuu = dr.ghuu; ghs2 = dr.ghs2;
end
if options.order > 2
if options_.order > 2
ghxxx = dr.ghxxx; ghxxu = dr.ghxxu; ghxuu = dr.ghxuu; ghxss = dr.ghxss; ghuuu = dr.ghuuu; ghuss = dr.ghuss;
end
end
@ -76,35 +76,35 @@ Yss = ys(dr.order_var); %steady state of model variables in DR order
%% out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]
if strcmp(outputflag,'perturbation_solution')
out = [Sigma_e(:); ghx(:); ghu(:)];
if options.order > 1
if options_.order > 1
out = [out; ghxx(:); ghxu(:); ghuu(:); ghs2(:);];
end
if options.order > 2
if options_.order > 2
out = [out; ghxxx(:); ghxxu(:); ghxuu(:); ghuuu(:); ghxss(:); ghuss(:)];
end
end
%% out = [Yss; vec(g1); vec(g2); vec(g3)]; of all endogenous variables, in DR order
if strcmp(outputflag,'dynamic_model')
[I,~] = find(M.lead_lag_incidence'); %I is used to evaluate dynamic model files
if options.order == 1
[~, g1] = feval([M.fname,'.dynamic'], ys(I), exo_steady_state', M.params, ys, 1);
[I,~] = find(M_.lead_lag_incidence'); %I is used to evaluate dynamic model files
if options_.order == 1
[~, g1] = feval([M_.fname,'.dynamic'], ys(I), exo_steady_state', M_.params, ys, 1);
out = [Yss; g1(:)];
elseif options.order == 2
[~, g1, g2] = feval([M.fname,'.dynamic'], ys(I), exo_steady_state', M.params, ys, 1);
elseif options_.order == 2
[~, g1, g2] = feval([M_.fname,'.dynamic'], ys(I), exo_steady_state', M_.params, ys, 1);
out = [Yss; g1(:); g2(:)];
elseif options.order == 3
[~, g1, g2, g3] = feval([M.fname,'.dynamic'], ys(I), exo_steady_state', M.params, ys, 1);
g3 = unfold_g3(g3, length(ys(I))+M.exo_nbr);
elseif options_.order == 3
[~, g1, g2, g3] = feval([M_.fname,'.dynamic'], ys(I), exo_steady_state', M_.params, ys, 1);
g3 = unfold_g3(g3, length(ys(I))+M_.exo_nbr);
out = [Yss; g1(:); g2(:); g3(:)];
end
end
%% out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')]; in DR order, where A and B are Kalman transition matrices
if strcmp(outputflag,'Kalman_Transition')
if options.order == 1
KalmanA = zeros(M.endo_nbr,M.endo_nbr);
KalmanA(:,M.nstatic+(1:M.nspred)) = ghx;
if options_.order == 1
KalmanA = zeros(M_.endo_nbr,M_.endo_nbr);
KalmanA(:,M_.nstatic+(1:M_.nspred)) = ghx;
Om = ghu*Sigma_e*transpose(ghu);
out = [Yss; KalmanA(:); dyn_vech(Om)];
else

28
matlab/histvalf.m
View File

@ -1,5 +1,5 @@
function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
%function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
function [endo_histval, exo_histval, exo_det_histval] = histvalf(M_, options_)
%function [endo_histval, exo_histval, exo_det_histval] = histvalf(M_, options_)
% Sets initial values for simulation using values contained in `fname`, a
% file possibly created by a call to `smoother2histval`
%
@ -30,33 +30,33 @@ function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
series = histvalf_initvalf('HISTVAL', M, options);
k = M.orig_maximum_lag - M.maximum_lag + 1;
series = histvalf_initvalf('HISTVAL', M_, options_);
k = M_.orig_maximum_lag - M_.maximum_lag + 1;
if ~isoctave && matlab_ver_less_than('9.7')
% Workaround for MATLAB bug described in dseries#45
% The solution is to avoid using the "end" keyword
myend = nobs(series);
endo_histval = series{M.endo_names{:}}.data(k:myend, :)';
endo_histval = series{M_.endo_names{:}}.data(k:myend, :)';
exo_histval = [];
if M.exo_nbr
exo_histval = series{M.exo_names{:}}.data(k:myend, :)';
if M_.exo_nbr
exo_histval = series{M_.exo_names{:}}.data(k:myend, :)';
end
exo_det_histval = [];
if M.exo_det_nbr
exo_det_histval = series{M.exo_names{:}}.data(k:myend, :)';
if M_.exo_det_nbr
exo_det_histval = series{M_.exo_names{:}}.data(k:myend, :)';
end
else
endo_histval = series{M.endo_names{:}}.data(k:end, :)';
endo_histval = series{M_.endo_names{:}}.data(k:end, :)';
exo_histval = [];
if M.exo_nbr
exo_histval = series{M.exo_names{:}}.data(k:end, :)';
if M_.exo_nbr
exo_histval = series{M_.exo_names{:}}.data(k:end, :)';
end
exo_det_histval = [];
if M.exo_det_nbr
exo_det_histval = series{M.exo_names{:}}.data(k:end, :)';
if M_.exo_det_nbr
exo_det_histval = series{M_.exo_names{:}}.data(k:end, :)';
end
end

34
matlab/identification_numerical_objective.m
View File

@ -1,5 +1,5 @@
function out = identification_numerical_objective(params, outputflag, estim_params, M, options, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state)
%function out = identification_numerical_objective(params, outputflag, %estim_params, M, options, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state)
function out = identification_numerical_objective(params, outputflag, estim_params_, M_, options_, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state)
% out = identification_numerical_objective(params, outputflag, estim_params_, M_, options_, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Objective function to compute numerically the Jacobians used for identification analysis
% Previously this function was called thet2tau.m
@ -8,9 +8,9 @@ function out = identification_numerical_objective(params, outputflag, estim_para
% params: [vector] parameter values at which to evaluate objective function
% stderr parameters come first, corr parameters second, model parameters third
% outputflag: [integer] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% options: [structure] storing the options
% estim_params_: [structure] storing the estimation information
% M_: [structure] storing the model information
% options_: [structure] storing the options
% indpmodel: [vector] Index of model parameters
% indpstderr: [vector] Index of stderr parameters
% indvar: [vector] Index of selected or observed variables
@ -35,7 +35,7 @@ function out = identification_numerical_objective(params, outputflag, estim_para
% * get_identification_jacobians.m (previously getJJ.m)
% -------------------------------------------------------------------------
% This function calls
% * [M.fname,'.dynamic']
% * [M_.fname,'.dynamic']
% * dyn_vech
% * resol
% * vec
@ -61,26 +61,26 @@ function out = identification_numerical_objective(params, outputflag, estim_para
%% Update stderr, corr and model parameters
%note that if no estimated_params_block is given, then all stderr and model parameters are selected but no corr parameters
if length(params) > length(indpmodel)
if isempty(indpstderr)==0 && isempty(estim_params.var_exo) %if there are stderr parameters but no estimated_params_block
if isempty(indpstderr)==0 && isempty(estim_params_.var_exo) %if there are stderr parameters but no estimated_params_block
%provide temporary necessary information for stderr parameters
estim_params.nvx = length(indpstderr);
estim_params.var_exo = indpstderr';
estim_params_.nvx = length(indpstderr);
estim_params_.var_exo = indpstderr';
end
if isempty(indpmodel)==0 && isempty(estim_params.param_vals) %if there are model parameters but no estimated_params_block
if isempty(indpmodel)==0 && isempty(estim_params_.param_vals) %if there are model parameters but no estimated_params_block
%provide temporary necessary information for model parameters
estim_params.np = length(indpmodel);
estim_params.param_vals = indpmodel';
estim_params_.np = length(indpmodel);
estim_params_.param_vals = indpmodel';
end
M = set_all_parameters(params,estim_params,M); %this function can only be used if there is some information in estim_params
M_ = set_all_parameters(params,estim_params_,M_); %this function can only be used if there is some information in estim_params_
else
%if there are only model parameters, we don't need to use set_all_parameters
M.params(indpmodel) = params;
M_.params(indpmodel) = params;
end
%% compute Kalman transition matrices and steady state with updated parameters
[dr,info,M.params] = compute_decision_rules(M,options,dr, steady_state, exo_steady_state, exo_det_steady_state);
options = rmfield(options,'options_ident');
pruned = pruned_state_space_system(M, options, dr, indvar, nlags, useautocorr, 0);
[dr,info,M_.params] = compute_decision_rules(M_,options_,dr, steady_state, exo_steady_state, exo_det_steady_state);
options_ = rmfield(options_,'options_ident');
pruned = pruned_state_space_system(M_, options_, dr, indvar, nlags, useautocorr, 0);
%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
if outputflag == 1

16
matlab/k_order_pert.m
View File

@ -1,4 +1,4 @@
function [dr,info] = k_order_pert(dr,M,options)
function [dr,info] = k_order_pert(dr,M_,options_)
% Compute decision rules using the k-order DLL from Dynare++
% Copyright © 2009-2023 Dynare Team
@ -20,18 +20,18 @@ function [dr,info] = k_order_pert(dr,M,options)
info = 0;
order = options.order;
order = options_.order;
if order>1 && options.loglinear
if order>1 && options_.loglinear
error('The loglinear-option currently only works at order 1')
end
if M.maximum_endo_lead == 0 && order>1
if M_.maximum_endo_lead == 0 && order>1
error(['2nd and 3rd order approximation not implemented for purely ' ...
'backward models'])
end
try
[dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,M,options);
[dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,M_,options_);
catch ME
disp(ME.message)
info(1)=9;
@ -44,7 +44,7 @@ for i = 0:order
gname = [ 'g_' num2str(i) ];
dr.(gname) = dynpp_derivs.(gname);
end
if options.pruning
if options_.pruning
dr.pruning = dynpp_derivs.pruning;
end
@ -54,8 +54,8 @@ end
dr.ghx = dyn_derivs.gy;
dr.ghu = dyn_derivs.gu;
if options.loglinear
k = find(dr.kstate(:,2) <= M.maximum_endo_lag+1);
if options_.loglinear
k = find(dr.kstate(:,2) <= M_.maximum_endo_lag+1);
klag = dr.kstate(k,[1 2]);
k1 = dr.order_var;
dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ...

10
matlab/metropolis_run_analysis.m
View File

@ -1,9 +1,9 @@
function metropolis_run_analysis(M,basetopt,j)
%function metropolis_run_analysis(M)
function metropolis_run_analysis(M_,basetopt,j)
%function metropolis_run_analysis(M_,basetopt,j
% analizes Metropolis runs
%
% INPUTS
% M: (struct) Model structure
% M_: (struct) Model structure
% basetopt: (struct) Estimated parameters structure
% j: (int) Index of estimated paramter
%
@ -30,7 +30,7 @@ function metropolis_run_analysis(M,basetopt,j)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
load([M.fname '/metropolis/' M.fname '_mh_history'])
load([M_.fname '/metropolis/' M_.fname '_mh_history'])
nblck = record.Nblck;
ndraws = sum(record.MhDraws(:,1));
@ -38,7 +38,7 @@ logPost = [];
params = [];
blck = 1;
for i=1:record.LastFileNumber
fname = [M.fname '/metropolis/' M.fname '_mh' int2str(i) '_blck' ...
fname = [M_.fname '/metropolis/' M_.fname '_mh' int2str(i) '_blck' ...
int2str(blck) '.mat'];
if exist(fname,'file')
o=load(fname);

8
matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
View File

@ -23,15 +23,15 @@ function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_functi
% maximum_lag [scalar] maximum lag present in the model
% T [scalar] number of simulation periods
% ny [scalar] number of endogenous variables
% i_cols [double] indices of variables appearing in M.lead_lag_incidence
% i_cols [double] indices of variables appearing in M_.lead_lag_incidence
% and that need to be passed to _dynamic-file
% i_cols_J1 [double] indices of contemporaneous and forward looking variables
% appearing in M.lead_lag_incidence
% appearing in M_.lead_lag_incidence
% i_cols_1 [double] indices of contemporaneous and forward looking variables in
% M.lead_lag_incidence in dynamic Jacobian (relevant in first period)
% M_.lead_lag_incidence in dynamic Jacobian (relevant in first period)
% i_cols_T [double] columns of dynamic Jacobian related to contemporaneous and backward-looking
% variables (relevant in last period)
% i_cols_j [double] indices of variables in M.lead_lag_incidence
% i_cols_j [double] indices of variables in M_.lead_lag_incidence
% in dynamic Jacobian (relevant in intermediate periods)
% eq_index [double] N*1 array, index vector describing residual mapping resulting
% from complementarity setup

80
matlab/perfect-foresight-models/private/initialize_stacked_problem.m
View File

@ -1,36 +1,36 @@
function [options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i_cols_0, i_cols_J0, dynamicmodel] = ...
initialize_stacked_problem(endogenousvariables, options, M, steadystate_y)
function [options_, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i_cols_0, i_cols_J0, dynamicmodel] = ...
initialize_stacked_problem(endogenousvariables, options_, M_, steadystate_y)
% Sets up the stacked perfect foresight problem for use with dynare_solve.m
%
% INPUTS
% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial guess).
% - options [struct] contains various options.
% - M [struct] contains a description of the model.
% - options_ [struct] contains various options.
% - M_ [struct] contains a description of the model.
% - steadystate_y [double] N*1 array, steady state for the endogenous variables.
%
% OUTPUTS
% - options [struct] contains various options.
% - options_ [struct] contains various options.
% - y0 [double] N*1 array, initial conditions for the endogenous variables
% - yT [double] N*1 array, terminal conditions for the endogenous variables
% - z [double] T*M array, paths for the exogenous variables.
% - i_cols [double] indices of variables appearing in M.lead_lag_incidence
% - i_cols [double] indices of variables appearing in M_.lead_lag_incidence
% and that need to be passed to _dynamic-file
% - i_cols_J1 [double] indices of contemporaneous and forward looking variables
% appearing in M.lead_lag_incidence
% appearing in M_.lead_lag_incidence
% - i_cols_T [double] columns of dynamic Jacobian related to
% contemporaneous and backward-looking
% variables (relevant in last period)
% - i_cols_j [double] indices of variables in M.lead_lag_incidence
% - i_cols_j [double] indices of variables in M_.lead_lag_incidence
% in dynamic Jacobian (relevant in intermediate periods)
% - i_cols_1 [double] indices of contemporaneous and forward looking variables in
% M.lead_lag_incidence in dynamic Jacobian (relevant in first period)
% - i_cols_0 [double] indices of contemporaneous variables in M.lead_lag_incidence in dynamic
% M_.lead_lag_incidence in dynamic Jacobian (relevant in first period)
% - i_cols_0 [double] indices of contemporaneous variables in M_.lead_lag_incidence in dynamic
% Jacobian (relevant in problems with periods=1)
% - i_cols_J0 [double] indices of contemporaneous variables appearing in M.lead_lag_incidence (relevant in problems with periods=1)
% - i_cols_J0 [double] indices of contemporaneous variables appearing in M_.lead_lag_incidence (relevant in problems with periods=1)
% - dynamicmodel [handle] function handle to _dynamic-file
% Copyright © 2015-2020 Dynare Team
% Copyright © 2015-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -47,53 +47,53 @@ function [options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
periods = options.periods;
if (options.solve_algo == 10)
if ~isfield(options.lmmcp,'lb')
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M,options.ramsey_policy);
if options.linear_approximation
periods = options_.periods;
if (options_.solve_algo == 10)
if ~isfield(options_.lmmcp,'lb')
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M_,options_.ramsey_policy);
if options_.linear_approximation
lb = lb - steadystate_y;
ub = ub - steadystate_y;
end
options.lmmcp.lb = repmat(lb,periods,1);
options.lmmcp.ub = repmat(ub,periods,1);
options_.lmmcp.lb = repmat(lb,periods,1);
options_.lmmcp.ub = repmat(ub,periods,1);
end
elseif (options.solve_algo == 11)
if ~isfield(options.mcppath,'lb')
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M,options.ramsey_policy);
if options.linear_approximation
elseif (options_.solve_algo == 11)
if ~isfield(options_.mcppath,'lb')
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M_,options_.ramsey_policy);
if options_.linear_approximation
lb = lb - steadystate_y;
ub = ub - steadystate_y;
end
options.mcppath.lb = repmat(lb,periods,1);
options.mcppath.ub = repmat(ub,periods,1);
options_.mcppath.lb = repmat(lb,periods,1);
options_.mcppath.ub = repmat(ub,periods,1);
end
end
if M.maximum_lag > 0
y0 = endogenousvariables(:, M.maximum_lag);
if M_.maximum_lag > 0
y0 = endogenousvariables(:, M_.maximum_lag);
else
y0 = NaN(M.endo_nbr, 1);
y0 = NaN(M_.endo_nbr, 1);
end
if M.maximum_lead > 0
yT = endogenousvariables(:, M.maximum_lag+periods+1);
if M_.maximum_lead > 0
yT = endogenousvariables(:, M_.maximum_lag+periods+1);
else
yT = NaN(M.endo_nbr, 1);
yT = NaN(M_.endo_nbr, 1);
end
z = endogenousvariables(:,M.maximum_lag+(1:periods));
illi = M.lead_lag_incidence';
z = endogenousvariables(:,M_.maximum_lag+(1:periods));
illi = M_.lead_lag_incidence';
[i_cols,~,i_cols_j] = find(illi(:));
if M.maximum_endo_lag == 0
i_cols = i_cols + M.endo_nbr;
if M_.maximum_endo_lag == 0
i_cols = i_cols + M_.endo_nbr;
end
illi = illi(:,(1+M.maximum_endo_lag):(1+M.maximum_endo_lag+M.maximum_endo_lead));
illi = illi(:,(1+M_.maximum_endo_lag):(1+M_.maximum_endo_lag+M_.maximum_endo_lead));
[i_cols_J1,~,i_cols_1] = find(illi(:));
i_cols_T = nonzeros(M.lead_lag_incidence(1:(1+M.maximum_endo_lag),:)');
i_cols_T = nonzeros(M_.lead_lag_incidence(1:(1+M_.maximum_endo_lag),:)');
if periods==1
i_cols_0 = nonzeros(M.lead_lag_incidence(1+M.maximum_endo_lag,:)');
i_cols_J0 = find(M.lead_lag_incidence(1+M.maximum_endo_lag,:)');
i_cols_0 = nonzeros(M_.lead_lag_incidence(1+M_.maximum_endo_lag,:)');
i_cols_J0 = find(M_.lead_lag_incidence(1+M_.maximum_endo_lag,:)');
else
i_cols_0 = [];
i_cols_J0 = [];
end
dynamicmodel = str2func([M.fname,'.dynamic']);
dynamicmodel = str2func([M_.fname,'.dynamic']);

70
matlab/perfect-foresight-models/sim1.m
View File

@ -1,16 +1,16 @@
function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M_, options_)
% Performs deterministic simulations with lead or lag of one period, using
% a basic Newton solver on sparse matrices.
% Uses perfect_foresight_problem DLL to construct the stacked problem.
%
% INPUTS
% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (initial condition + initial guess + terminal condition).
% - endogenousvariables [double] N*(T+M_.maximum_lag+M_.maximum_lead) array, paths for the endogenous variables (initial condition + initial guess + terminal condition).
% - exogenousvariables [double] T*M array, paths for the exogenous variables.
% - steadystate [double] N*1 array, steady state for the endogenous variables.
% - M [struct] contains a description of the model.
% - options [struct] contains various options.
% - M_ [struct] contains a description of the model.
% - options_ [struct] contains various options.
% OUTPUTS
% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model).
% - endogenousvariables [double] N*(T+M_.maximum_lag+M_.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model).
% - success [logical] Whether a solution was found
% - err [double] ∞-norm of the residual
% - iter [integer] Number of iterations
@ -32,25 +32,25 @@ function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, e
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
verbose = options.verbosity && ~options.noprint;
verbose = options_.verbosity && ~options_.noprint;
ny = M.endo_nbr;
periods = options.periods;
vperiods = periods*ones(1,options.simul.maxit);
ny = M_.endo_nbr;
periods = options_.periods;
vperiods = periods*ones(1,options_.simul.maxit);
if M.maximum_lag > 0
y0 = endogenousvariables(:, M.maximum_lag);
if M_.maximum_lag > 0
y0 = endogenousvariables(:, M_.maximum_lag);
else
y0 = NaN(ny, 1);
end
if M.maximum_lead > 0
yT = endogenousvariables(:, M.maximum_lag+periods+1);
if M_.maximum_lead > 0
yT = endogenousvariables(:, M_.maximum_lag+periods+1);
else
yT = NaN(ny, 1);
end
y = reshape(endogenousvariables(:, M.maximum_lag+(1:periods)), ny*periods, 1);
y = reshape(endogenousvariables(:, M_.maximum_lag+(1:periods)), ny*periods, 1);
stop = false;
@ -63,13 +63,13 @@ end
h1 = clock;
for iter = 1:options.simul.maxit
for iter = 1:options_.simul.maxit
h2 = clock;
[res, A] = perfect_foresight_problem(y, y0, yT, exogenousvariables, M.params, steadystate, periods, M, options);
[res, A] = perfect_foresight_problem(y, y0, yT, exogenousvariables, M_.params, steadystate, periods, M_, options_);
% A is the stacked Jacobian with period x equations alongs the rows and
% periods times variables (in declaration order) along the columns
if options.debug && iter==1
if options_.debug && iter==1
[row,col]=find(A);
row=setdiff(1:periods*ny,row);
column=setdiff(1:periods*ny,col);
@ -91,14 +91,14 @@ for iter = 1:options.simul.maxit
end
end
end
if options.endogenous_terminal_period && iter > 1
if options_.endogenous_terminal_period && iter > 1
for it = 1:periods
if max(abs(res((it-1)*ny+(1:ny)))) < options.dynatol.f/1e7
if max(abs(res((it-1)*ny+(1:ny)))) < options_.dynatol.f/1e7
if it < periods
res = res(1:(it*ny));
A = A(1:(it*ny), 1:(it*ny));
yT = y(it*ny+(1:ny));
endogenousvariables(:, M.maximum_lag+((it+1):periods)) = reshape(y(it*ny+(1:(ny*(periods-it)))), ny, periods-it);
endogenousvariables(:, M_.maximum_lag+((it+1):periods)) = reshape(y(it*ny+(1:(ny*(periods-it)))), ny, periods-it);
y = y(1:(it*ny));
periods = it;
end
@ -109,7 +109,7 @@ for iter = 1:options.simul.maxit
end
err = max(abs(res));
if options.debug
if options_.debug
fprintf('\nLargest absolute residual at iteration %d: %10.3f\n',iter,err);
if any(isnan(res)) || any(isinf(res)) || any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables)))
fprintf('\nWARNING: NaN or Inf detected in the residuals or endogenous variables.\n');
@ -119,28 +119,28 @@ for iter = 1:options.simul.maxit
if verbose
fprintf('Iter: %d,\t err. = %g,\t time = %g\n', iter, err, etime(clock,h2));
end
if err < options.dynatol.f
if err < options_.dynatol.f
stop = true;
break
end
if options.simul.robust_lin_solve
dy = -lin_solve_robust(A, res, verbose, options);
if options_.simul.robust_lin_solve
dy = -lin_solve_robust(A, res, verbose, options_);
else
dy = -lin_solve(A, res, verbose);
end
if any(isnan(dy)) || any(isinf(dy))
if verbose
display_critical_variables(reshape(dy,[ny periods])', M, options.noprint);
display_critical_variables(reshape(dy,[ny periods])', M_, options_.noprint);
end
end
y = y + dy;
end
endogenousvariables(:, M.maximum_lag+(1:periods)) = reshape(y, ny, periods);
endogenousvariables(:, M_.maximum_lag+(1:periods)) = reshape(y, ny, periods);
if options.endogenous_terminal_period
periods = options.periods;
err = evaluate_max_dynamic_residual(str2func([M.fname,'.dynamic']), endogenousvariables, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_endo_lag, M.lead_lag_incidence);
if options_.endogenous_terminal_period
periods = options_.periods;
err = evaluate_max_dynamic_residual(str2func([M_.fname,'.dynamic']), endogenousvariables, exogenousvariables, M_.params, steadystate, periods, ny, M_.maximum_endo_lag, M_.lead_lag_incidence);
end
if stop
@ -152,7 +152,7 @@ if stop
skipline()
fprintf('Total time of simulation: %g.\n', etime(clock,h1))
disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.')
display_critical_variables(reshape(dy,[ny periods])', M, options.noprint);
display_critical_variables(reshape(dy,[ny periods])', M_, options_.noprint);
disp('There is most likely something wrong with your model. Try model_diagnostics or another simulation method.')
printline(105)
end
@ -191,7 +191,7 @@ if relres > 1e-6 && verbose
fprintf('WARNING : Failed to find a solution to the linear system.\n');
end
function [ x, flag, relres ] = lin_solve_robust(A, b ,verbose, options)
function [ x, flag, relres ] = lin_solve_robust(A, b ,verbose, options_)
if norm(b) < sqrt(eps) % then x = 0 is a solution
x = 0;
flag = 0;
@ -206,7 +206,7 @@ if flag == 0
return
end
if ~options.noprint
if ~options_.noprint
disp( relres );
end
@ -261,7 +261,7 @@ if flag ~= 0 && verbose
fprintf('WARNING : Failed to find a solution to the linear system\n');
end
function display_critical_variables(dyy, M, noprint)
function display_critical_variables(dyy, M_, noprint)
if noprint
return
@ -269,14 +269,14 @@ end
if any(isnan(dyy))
indx = find(any(isnan(dyy)));
endo_names= M.endo_names(indx);
endo_names= M_.endo_names(indx);
disp('Last iteration provided NaN for the following variables:')
fprintf('%s, ', endo_names{:}),
fprintf('\n'),
end
if any(isinf(dyy))
indx = find(any(isinf(dyy)));
endo_names = M.endo_names(indx);
endo_names = M_.endo_names(indx);
disp('Last iteration diverged (Inf) for the following variables:')
fprintf('%s, ', endo_names{:}),
fprintf('\n'),

32
matlab/perfect-foresight-models/sim1_lbj.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, success, err, iter] = sim1_lbj(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success, err, iter] = sim1_lbj(endogenousvariables, exogenousvariables, steadystate, M_, options_)
% Performs deterministic simulations with lead or lag on one period using the historical LBJ algorithm
%
@ -32,7 +32,7 @@ function [endogenousvariables, success, err, iter] = sim1_lbj(endogenousvariable
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
lead_lag_incidence = M.lead_lag_incidence;
lead_lag_incidence = M_.lead_lag_incidence;
ny = size(endogenousvariables,1);
nyp = nnz(lead_lag_incidence(1,:));
@ -51,17 +51,17 @@ iz = [1:ny+nyp+nyf];
function [r, g1] = bytecode_wrapper(y, xpath, params, ys, it_)
ypath = NaN(ny, 3);
ypath(find(lead_lag_incidence')) = y;
[r, s] = bytecode('dynamic', 'evaluate', M, options, ypath, xpath(it_+(-1:1), :), params, ys, 1);
[r, s] = bytecode('dynamic', 'evaluate', M_, options_, ypath, xpath(it_+(-1:1), :), params, ys, 1);
g1 = s.g1;
end
if options.bytecode
if options_.bytecode
dynamicmodel = @bytecode_wrapper;
else
dynamicmodel = str2func(sprintf('%s.dynamic', M.fname));
dynamicmodel = str2func(sprintf('%s.dynamic', M_.fname));
end
verbose = options.verbosity;
verbose = options_.verbosity;
if verbose
printline(56)
@ -69,37 +69,37 @@ if verbose
skipline()
end
it_init = M.maximum_lag+1;
it_init = M_.maximum_lag+1;
h1 = clock;
for iter = 1:options.simul.maxit
for iter = 1:options_.simul.maxit
h2 = clock;
c = zeros(ny*options.periods, nrc);
c = zeros(ny*options_.periods, nrc);
it_ = it_init;
z = [endogenousvariables(iyp,it_-1) ; endogenousvariables(:,it_) ; endogenousvariables(iyf,it_+1)];
[d1, jacobian] = dynamicmodel(z, exogenousvariables, M.params, steadystate, it_);
[d1, jacobian] = dynamicmodel(z, exogenousvariables, M_.params, steadystate, it_);
jacobian = [jacobian(:,iz), -d1];
ic = [1:ny];
icp = iyp;
c (ic,:) = jacobian(:,is)\jacobian(:,isf1);
for it_ = it_init+(1:options.periods-1)
for it_ = it_init+(1:options_.periods-1)
z = [endogenousvariables(iyp,it_-1) ; endogenousvariables(:,it_) ; endogenousvariables(iyf,it_+1)];
[d1, jacobian] = dynamicmodel(z, exogenousvariables, M.params, steadystate, it_);
[d1, jacobian] = dynamicmodel(z, exogenousvariables, M_.params, steadystate, it_);
jacobian = [jacobian(:,iz), -d1];
jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:);
ic = ic + ny;
icp = icp + ny;
c (ic,:) = jacobian(:,is)\jacobian(:,isf1);
end
c = bksup1(c, ny, nrc, iyf, options.periods);
c = reshape(c, ny, options.periods);
endogenousvariables(:,it_init+(0:options.periods-1)) = endogenousvariables(:,it_init+(0:options.periods-1))+c;
c = bksup1(c, ny, nrc, iyf, options_.periods);
c = reshape(c, ny, options_.periods);
endogenousvariables(:,it_init+(0:options_.periods-1)) = endogenousvariables(:,it_init+(0:options_.periods-1))+c;
err = max(max(abs(c)));
if verbose
str = sprintf('Iter: %s,\t err. = %s, \t time = %s', num2str(iter), num2str(err), num2str(etime(clock, h2)));
disp(str);
end
if err < options.dynatol.f
if err < options_.dynatol.f
stop = true;
if verbose
skipline()

54
matlab/perfect-foresight-models/sim1_linear.m
View File

@ -1,5 +1,5 @@
function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M_, options_)
% [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M_, options_)
% Solves a linear approximation of a perfect foresight model using sparse matrix.
%
% INPUTS
@ -7,8 +7,8 @@ function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables,
% - exogenousvariables [double] T*M array, paths for the exogenous variables.
% - steadystate_y [double] N*1 array, steady state for the endogenous variables.
% - steadystate_x [double] M*1 array, steady state for the exogenous variables.
% - M [struct] contains a description of the model.
% - options [struct] contains various options.
% - M_ [struct] contains a description of the model.
% - options_ [struct] contains various options.
%
% OUTPUTS
% - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
@ -21,7 +21,7 @@ function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables,
% - T is the number of periods (including initial and/or terminal conditions).
%
% REMARKS
% - The structure `M` describing the structure of the model, must contain the
% - The structure `M_` describing the structure of the model, must contain the
% following informations:
% + lead_lag_incidence, incidence matrix (given by the preprocessor).
% + endo_nbr, number of endogenous variables (including aux. variables).
@ -31,7 +31,7 @@ function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables,
% + params, values of model's parameters.
% + fname, name of the model.
% + NNZDerivatives, number of non zero elements in the jacobian of the dynamic model.
% - The structure `options`, must contain the following options:
% - The structure `options_`, must contain the following options:
% + verbosity, controls the quantity of information displayed.
% + periods, the number of periods in the perfect foresight model.
% + debug.
@ -56,15 +56,15 @@ function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
verbose = options.verbosity;
verbose = options_.verbosity;
lead_lag_incidence = M.lead_lag_incidence;
lead_lag_incidence = M_.lead_lag_incidence;
ny = M.endo_nbr;
nx = M.exo_nbr;
ny = M_.endo_nbr;
nx = M_.exo_nbr;
maximum_lag = M.maximum_lag;
max_lag = M.maximum_endo_lag;
maximum_lag = M_.maximum_lag;
max_lag = M_.maximum_endo_lag;
nyp = nnz(lead_lag_incidence(1,:));
ny0 = nnz(lead_lag_incidence(2,:));
@ -72,9 +72,9 @@ nyf = nnz(lead_lag_incidence(3,:));
nd = nyp+ny0+nyf; % size of y (first argument passed to the dynamic file).
periods = options.periods;
periods = options_.periods;
params = M.params;
params = M_.params;
% Indices in A.
ip = find(lead_lag_incidence(1,:)');
@ -108,32 +108,32 @@ if verbose
skipline()
end
dynamicmodel = str2func([M.fname,'.dynamic']);
dynamicmodel = str2func([M_.fname,'.dynamic']);
z = steadystate_y([ip; ic; in]);
x = repmat(transpose(steadystate_x), 1+M.maximum_exo_lag+M.maximum_exo_lead, 1);
x = repmat(transpose(steadystate_x), 1+M_.maximum_exo_lag+M_.maximum_exo_lead, 1);
% Evaluate the Jacobian of the dynamic model at the deterministic steady state.
[d1, jacobian] = dynamicmodel(z, x, params, steadystate_y, M.maximum_exo_lag+1);
if options.debug
[d1, jacobian] = dynamicmodel(z, x, params, steadystate_y, M_.maximum_exo_lag+1);
if options_.debug
column=find(all(jacobian==0,1));
if ~isempty(column)
fprintf('The dynamic Jacobian is singular. The problem derives from:\n')
for iter=1:length(column)
[var_row,var_index]=find(M.lead_lag_incidence==column(iter));
[var_row,var_index]=find(M_.lead_lag_incidence==column(iter));
if var_row==2
fprintf('The derivative with respect to %s being 0 for all equations.\n',M.endo_names{var_index})
fprintf('The derivative with respect to %s being 0 for all equations.\n',M_.endo_names{var_index})
elseif var_row==1
fprintf('The derivative with respect to the lag of %s being 0 for all equations.\n',M.endo_names{var_index})
fprintf('The derivative with respect to the lag of %s being 0 for all equations.\n',M_.endo_names{var_index})
elseif var_row==3
fprintf('The derivative with respect to the lead of %s being 0 for all equations.\n',M.endo_names{var_index})
fprintf('The derivative with respect to the lead of %s being 0 for all equations.\n',M_.endo_names{var_index})
end
end
end
end
% Check that the dynamic model was evaluated at the steady state.
if ~options.steadystate.nocheck && max(abs(d1))>options.solve_tolf
if ~options_.steadystate.nocheck && max(abs(d1))>options_.solve_tolf
error('Jacobian is not evaluated at the steady state!')
end
@ -155,7 +155,7 @@ iv = 1:length(vv);
res = zeros(periods*ny, 1);
% Initialize the sparse Jacobian.
iA = zeros(periods*M.NNZDerivatives(1), 3);
iA = zeros(periods*M_.NNZDerivatives(1), 3);
h2 = clock;
i_rows = (1:ny)';
@ -190,7 +190,7 @@ end
% Evaluation of the maximum residual at the initial guess (steady state for the endogenous variables).
err = max(abs(res));
if options.debug
if options_.debug
fprintf('\nLargest absolute residual at iteration %d: %10.3f\n', 1, err);
if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y))
fprintf('\nWARNING: NaN or Inf detected in the residuals or endogenous variables.\n');
@ -239,7 +239,7 @@ end
if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) || ~isreal(res) || ~isreal(Y)
success = false; % NaN or Inf occurred
endogenousvariables = reshape(Y, ny, periods+maximum_lag+M.maximum_lead);
endogenousvariables = reshape(Y, ny, periods+maximum_lag+M_.maximum_lead);
if verbose
skipline()
if ~isreal(res) || ~isreal(Y)
@ -251,7 +251,7 @@ if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) || ~isre
end
else
success = true; % Convergency obtained.
endogenousvariables = bsxfun(@plus, reshape(Y, ny, periods+maximum_lag+M.maximum_lead), steadystate_y);
endogenousvariables = bsxfun(@plus, reshape(Y, ny, periods+maximum_lag+M_.maximum_lead), steadystate_y);
end
if verbose

44
matlab/perfect-foresight-models/sim1_purely_backward.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, success] = sim1_purely_backward(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success] = sim1_purely_backward(endogenousvariables, exogenousvariables, steadystate, M_, options_)
% Performs deterministic simulation of a purely backward model
@ -19,41 +19,41 @@ function [endogenousvariables, success] = sim1_purely_backward(endogenousvariabl
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if ismember(options.solve_algo, [12,14])
[funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
if ismember(options_.solve_algo, [12,14])
[funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_);
else
dynamic_resid = str2func([M.fname '.sparse.dynamic_resid']);
dynamic_g1 = str2func([M.fname '.sparse.dynamic_g1']);
dynamic_resid = str2func([M_.fname '.sparse.dynamic_resid']);
dynamic_g1 = str2func([M_.fname '.sparse.dynamic_g1']);
end
function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
% NB: do as few computations as possible inside this function, since it is
% called a very large number of times
y_dynamic(feedback_vars_idx) = z;
[~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
[~, ~, r, J] = feval(func, y_dynamic, x, M_.params, steadystate, ...
sparse_rowval, sparse_colval, sparse_colptr, T);
end
success = true;
for it = M.maximum_lag + (1:options.periods)
for it = M_.maximum_lag + (1:options_.periods)
y = endogenousvariables(:,it-1); % Values at previous period, also used as guess value for current period
x = exogenousvariables(it,:);
if ismember(options.solve_algo, [12,14])
T = NaN(M.block_structure.dyn_tmp_nbr);
y_dynamic = [y; y; NaN(M.endo_nbr, 1)];
for blk = 1:length(M.block_structure.block)
sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
if M.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
if ismember(options_.solve_algo, [12,14])
T = NaN(M_.block_structure.dyn_tmp_nbr);
y_dynamic = [y; y; NaN(M_.endo_nbr, 1)];
for blk = 1:length(M_.block_structure.block)
sparse_rowval = M_.block_structure.block(blk).g1_sparse_rowval;
sparse_colval = M_.block_structure.block(blk).g1_sparse_colval;
sparse_colptr = M_.block_structure.block(blk).g1_sparse_colptr;
if M_.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
[z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
options.simul.maxit, options.dynatol.f, ...
options.dynatol.x, options, ...
options_.simul.maxit, options_.dynatol.f, ...
options_.dynatol.x, options_, ...
feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
if check
success = false;
if options.debug
if options_.debug
dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
end
end
@ -61,15 +61,15 @@ for it = M.maximum_lag + (1:options.periods)
end
%% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
%% Also update the temporary terms vector.
[y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
[y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M_.params, ...
steadystate, sparse_rowval, sparse_colval, ...
sparse_colptr, T);
end
endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
endogenousvariables(:,it) = y_dynamic(M_.endo_nbr+(1:M_.endo_nbr));
else
[tmp, check, ~, ~, errorcode] = dynare_solve(@dynamic_backward_model_for_simulation, y, ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, dynamic_resid, dynamic_g1, y, x, M.params, steadystate, M.dynamic_g1_sparse_rowval, M.dynamic_g1_sparse_colval, M.dynamic_g1_sparse_colptr);
options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
options_, dynamic_resid, dynamic_g1, y, x, M_.params, steadystate, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, M_.dynamic_g1_sparse_colptr);
if check
success = false;
dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in period %i', errorcode, it)

46
matlab/perfect-foresight-models/sim1_purely_forward.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, success] = sim1_purely_forward(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success] = sim1_purely_forward(endogenousvariables, exogenousvariables, steadystate, M_, options_)
% Performs deterministic simulation of a purely forward model
% Copyright © 2012-2023 Dynare Team
@ -18,41 +18,41 @@ function [endogenousvariables, success] = sim1_purely_forward(endogenousvariable
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if ismember(options.solve_algo, [12,14])
[funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
if ismember(options_.solve_algo, [12,14])
[funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_);
else
dynamic_resid = str2func([M.fname '.sparse.dynamic_resid']);
dynamic_g1 = str2func([M.fname '.sparse.dynamic_g1']);
dynamic_resid = str2func([M_.fname '.sparse.dynamic_resid']);
dynamic_g1 = str2func([M_.fname '.sparse.dynamic_g1']);
end
function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
% NB: do as few computations as possible inside this function, since it is
% called a very large number of times
y_dynamic(feedback_vars_idx) = z;
[~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
[~, ~, r, J] = feval(func, y_dynamic, x, M_.params, steadystate, ...
sparse_rowval, sparse_colval, sparse_colptr, T);
end
success = true;
for it = options.periods:-1:1
for it = options_.periods:-1:1
yf = endogenousvariables(:,it+1); % Values at next period, also used as guess value for current period
x = exogenousvariables(it,:);
if ismember(options.solve_algo, [12,14])
T = NaN(M.block_structure.dyn_tmp_nbr);
y_dynamic = [NaN(M.endo_nbr, 1); yf; yf];
for blk = 1:length(M.block_structure.block)
sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
if M.block_structure.block(blk).Simulation_Type ~= 2 % Not an evaluate backward block
if ismember(options_.solve_algo, [12,14])
T = NaN(M_.block_structure.dyn_tmp_nbr);
y_dynamic = [NaN(M_.endo_nbr, 1); yf; yf];
for blk = 1:length(M_.block_structure.block)
sparse_rowval = M_.block_structure.block(blk).g1_sparse_rowval;
sparse_colval = M_.block_structure.block(blk).g1_sparse_colval;
sparse_colptr = M_.block_structure.block(blk).g1_sparse_colptr;
if M_.block_structure.block(blk).Simulation_Type ~= 2 % Not an evaluate backward block
[z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
options.simul.maxit, options.dynatol.f, ...
options.dynatol.x, options, ...
options_.simul.maxit, options_.dynatol.f, ...
options_.dynatol.x, options_, ...
feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
if check
success = false;
if options.debug
if options_.debug
dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
end
end
@ -60,21 +60,21 @@ for it = options.periods:-1:1
end
%% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
%% Also update the temporary terms vector.
[y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
[y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M_.params, ...
steadystate, sparse_rowval, sparse_colval, ...
sparse_colptr, T);
end
endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
endogenousvariables(:,it) = y_dynamic(M_.endo_nbr+(1:M_.endo_nbr));
else
[tmp, check, ~, ~, errorcode] = dynare_solve(@dynamic_forward_model_for_simulation, yf, ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, dynamic_resid, dynamic_g1, yf, x, M.params, steadystate, M.dynamic_g1_sparse_rowval, M.dynamic_g1_sparse_colval, M.dynamic_g1_sparse_colptr);
options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
options_, dynamic_resid, dynamic_g1, yf, x, M_.params, steadystate, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, M_.dynamic_g1_sparse_colptr);
if check
success = false;
dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
break
end
endogenousvariables(:,it) = tmp(1:M.endo_nbr);
endogenousvariables(:,it) = tmp(1:M_.endo_nbr);
end
end

46
matlab/perfect-foresight-models/sim1_purely_static.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, success] = sim1_purely_static(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success] = sim1_purely_static(endogenousvariables, exogenousvariables, steadystate, M_, options_)
% Performs deterministic simulation of a purely static model
@ -19,18 +19,18 @@ function [endogenousvariables, success] = sim1_purely_static(endogenousvariables
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if ismember(options.solve_algo, [12,14])
[funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
if ismember(options_.solve_algo, [12,14])
[funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_);
else
dynamic_resid = str2func([M.fname '.sparse.dynamic_resid']);
dynamic_g1 = str2func([M.fname '.sparse.dynamic_g1']);
dynamic_resid = str2func([M_.fname '.sparse.dynamic_resid']);
dynamic_g1 = str2func([M_.fname '.sparse.dynamic_g1']);
end
function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
% NB: do as few computations as possible inside this function, since it is
% called a very large number of times
y_dynamic(feedback_vars_idx) = z;
[~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
[~, ~, r, J] = feval(func, y_dynamic, x, M_.params, steadystate, ...
sparse_rowval, sparse_colval, sparse_colptr, T);
end
@ -38,23 +38,23 @@ success = true;
y = endogenousvariables(:,1);
for it = 1:options.periods
for it = 1:options_.periods
x = exogenousvariables(it,:);
if ismember(options.solve_algo, [12,14])
T = NaN(M.block_structure.dyn_tmp_nbr);
y_dynamic = [NaN(M.endo_nbr, 1); y; NaN(M.endo_nbr, 1)];
for blk = 1:length(M.block_structure.block)
sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
if M.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
if ismember(options_.solve_algo, [12,14])
T = NaN(M_.block_structure.dyn_tmp_nbr);
y_dynamic = [NaN(M_.endo_nbr, 1); y; NaN(M_.endo_nbr, 1)];
for blk = 1:length(M_.block_structure.block)
sparse_rowval = M_.block_structure.block(blk).g1_sparse_rowval;
sparse_colval = M_.block_structure.block(blk).g1_sparse_colval;
sparse_colptr = M_.block_structure.block(blk).g1_sparse_colptr;
if M_.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
[z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
options.simul.maxit, options.dynatol.f, ...
options.dynatol.x, options, ...
options_.simul.maxit, options_.dynatol.f, ...
options_.dynatol.x, options_, ...
feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
if check
success = false;
if options.debug
if options_.debug
dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
end
end
@ -62,18 +62,18 @@ for it = 1:options.periods
end
%% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
%% Also update the temporary terms vector.
[y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
[y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M_.params, ...
steadystate, sparse_rowval, sparse_colval, ...
sparse_colptr, T);
end
endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
endogenousvariables(:,it) = y_dynamic(M_.endo_nbr+(1:M_.endo_nbr));
else
[tmp, check, ~, ~, errorcode] = dynare_solve(@dynamic_static_model_for_simulation, y, ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, dynamic_resid, dynamic_g1, x, M.params, steadystate, M.dynamic_g1_sparse_rowval, M.dynamic_g1_sparse_colval, M.dynamic_g1_sparse_colptr);
options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
options_, dynamic_resid, dynamic_g1, x, M_.params, steadystate, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, M_.dynamic_g1_sparse_colptr);
if check
success = false;
if options.debug
if options_.debug
dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
end
end

38
matlab/perfect-foresight-models/solve_stacked_linear_problem.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, success] = solve_stacked_linear_problem(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
function [endogenousvariables, success] = solve_stacked_linear_problem(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M_, options_)
% Copyright © 2015-2023 Dynare Team
%
@ -17,40 +17,40 @@ function [endogenousvariables, success] = solve_stacked_linear_problem(endogenou
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i_cols_0, i_cols_J0, dynamicmodel] = ...
initialize_stacked_problem(endogenousvariables, options, M, steadystate_y);
[options_, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i_cols_0, i_cols_J0, dynamicmodel] = ...
initialize_stacked_problem(endogenousvariables, options_, M_, steadystate_y);
ip = find(M.lead_lag_incidence(1,:)');
ic = find(M.lead_lag_incidence(2,:)');
in = find(M.lead_lag_incidence(3,:)');
ip = find(M_.lead_lag_incidence(1,:)');
ic = find(M_.lead_lag_incidence(2,:)');
in = find(M_.lead_lag_incidence(3,:)');
% Evaluate the Jacobian of the dynamic model at the deterministic steady state.
[d1,jacobian] = dynamicmodel(steadystate_y([ip; ic; in]), transpose(steadystate_x), M.params, steadystate_y, 1);
[d1,jacobian] = dynamicmodel(steadystate_y([ip; ic; in]), transpose(steadystate_x), M_.params, steadystate_y, 1);
% Check that the dynamic model was evaluated at the steady state.
if ~options.steadystate.nocheck && max(abs(d1))>1e-12
if ~options_.steadystate.nocheck && max(abs(d1))>1e-12
error('Jacobian is not evaluated at the steady state!')
end
nyp = nnz(M.lead_lag_incidence(1,:));
ny0 = nnz(M.lead_lag_incidence(2,:));
nyf = nnz(M.lead_lag_incidence(3,:));
nyp = nnz(M_.lead_lag_incidence(1,:));
ny0 = nnz(M_.lead_lag_incidence(2,:));
nyf = nnz(M_.lead_lag_incidence(3,:));
nd = nyp+ny0+nyf; % size of y (first argument passed to the dynamic file).
jexog = transpose(nd+(1:M.exo_nbr));
jexog = transpose(nd+(1:M_.exo_nbr));
jendo = transpose(1:nd);
z = bsxfun(@minus, z, steadystate_y);
x = bsxfun(@minus, exogenousvariables, steadystate_x');
[y, check, ~, ~, errorcode] = dynare_solve(@linear_perfect_foresight_problem, z(:), ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, ...
options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
options_, ...
jacobian, y0-steadystate_y, yT-steadystate_y, ...
x, M.params, steadystate_y, ...
M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
x, M_.params, steadystate_y, ...
M_.maximum_lag, options_.periods, M_.endo_nbr, i_cols, ...
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
jendo, jexog);
if all(imag(y)<.1*options.dynatol.x)
if all(imag(y)<.1*options_.dynatol.x)
if ~isreal(y)
y = real(y);
end
@ -58,10 +58,10 @@ else
check = 1;
end
endogenousvariables = [y0 bsxfun(@plus,reshape(y,M.endo_nbr,options.periods), steadystate_y) yT];
endogenousvariables = [y0 bsxfun(@plus,reshape(y,M_.endo_nbr,options_.periods), steadystate_y) yT];
success = ~check;
if ~success && options.debug
if ~success && options_.debug
dprintf('solve_stacked_linear_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
end

52
matlab/perfect-foresight-models/solve_stacked_problem.m
View File

@ -1,13 +1,13 @@
function [endogenousvariables, success, maxerror] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success, maxerror] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M_, options_)
% [endogenousvariables, success, maxerror] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M_, options_)
% Solves the perfect foresight model using dynare_solve
%
% INPUTS
% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial guess).
% - exogenousvariables [double] T*M array, paths for the exogenous variables.
% - steadystate [double] N*1 array, steady state for the endogenous variables.
% - M [struct] contains a description of the model.
% - options [struct] contains various options.
% - M_ [struct] contains a description of the model.
% - options_ [struct] contains various options.
%
% OUTPUTS
% - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
@ -31,39 +31,39 @@ function [endogenousvariables, success, maxerror] = solve_stacked_problem(endoge
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i_cols_0, i_cols_J0, dynamicmodel] = ...
initialize_stacked_problem(endogenousvariables, options, M, steadystate);
[options_, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i_cols_0, i_cols_J0, dynamicmodel] = ...
initialize_stacked_problem(endogenousvariables, options_, M_, steadystate);
if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity problem
[lb,ub,eq_index] = get_complementarity_conditions(M,options.ramsey_policy);
if options.linear_approximation
if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementarity problem
[lb,ub,eq_index] = get_complementarity_conditions(M_,options_.ramsey_policy);
if options_.linear_approximation
lb = lb - steadystate_y;
ub = ub - steadystate_y;
end
if options.solve_algo == 10
options.lmmcp.lb = repmat(lb,options.periods,1);
options.lmmcp.ub = repmat(ub,options.periods,1);
elseif options.solve_algo == 11
options.mcppath.lb = repmat(lb,options.periods,1);
options.mcppath.ub = repmat(ub,options.periods,1);
if options_.solve_algo == 10
options_.lmmcp.lb = repmat(lb,options_.periods,1);
options_.lmmcp.ub = repmat(ub,options_.periods,1);
elseif options_.solve_algo == 11
options_.mcppath.lb = repmat(lb,options_.periods,1);
options_.mcppath.ub = repmat(ub,options_.periods,1);
end
[y, check, res, ~, errorcode] = dynare_solve(@perfect_foresight_mcp_problem, z(:), ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, ...
options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
options_, ...
dynamicmodel, y0, yT, ...
exogenousvariables, M.params, steadystate, ...
M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
exogenousvariables, M_.params, steadystate, ...
M_.maximum_lag, options_.periods, M_.endo_nbr, i_cols, ...
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
eq_index);
eq_to_ignore=find(isfinite(lb) | isfinite(ub));
else
[y, check, res, ~, errorcode] = dynare_solve(@perfect_foresight_problem, z(:), ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, y0, yT, exogenousvariables, M.params, steadystate, options.periods, M, options);
options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
options_, y0, yT, exogenousvariables, M_.params, steadystate, options_.periods, M_, options_);
end
if all(imag(y)<.1*options.dynatol.x)
if all(imag(y)<.1*options_.dynatol.x)
if ~isreal(y)
y = real(y);
end
@ -71,15 +71,15 @@ else
check = 1;
end
endogenousvariables(:, M.maximum_lag+(1:options.periods)) = reshape(y, M.endo_nbr, options.periods);
endogenousvariables(:, M_.maximum_lag+(1:options_.periods)) = reshape(y, M_.endo_nbr, options_.periods);
residuals=zeros(size(endogenousvariables));
residuals(:, M.maximum_lag+(1:options.periods)) = reshape(res, M.endo_nbr, options.periods);
if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity problem
residuals(:, M_.maximum_lag+(1:options_.periods)) = reshape(res, M_.endo_nbr, options_.periods);
if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementarity problem
residuals(eq_to_ignore,bsxfun(@le, endogenousvariables(eq_to_ignore,:), lb(eq_to_ignore)+eps) | bsxfun(@ge,endogenousvariables(eq_to_ignore,:),ub(eq_to_ignore)-eps))=0;
end
maxerror = max(max(abs(residuals)));
success = ~check;
if ~success && options.debug
if ~success && options_.debug
dprintf('solve_stacked_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
end

62
matlab/pruned_state_space_system.m
View File

@ -1,4 +1,4 @@
function pruned_state_space = pruned_state_space_system(M, options, dr, indy, nlags, useautocorr, compute_derivs)
function pruned_state_space = pruned_state_space_system(M_, options_, dr, indy, nlags, useautocorr, compute_derivs)
% Set up the pruned state space ABCD representation:
% z = c + A*z(-1) + B*inov
% y = ys + d + C*z(-1) + D*inov
@ -10,8 +10,8 @@ function pruned_state_space = pruned_state_space_system(M, options, dr, indy, nl
% Econometrics and Statistics, Volume 6, Pages 44-56.
% =========================================================================
% INPUTS
% M: [structure] storing the model information
% options: [structure] storing the options
% M_: [structure] storing the model information
% options_: [structure] storing the options
% dr: [structure] storing the results from perturbation approximation
% indy: [vector] index of control variables in DR order
% nlags: [integer] number of lags in autocovariances and autocorrelations
@ -244,19 +244,19 @@ function pruned_state_space = pruned_state_space_system(M, options, dr, indy, nl
persistent QPu COMBOS4 Q6Pu COMBOS6 K_u_xx K_u_ux K_xx_x
%% Auxiliary indices and objects
order = options.order;
if isempty(options.qz_criterium)
order = options_.order;
if isempty(options_.qz_criterium)
% set default value for qz_criterium: if there are no unit roots one can use 1.0
% If they are possible, you may have have multiple unit roots and the accuracy
% decreases when computing the eigenvalues in lyapunov_symm. Hence, we normally use 1+1e-6
% Note that unit roots are only possible at first-order, at higher order we set it to 1
options.qz_criterium = 1+1e-6;
options_.qz_criterium = 1+1e-6;
end
indx = [M.nstatic+(1:M.nspred) M.endo_nbr+(1:size(dr.ghx,2)-M.nspred)]';
indx = [M_.nstatic+(1:M_.nspred) M_.endo_nbr+(1:size(dr.ghx,2)-M_.nspred)]';
if isempty(indy)
indy = (1:M.endo_nbr)'; %by default select all variables in DR order
indy = (1:M_.endo_nbr)'; %by default select all variables in DR order
end
u_nbr = M.exo_nbr;
u_nbr = M_.exo_nbr;
x_nbr = length(indx);
y_nbr = length(indy);
Yss = dr.ys(dr.order_var);
@ -264,7 +264,7 @@ hx = dr.ghx(indx,:);
gx = dr.ghx(indy,:);
hu = dr.ghu(indx,:);
gu = dr.ghu(indy,:);
E_uu = M.Sigma_e; %this is E[u*u']
E_uu = M_.Sigma_e; %this is E[u*u']
if compute_derivs
stderrparam_nbr = length(dr.derivs.indpstderr);
@ -350,7 +350,7 @@ end
% Auxiliary state vector z is defined by: z = [xf]
% Auxiliary innovations vector inov is defined by: inov = [u]
z_nbr = x_nbr;
inov_nbr = M.exo_nbr;
inov_nbr = M_.exo_nbr;
A = hx;
B = hu;
c = zeros(x_nbr,1);
@ -364,9 +364,9 @@ E_xf = zeros(x_nbr,1);
lyapunov_symm_method = 1; %method=1 to initialize persistent variables
[Var_z,Schur_u] = lyapunov_symm(A, Om_z,... %at first-order this algorithm is well established and also used in th_autocovariances.m
options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold,...
options_.lyapunov_fixed_point_tol, options_.qz_criterium, options_.lyapunov_complex_threshold,...
lyapunov_symm_method,...
options.debug); %we use Schur_u to take care of (possible) nonstationary VAROBS variables in moment computations
options_.debug); %we use Schur_u to take care of (possible) nonstationary VAROBS variables in moment computations
%find stationary vars
stationary_vars = (1:y_nbr)';
if ~isempty(Schur_u)
@ -375,7 +375,7 @@ if ~isempty(Schur_u)
end
%base this only on first order, because if first-order is stable so are the higher-order pruned systems
x = abs(gx*Schur_u);
stationary_vars = find(all(x < options.schur_vec_tol,2));
stationary_vars = find(all(x < options_.schur_vec_tol,2));
end
if compute_derivs == 1
@ -394,21 +394,21 @@ if compute_derivs == 1
if jp1 <= stderrparam_nbr+corrparam_nbr
dOm_z_jp1 = B*dVarinov(:,:,jp1)*B';
dVar_z(:,:,jp1) = lyapunov_symm(A, dOm_z_jp1,...
options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,...
options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,...
lyapunov_symm_method,...
options.debug);
options_.debug);
else
dOm_z_jp1 = dB(:,:,jp1)*Varinov*B' + B*Varinov*dB(:,:,jp1)';
dVar_z(:,:,jp1) = lyapunov_symm(A, dA(:,:,jp1)*Var_z*A' + A*Var_z*dA(:,:,jp1)' + dOm_z_jp1,...
options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,...
options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,...
lyapunov_symm_method,...
options.debug);
options_.debug);
end
end
end
if order > 1
options.qz_criterium = 1; %pruned state space has no unit roots
options_.qz_criterium = 1; %pruned state space has no unit roots
% Some common and useful objects for order > 1
E_xfxf = Var_z;
if compute_derivs
@ -491,13 +491,13 @@ if order > 1
Om_z = B*Varinov*transpose(B);
lyapunov_symm_method = 1; %method=1 to initialize persistent variables (if errorflag)
[Var_z, errorflag] = disclyap_fast(A,Om_z,options.lyapunov_doubling_tol);
[Var_z, errorflag] = disclyap_fast(A,Om_z,options_.lyapunov_doubling_tol);
if errorflag %use Schur-based method
fprintf('PRUNED_STATE_SPACE_SYSTEM: error flag in disclyap_fast at order=2, use lyapunov_symm\n');
Var_z = lyapunov_symm(A,Om_z,...
options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,...
options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,...
lyapunov_symm_method,...
options.debug);
options_.debug);
lyapunov_symm_method = 2; %in the following we can reuse persistent variables
end
% Make sure some stuff is zero due to Gaussianity
@ -575,12 +575,12 @@ if order > 1
dE_xs(:,jp2) = invIx_hx*( dhx_jp2*E_xs + 1/2*dhxx_jp2*E_xfxf(:) + 1/2*hxx*dE_xfxf_jp2(:) + dc(id_z2_xs,jp2) );
dOm_z_jp2 = dB(:,:,jp2)*Varinov*B' + B*dVarinov(:,:,jp2)*B' + B*Varinov*dB(:,:,jp2)';
[dVar_z(:,:,jp2), errorflag] = disclyap_fast(A, dA(:,:,jp2)*Var_z*A' + A*Var_z*dA(:,:,jp2)' + dOm_z_jp2, options.lyapunov_doubling_tol);
[dVar_z(:,:,jp2), errorflag] = disclyap_fast(A, dA(:,:,jp2)*Var_z*A' + A*Var_z*dA(:,:,jp2)' + dOm_z_jp2, options_.lyapunov_doubling_tol);
if errorflag
dVar_z(:,:,jp2) = lyapunov_symm(A, dA(:,:,jp2)*Var_z*A' + A*Var_z*dA(:,:,jp2)' + dOm_z_jp2,...
options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,...
options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,...
lyapunov_symm_method,...
options.debug);
options_.debug);
if lyapunov_symm_method == 1
lyapunov_symm_method = 2; %now we can reuse persistent schur
end
@ -817,13 +817,13 @@ if order > 1
Om_z = B*Varinov*transpose(B) + Binovzlag1A + transpose(Binovzlag1A);
lyapunov_symm_method = 1; %method=1 to initialize persistent variables
[Var_z, errorflag] = disclyap_fast(A,Om_z,options.lyapunov_doubling_tol);
[Var_z, errorflag] = disclyap_fast(A,Om_z,options_.lyapunov_doubling_tol);
if errorflag %use Schur-based method
fprintf('PRUNED_STATE_SPACE_SYSTEM: error flag in disclyap_fast at order=3, use lyapunov_symm\n');
Var_z = lyapunov_symm(A,Om_z,...
options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,...
options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,...
lyapunov_symm_method,...
options.debug);
options_.debug);
lyapunov_symm_method = 2; %we can now make use of persistent variables from shur
end
%make sure some stuff is zero due to Gaussianity
@ -998,12 +998,12 @@ if order > 1
dBinovzlag1A_jp3 = dB(:,:,jp3)*E_inovzlag1*transpose(A) + B*dE_inovzlag1(:,:,jp3)*transpose(A) + B*E_inovzlag1*transpose(dA(:,:,jp3));
dOm_z_jp3 = dB(:,:,jp3)*Varinov*transpose(B) + B*dVarinov(:,:,jp3)*transpose(B) + B*Varinov*transpose(dB(:,:,jp3)) + dBinovzlag1A_jp3 + transpose(dBinovzlag1A_jp3);
[dVar_z(:,:,jp3), errorflag] = disclyap_fast(A, dA(:,:,jp3)*Var_z*A' + A*Var_z*dA(:,:,jp3)' + dOm_z_jp3, options.lyapunov_doubling_tol);
[dVar_z(:,:,jp3), errorflag] = disclyap_fast(A, dA(:,:,jp3)*Var_z*A' + A*Var_z*dA(:,:,jp3)' + dOm_z_jp3, options_.lyapunov_doubling_tol);
if errorflag
dVar_z(:,:,jp3) = lyapunov_symm(A, dA(:,:,jp3)*Var_z*A' + A*Var_z*dA(:,:,jp3)' + dOm_z_jp3,...
options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,...
options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,...
lyapunov_symm_method,...
options.debug);
options_.debug);
if lyapunov_symm_method == 1
lyapunov_symm_method = 2; %now we can reuse persistent schur
end

56
matlab/set_all_parameters.m
View File

@ -1,7 +1,7 @@
function M = set_all_parameters(xparam1,estim_params,M)
function M_ = set_all_parameters(xparam1,estim_params_,M_)
%@info:
%! @deftypefn {Function File} {@var{M} =} dseries (@var{xparams1},@var{estim_params},@var{M})
%! @deftypefn {Function File} {@var{M_} =} dseries (@var{xparams1},@var{estim_params_},@var{M_})
%! @anchor{set_all_parameters}
%! @sp 1
%! Update parameter values (deep parameters and covariance matrices).
@ -11,16 +11,16 @@ function M = set_all_parameters(xparam1,estim_params,M)
%! @table @ @var
%! @item xparam1
%! N*1 vector of doubles, the values of the N estimated parameters.
%! @item estim_params
%! @item estim_params_
%! Dynare structure describing the estimated parameters.
%! @item M
%! @item M_
%! Dynare structure describing the model.
%! @end table
%! @sp 1
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item M
%! @item M_
%! Dynare structure describing the model, with updated parameters and covariances matrices.
%! @end table
%! @sp 2
@ -50,21 +50,21 @@ function M = set_all_parameters(xparam1,estim_params,M)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
nvx = estim_params.nvx;
ncx = estim_params.ncx;
nvn = estim_params.nvn;
ncn = estim_params.ncn;
np = estim_params.np;
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
np = estim_params_.np;
if nvx || ncx
Sigma_e = M.Sigma_e;
Correlation_matrix = M.Correlation_matrix;
Sigma_e = M_.Sigma_e;
Correlation_matrix = M_.Correlation_matrix;
end
H = M.H;
Correlation_matrix_ME = M.Correlation_matrix_ME;
H = M_.H;
Correlation_matrix_ME = M_.Correlation_matrix_ME;
% setting shocks variance on the diagonal of Covariance matrix; used later
% for updating covariances
if nvx
var_exo = estim_params.var_exo;
var_exo = estim_params_.var_exo;
for i=1:nvx
k =var_exo(i,1);
Sigma_e(k,k) = xparam1(i)^2;
@ -77,7 +77,7 @@ offset = nvx;
% for updating covariances
if nvn
for i=1:nvn
k = estim_params.nvn_observable_correspondence(i,1);
k = estim_params_.nvn_observable_correspondence(i,1);
H(k,k) = xparam1(i+offset)^2;
end
end
@ -87,7 +87,7 @@ offset = nvx+nvn;
% setting shocks covariances
if ncx
corrx = estim_params.corrx;
corrx = estim_params_.corrx;
for i=1:ncx
k1 = corrx(i,1);
k2 = corrx(i,2);
@ -100,7 +100,7 @@ offset = nvx+nvn+ncx;
% setting measurement error covariances
if ncn
corrn_observable_correspondence = estim_params.corrn_observable_correspondence;
corrn_observable_correspondence = estim_params_.corrn_observable_correspondence;
for i=1:ncn
k1 = corrn_observable_correspondence(i,1);
k2 = corrn_observable_correspondence(i,2);
@ -114,29 +114,29 @@ offset = nvx+ncx+nvn+ncn;
% setting structural parameters
%
if np
M.params(estim_params.param_vals(:,1)) = xparam1(offset+1:end);
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end
% updating matrices in M
% updating matrices in M_
if nvx || ncx
%build covariance matrix from correlation matrix and variances already on
%diagonal
Sigma_e = diag(sqrt(diag(Sigma_e)))*Correlation_matrix*diag(sqrt(diag(Sigma_e)));
%if calibrated covariances, set them now to their stored value
if isfield(estim_params,'calibrated_covariances')
Sigma_e(estim_params.calibrated_covariances.position)=estim_params.calibrated_covariances.cov_value;
if isfield(estim_params_,'calibrated_covariances')
Sigma_e(estim_params_.calibrated_covariances.position)=estim_params_.calibrated_covariances.cov_value;
end
M.Sigma_e = Sigma_e;
M.Correlation_matrix=Correlation_matrix;
M_.Sigma_e = Sigma_e;
M_.Correlation_matrix=Correlation_matrix;
end
if nvn || ncn
%build covariance matrix from correlation matrix and variances already on
%diagonal
H = diag(sqrt(diag(H)))*Correlation_matrix_ME*diag(sqrt(diag(H)));
%if calibrated covariances, set them now to their stored value
if isfield(estim_params,'calibrated_covariances_ME')
H(estim_params.calibrated_covariances_ME.position)=estim_params.calibrated_covariances_ME.cov_value;
if isfield(estim_params_,'calibrated_covariances_ME')
H(estim_params_.calibrated_covariances_ME.position)=estim_params_.calibrated_covariances_ME.cov_value;
end
M.H = H;
M.Correlation_matrix_ME=Correlation_matrix_ME;
M_.H = H;
M_.Correlation_matrix_ME=Correlation_matrix_ME;
end

9
matlab/set_local_param_value.m
View File

@ -1,6 +1,7 @@
function M = set_local_param_value(pname,value,M)
function M_ = set_local_param_value(pname,value,M_)
% M_ = set_local_param_value(pname,value,M_)
% Copyright © 2017-2018 Dynare Team
% Copyright © 2017-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -18,10 +19,10 @@ function M = set_local_param_value(pname,value,M)
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
i = strmatch(pname, M.param_names, 'exact');
i = strmatch(pname, M_.param_names, 'exact');
if isempty(i)
error(['Parameter name ' pname ' doesn''t exist'])
end
M.params(i) = value;
M_.params(i) = value;

60
matlab/solve_one_boundary.m
View File

@ -1,5 +1,5 @@
function [y, T, success, max_res, iter] = solve_one_boundary(fh, y, x, params, steady_state, T, ...
y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo, is_forward, is_dynamic, verbose, M, options)
y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo, is_forward, is_dynamic, verbose, M_, options_)
% Computes the deterministic simulation or the steady state for a block of equations containing
% only lags or only leads (but not both).
%
@ -27,6 +27,8 @@ function [y, T, success, max_res, iter] = solve_one_boundary(fh, y, x, params, s
% If false, the block is solved backward
% is_dynamic [logical] Whether this is a deterministic simulation
% verbose [logical] Whether iterations are to be printed
% M_ [structure] storing the model information
% options_ [structure] storing the options
%
% OUTPUTS
% y [matrix] All endogenous variables of the model
@ -78,15 +80,15 @@ for it_=start:incr:finish
g1=spalloc( Blck_size, Blck_size, nze);
while ~(cvg || iter>maxit_)
if is_dynamic
[yy, T(:, it_), r, g1] = fh(dynendo(y, it_, M), x(it_, :), params, steady_state, ...
M.block_structure.block(Block_Num).g1_sparse_rowval, ...
M.block_structure.block(Block_Num).g1_sparse_colval, ...
M.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
y(:, it_) = yy(M.endo_nbr+(1:M.endo_nbr));
[yy, T(:, it_), r, g1] = fh(dynendo(y, it_, M_), x(it_, :), params, steady_state, ...
M_.block_structure.block(Block_Num).g1_sparse_rowval, ...
M_.block_structure.block(Block_Num).g1_sparse_colval, ...
M_.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
y(:, it_) = yy(M_.endo_nbr+(1:M_.endo_nbr));
else
[y, T, r, g1] = fh(y, x, params, M.block_structure_stat.block(Block_Num).g1_sparse_rowval, ...
M.block_structure_stat.block(Block_Num).g1_sparse_colval, ...
M.block_structure_stat.block(Block_Num).g1_sparse_colptr, T);
[y, T, r, g1] = fh(y, x, params, M_.block_structure_stat.block(Block_Num).g1_sparse_rowval, ...
M_.block_structure_stat.block(Block_Num).g1_sparse_colval, ...
M_.block_structure_stat.block(Block_Num).g1_sparse_colptr, T);
end
if ~isreal(r)
max_res=(-(max(max(abs(r))))^2)^0.5;
@ -96,9 +98,9 @@ for it_=start:incr:finish
if verbose
disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)])
if is_dynamic
disp([char(M.endo_names{y_index_eq}) repmat(' ', numel(y_index_eq), 1) num2str([y(y_index_eq, it_) r g1])])
disp([char(M_.endo_names{y_index_eq}) repmat(' ', numel(y_index_eq), 1) num2str([y(y_index_eq, it_) r g1])])
else
disp([char(M.endo_names{y_index_eq}) repmat(' ', numel(y_index_eq), 1) num2str([y(y_index_eq) r g1])])
disp([char(M_.endo_names{y_index_eq}) repmat(' ', numel(y_index_eq), 1) num2str([y(y_index_eq) r g1])])
end
end
if ~isreal(max_res) || isnan(max_res)
@ -182,15 +184,15 @@ for it_=start:incr:finish
p = -g1\r ;
[ya,f,r,check]=lnsrch1(ya,f,g,p,stpmax, ...
@lnsrch1_wrapper_one_boundary,nn, ...
nn, options.solve_tolx, y_index_eq, fh, Block_Num, y, x, params, steady_state, T(:, it_), it_, M);
nn, options_.solve_tolx, y_index_eq, fh, Block_Num, y, x, params, steady_state, T(:, it_), it_, M_);
dx = ya - y(y_index_eq, it_);
y(y_index_eq, it_) = ya;
%% Recompute temporary terms, since they are not given as output of lnsrch1
[~, T(:, it_)] = fh(dynendo(y, it_, M), x(it_, :), params, steady_state, ...
M.block_structure.block(Block_Num).g1_sparse_rowval, ...
M.block_structure.block(Block_Num).g1_sparse_colval, ...
M.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0 || stack_solve_algo==6)) || (~is_dynamic && options.solve_algo==6)
[~, T(:, it_)] = fh(dynendo(y, it_, M_), x(it_, :), params, steady_state, ...
M_.block_structure.block(Block_Num).g1_sparse_rowval, ...
M_.block_structure.block(Block_Num).g1_sparse_colval, ...
M_.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0 || stack_solve_algo==6)) || (~is_dynamic && options_.solve_algo==6)
if verbose && ~is_dynamic
disp('steady: Sparse LU ')
end
@ -201,7 +203,7 @@ for it_=start:incr:finish
else
y(y_index_eq) = ya;
end
elseif (stack_solve_algo==2 && is_dynamic) || (options.solve_algo==7 && ~is_dynamic)
elseif (stack_solve_algo==2 && is_dynamic) || (options_.solve_algo==7 && ~is_dynamic)
flag1=1;
if verbose && ~is_dynamic
disp('steady: GMRES ')
@ -230,7 +232,7 @@ for it_=start:incr:finish
end
end
end
elseif (stack_solve_algo==3 && is_dynamic) || (options.solve_algo==8 && ~is_dynamic)
elseif (stack_solve_algo==3 && is_dynamic) || (options_.solve_algo==8 && ~is_dynamic)
flag1=1;
if verbose && ~is_dynamic
disp('steady: BiCGStab')
@ -263,7 +265,7 @@ for it_=start:incr:finish
if is_dynamic
error(['options_.stack_solve_algo = ' num2str(stack_solve_algo) ' not implemented'])
else
error(['options_.solve_algo = ' num2str(options.solve_algo) ' not implemented'])
error(['options_.solve_algo = ' num2str(options_.solve_algo) ' not implemented'])
end
end
iter=iter+1;
@ -286,20 +288,20 @@ end
success = true;
function y3n = dynendo(y, it_, M)
function y3n = dynendo(y, it_, M_)
if it_ > 1 && it_ < size(y, 2)
y3n = reshape(y(:, it_+(-1:1)), 3*M.endo_nbr, 1);
y3n = reshape(y(:, it_+(-1:1)), 3*M_.endo_nbr, 1);
elseif it_ > 1 % Purely backward model (in last period)
y3n = [ reshape(y(:, it_+(-1:0)), 2*M.endo_nbr, 1); NaN(M_.endo_nbr, 1) ];
y3n = [ reshape(y(:, it_+(-1:0)), 2*M_.endo_nbr, 1); NaN(M_.endo_nbr, 1) ];
elseif it_ < size(y, 2) % Purely forward model (in first period)
y3n = [ NaN(M_.endo_nbr, 1); reshape(y(:, it_+(0:1)), 2*M.endo_nbr, 1) ];
y3n = [ NaN(M_.endo_nbr, 1); reshape(y(:, it_+(0:1)), 2*M_.endo_nbr, 1) ];
else % Static model
y3n = [ NaN(M_.endo_nbr, 1); y(:, it_); NaN(M_.endo_nbr, 1) ]
end
function r = lnsrch1_wrapper_one_boundary(ya, y_index, fh, Block_Num, y, x, params, steady_state, T, it_, M)
function r = lnsrch1_wrapper_one_boundary(ya, y_index, fh, Block_Num, y, x, params, steady_state, T, it_, M_)
y(y_index, it_) = ya;
[~, ~, r] = fh(dynendo(y, it_, M), x(it_, :), params, steady_state, ...
M.block_structure.block(Block_Num).g1_sparse_rowval, ...
M.block_structure.block(Block_Num).g1_sparse_colval, ...
M.block_structure.block(Block_Num).g1_sparse_colptr, T);
[~, ~, r] = fh(dynendo(y, it_, M_), x(it_, :), params, steady_state, ...
M_.block_structure.block(Block_Num).g1_sparse_rowval, ...
M_.block_structure.block(Block_Num).g1_sparse_colval, ...
M_.block_structure.block(Block_Num).g1_sparse_colptr, T);

33
matlab/solve_two_boundaries.m
View File

@ -1,4 +1,4 @@
function [y, T, success, max_res, iter] = solve_two_boundaries(fh, y, x, params, steady_state, T, y_index, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo,options,M)
function [y, T, success, max_res, iter] = solve_two_boundaries(fh, y, x, params, steady_state, T, y_index, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo,options_,M_)
% Computes the deterministic simulation of a block of equation containing
% both lead and lag variables using relaxation methods
%
@ -27,7 +27,8 @@ function [y, T, success, max_res, iter] = solve_two_boundaries(fh, y, x, params,
% - 2 GMRES
% - 3 BicGStab
% - 4 Optimal path length
% M [structure] Model description
% options_ [structure] storing the options
% M_ [structure] Model description
%
% OUTPUTS
% y [matrix] All endogenous variables of the model
@ -56,7 +57,7 @@ function [y, T, success, max_res, iter] = solve_two_boundaries(fh, y, x, params,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
verbose = options.verbosity;
verbose = options_.verbosity;
cvg=false;
iter=0;
@ -76,11 +77,11 @@ while ~(cvg || iter>maxit_)
r = NaN(Blck_size, periods);
g1a = spalloc(Blck_size*periods, Blck_size*periods, nze*periods);
for it_ = y_kmin+(1:periods)
[yy, T(:, it_), r(:, it_-y_kmin), g1]=fh(dynendo(y, it_, M), x(it_, :), params, steady_state, ...
M.block_structure.block(Block_Num).g1_sparse_rowval, ...
M.block_structure.block(Block_Num).g1_sparse_colval, ...
M.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
y(:, it_) = yy(M.endo_nbr+(1:M.endo_nbr));
[yy, T(:, it_), r(:, it_-y_kmin), g1]=fh(dynendo(y, it_, M_), x(it_, :), params, steady_state, ...
M_.block_structure.block(Block_Num).g1_sparse_rowval, ...
M_.block_structure.block(Block_Num).g1_sparse_colval, ...
M_.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
y(:, it_) = yy(M_.endo_nbr+(1:M_.endo_nbr));
if periods == 1
g1a = g1(:, Blck_size+(1:Blck_size));
elseif it_ == y_kmin+1
@ -110,7 +111,7 @@ while ~(cvg || iter>maxit_)
if iter>0
if ~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1)
if verbose && ~isreal(max_res)
disp(['Variable ' M.endo_names{max_indx} ' (' int2str(max_indx) ') returns an undefined value']);
disp(['Variable ' M_.endo_names{max_indx} ' (' int2str(max_indx) ') returns an undefined value']);
end
if isnan(max_res)
detJ=det(g1aa);
@ -308,7 +309,7 @@ while ~(cvg || iter>maxit_)
g = (ra'*g1a)';
f = 0.5*ra'*ra;
p = -g1a\ra;
[yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,@lnsrch1_wrapper_two_boundaries,nn,nn, options.solve_tolx, fh, Block_Num, y, y_index,x, params, steady_state, T, periods, Blck_size, M);
[yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,@lnsrch1_wrapper_two_boundaries,nn,nn, options_.solve_tolx, fh, Block_Num, y, y_index,x, params, steady_state, T, periods, Blck_size, M_);
dx = ya - yn;
y(y_index, y_kmin+(1:periods))=reshape(yn',length(y_index),periods);
end
@ -331,14 +332,14 @@ end
success = true;
function y3n = dynendo(y, it_, M)
y3n = reshape(y(:, it_+(-1:1)), 3*M.endo_nbr, 1);
function y3n = dynendo(y, it_, M_)
y3n = reshape(y(:, it_+(-1:1)), 3*M_.endo_nbr, 1);
function ra = lnsrch1_wrapper_two_boundaries(ya, fh, Block_Num, y, y_index, x, ...
params, steady_state, T, periods, ...
y_size, M)
y(y_index, M.maximum_lag+(1:periods)) = reshape(ya',length(y_index),periods);
y_size, M_)
y(y_index, M_.maximum_lag+(1:periods)) = reshape(ya',length(y_index),periods);
ra = NaN(periods*y_size, 1);
for it_ = M.maximum_lag+(1:periods)
[~, ~, ra((it_-M.maximum_lag-1)*y_size+(1:y_size)), g1] = fh(dynendo(y, it_, M), x(it_, :), params, steady_state, M.block_structure.block(Block_Num).g1_sparse_rowval, M.block_structure.block(Block_Num).g1_sparse_colval, M.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
for it_ = M_.maximum_lag+(1:periods)
[~, ~, ra((it_-M_.maximum_lag-1)*y_size+(1:y_size)), g1] = fh(dynendo(y, it_, M_), x(it_, :), params, steady_state, M_.block_structure.block(Block_Num).g1_sparse_rowval, M_.block_structure.block(Block_Num).g1_sparse_colval, M_.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
end

15
matlab/utilities/general/get_all_variables_but_lagged_leaded_exogenous.m
View File

@ -1,15 +1,16 @@
function k = get_all_variables_but_lagged_leaded_exogenous(M)
function k = get_all_variables_but_lagged_leaded_exogenous(M_)
% k = get_all_variables_but_lagged_leaded_exogenous(M_)
% returns indices of all endogenous variables in declaration order except
% lagged and leaded exogenous
%
% INPUT
% M: model structure
% M_: model structure
%
% OUTPUT
% k: vector of variable indices
%
% Copyright © 2011-2017 Dynare Team
% Copyright © 2011-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -26,10 +27,10 @@ function k = get_all_variables_but_lagged_leaded_exogenous(M)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if isempty(M.aux_vars)
k = 1:M.endo_nbr;
if isempty(M_.aux_vars)
k = 1:M_.endo_nbr;
else
type = [M.aux_vars.type];
k = [1:M.orig_endo_nbr, M.orig_endo_nbr ...
type = [M_.aux_vars.type];
k = [1:M_.orig_endo_nbr, M_.orig_endo_nbr ...
+ find((type ~= 2) & (type ~= 3))];
end

11
matlab/utilities/general/indices_lagged_leaded_exogenous_variables.m
View File

@ -1,10 +1,11 @@
function [il,l1,ik,k1] = indices_lagged_leaded_exogenous_variables(k,M)
function [il,l1,ik,k1] = indices_lagged_leaded_exogenous_variables(k,M_)
% [il,l1,ik,k1] = indices_lagged_leaded_exogenous_variables(k,M_)
% returns indices of all endogenous variables split between auxiliary
% variables for lagged or leaded exogenous variables and all other ones
%
% INPUT
% k: vector of endogenous variables ID
% M: model structure
% M_: model structure
%
% OUTPUT
% il: indices of lagged or leaded variable in vector k
@ -31,14 +32,14 @@ function [il,l1,ik,k1] = indices_lagged_leaded_exogenous_variables(k,M)
il = [];
l1 = [];
if isempty(M.aux_vars)
if isempty(M_.aux_vars)
ik = 1:length(k);
k1 = k;
else
ik = [];
k1 = [];
orig_endo_nbr = M.orig_endo_nbr;
type = [M.aux_vars.type];
orig_endo_nbr = M_.orig_endo_nbr;
type = [M_.aux_vars.type];
for j=1:length(k)
if (k(j) > orig_endo_nbr)
ty = type(k(j) - orig_endo_nbr);

8
matlab/utilities/general/log_variable.m
View File

@ -1,11 +1,11 @@
function y = log_variable(ivar,x,M)
function y = log_variable(ivar,x,M_)
% returns the log of an endogenous variables except if
% it is a lagged/leaded exogenous variable
%
% INPUT
% ivar: vector of variable indices (in declaration order)
% x: vector of variables value
% M: model structure
% M_: model structure
%
% OUTPUT
% y: log of the selected variables if there are not auxiliary variables
@ -29,8 +29,8 @@ function y = log_variable(ivar,x,M)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
orig_endo_nbr = M.orig_endo_nbr;
aux_vars = M.aux_vars;
orig_endo_nbr = M_.orig_endo_nbr;
aux_vars = M_.aux_vars;
y = zeros(length(ivar),1);
for i = 1:length(ivar)