diff --git a/matlab/AIM_first_order_solver.m b/matlab/AIM_first_order_solver.m
index 026c2619e..ce149e1a6 100644
--- a/matlab/AIM_first_order_solver.m
+++ b/matlab/AIM_first_order_solver.m
@@ -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));
diff --git a/matlab/check_for_calibrated_covariances.m b/matlab/check_for_calibrated_covariances.m
index 276a1c932..e07b89dba 100644
--- a/matlab/check_for_calibrated_covariances.m
+++ b/matlab/check_for_calibrated_covariances.m
@@ -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 .
-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;
diff --git a/matlab/dyn_ramsey_static.m b/matlab/dyn_ramsey_static.m
index 32d399b1d..ac4b4d23c 100644
--- a/matlab/dyn_ramsey_static.m
+++ b/matlab/dyn_ramsey_static.m
@@ -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 .
-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;
diff --git a/matlab/dyn_second_order_solver.m b/matlab/dyn_second_order_solver.m
index 7e971bff8..9a4456428 100644
--- a/matlab/dyn_second_order_solver.m
+++ b/matlab/dyn_second_order_solver.m
@@ -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;
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index d1518c891..4932e3042 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -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
diff --git a/matlab/evaluate_static_model.m b/matlab/evaluate_static_model.m
index b05c65d41..e92931ca0 100644
--- a/matlab/evaluate_static_model.m
+++ b/matlab/evaluate_static_model.m
@@ -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 .
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
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 2bfa4a042..e5a6c7a47 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -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 .
-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
diff --git a/matlab/evaluate_steady_state_file.m b/matlab/evaluate_steady_state_file.m
index d13ac7091..14ef63a93 100644
--- a/matlab/evaluate_steady_state_file.m
+++ b/matlab/evaluate_steady_state_file.m
@@ -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 .
-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
diff --git a/matlab/expand_group.m b/matlab/expand_group.m
index 8e54f995a..ba953230a 100644
--- a/matlab/expand_group.m
+++ b/matlab/expand_group.m
@@ -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_});
diff --git a/matlab/get_identification_jacobians.m b/matlab/get_identification_jacobians.m
index d04d49266..c60c0d1af 100644
--- a/matlab/get_identification_jacobians.m
+++ b/matlab/get_identification_jacobians.m
@@ -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;
diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/get_perturbation_params_derivs.m
index 70540b1a3..fab20ef03 100644
--- a/matlab/get_perturbation_params_derivs.m
+++ b/matlab/get_perturbation_params_derivs.m
@@ -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 .
% =========================================================================
-% 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
diff --git a/matlab/get_perturbation_params_derivs_numerical_objective.m b/matlab/get_perturbation_params_derivs_numerical_objective.m
index a4540250e..efb997309 100644
--- a/matlab/get_perturbation_params_derivs_numerical_objective.m
+++ b/matlab/get_perturbation_params_derivs_numerical_objective.m
@@ -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
diff --git a/matlab/histvalf.m b/matlab/histvalf.m
index 112f57b8a..f4af96faf 100644
--- a/matlab/histvalf.m
+++ b/matlab/histvalf.m
@@ -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 .
-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
diff --git a/matlab/identification_numerical_objective.m b/matlab/identification_numerical_objective.m
index 4a0683fde..548e68784 100644
--- a/matlab/identification_numerical_objective.m
+++ b/matlab/identification_numerical_objective.m
@@ -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
diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m
index a9dc60ec9..28946e002 100644
--- a/matlab/k_order_pert.m
+++ b/matlab/k_order_pert.m
@@ -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.* ...
diff --git a/matlab/metropolis_run_analysis.m b/matlab/metropolis_run_analysis.m
index 4fd11ce99..0baa03ab2 100644
--- a/matlab/metropolis_run_analysis.m
+++ b/matlab/metropolis_run_analysis.m
@@ -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 .
-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);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
index f465588ed..13bb41205 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
@@ -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
diff --git a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
index 3f3a00648..333805e41 100644
--- a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
+++ b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
@@ -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 .
-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']);
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index f9dafab9a..9d0b45dcf 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -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 .
-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'),
diff --git a/matlab/perfect-foresight-models/sim1_lbj.m b/matlab/perfect-foresight-models/sim1_lbj.m
index 8ebb24e83..0ad115a2e 100644
--- a/matlab/perfect-foresight-models/sim1_lbj.m
+++ b/matlab/perfect-foresight-models/sim1_lbj.m
@@ -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 .
-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()
diff --git a/matlab/perfect-foresight-models/sim1_linear.m b/matlab/perfect-foresight-models/sim1_linear.m
index f6a67759b..0cf316034 100644
--- a/matlab/perfect-foresight-models/sim1_linear.m
+++ b/matlab/perfect-foresight-models/sim1_linear.m
@@ -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 .
-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
diff --git a/matlab/perfect-foresight-models/sim1_purely_backward.m b/matlab/perfect-foresight-models/sim1_purely_backward.m
index 904b6e647..2e2f50b32 100644
--- a/matlab/perfect-foresight-models/sim1_purely_backward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_backward.m
@@ -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 .
-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)
diff --git a/matlab/perfect-foresight-models/sim1_purely_forward.m b/matlab/perfect-foresight-models/sim1_purely_forward.m
index 58778ec57..b49d4180d 100644
--- a/matlab/perfect-foresight-models/sim1_purely_forward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_forward.m
@@ -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 .
-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
diff --git a/matlab/perfect-foresight-models/sim1_purely_static.m b/matlab/perfect-foresight-models/sim1_purely_static.m
index 76913d740..1500485da 100644
--- a/matlab/perfect-foresight-models/sim1_purely_static.m
+++ b/matlab/perfect-foresight-models/sim1_purely_static.m
@@ -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 .
-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
diff --git a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
index c937fdfa8..3e104717a 100644
--- a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
@@ -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 .
-[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
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
index f9c2c4fcb..2a5d0bbcd 100644
--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -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 .
-[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
diff --git a/matlab/pruned_state_space_system.m b/matlab/pruned_state_space_system.m
index 6d7159791..3f1f51a14 100644
--- a/matlab/pruned_state_space_system.m
+++ b/matlab/pruned_state_space_system.m
@@ -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
diff --git a/matlab/set_all_parameters.m b/matlab/set_all_parameters.m
index a69e95d53..c97f98aa2 100644
--- a/matlab/set_all_parameters.m
+++ b/matlab/set_all_parameters.m
@@ -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 .
-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
\ No newline at end of file
diff --git a/matlab/set_local_param_value.m b/matlab/set_local_param_value.m
index 28cbdfcc7..d2e608884 100644
--- a/matlab/set_local_param_value.m
+++ b/matlab/set_local_param_value.m
@@ -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 .
-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;
diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m
index 0d49b58f5..876876d8c 100644
--- a/matlab/solve_one_boundary.m
+++ b/matlab/solve_one_boundary.m
@@ -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);
diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m
index 416d6fbb7..b5daed82d 100644
--- a/matlab/solve_two_boundaries.m
+++ b/matlab/solve_two_boundaries.m
@@ -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 .
-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_resa1)
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
diff --git a/matlab/utilities/general/get_all_variables_but_lagged_leaded_exogenous.m b/matlab/utilities/general/get_all_variables_but_lagged_leaded_exogenous.m
index e34389ae0..3c3dee6a3 100644
--- a/matlab/utilities/general/get_all_variables_but_lagged_leaded_exogenous.m
+++ b/matlab/utilities/general/get_all_variables_but_lagged_leaded_exogenous.m
@@ -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 .
-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
\ No newline at end of file
diff --git a/matlab/utilities/general/indices_lagged_leaded_exogenous_variables.m b/matlab/utilities/general/indices_lagged_leaded_exogenous_variables.m
index 6683bd3b5..998e430c6 100644
--- a/matlab/utilities/general/indices_lagged_leaded_exogenous_variables.m
+++ b/matlab/utilities/general/indices_lagged_leaded_exogenous_variables.m
@@ -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);
diff --git a/matlab/utilities/general/log_variable.m b/matlab/utilities/general/log_variable.m
index c81325f09..45cdc74fb 100644
--- a/matlab/utilities/general/log_variable.m
+++ b/matlab/utilities/general/log_variable.m
@@ -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 .
-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)