From 0561200f1c9b67ff49366831444432d7594cc56b Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 25 Oct 2023 17:29:55 +0200 Subject: [PATCH] Implement even more naming consistency Mostly removes M for M_ --- matlab/AIM_first_order_solver.m | 12 +- matlab/check_for_calibrated_covariances.m | 78 ++--- matlab/dyn_ramsey_static.m | 70 ++--- matlab/dyn_second_order_solver.m | 42 +-- matlab/dynare_estimation_init.m | 34 +- matlab/evaluate_static_model.m | 18 +- matlab/evaluate_steady_state.m | 296 +++++++++--------- matlab/evaluate_steady_state_file.m | 32 +- matlab/expand_group.m | 46 +-- matlab/get_identification_jacobians.m | 36 +-- matlab/get_perturbation_params_derivs.m | 124 ++++---- ...bation_params_derivs_numerical_objective.m | 46 +-- matlab/histvalf.m | 28 +- matlab/identification_numerical_objective.m | 34 +- matlab/k_order_pert.m | 16 +- matlab/metropolis_run_analysis.m | 10 +- .../perfect_foresight_mcp_problem.m | 8 +- .../private/initialize_stacked_problem.m | 80 ++--- matlab/perfect-foresight-models/sim1.m | 70 ++--- matlab/perfect-foresight-models/sim1_lbj.m | 32 +- matlab/perfect-foresight-models/sim1_linear.m | 54 ++-- .../sim1_purely_backward.m | 44 +-- .../sim1_purely_forward.m | 46 +-- .../sim1_purely_static.m | 46 +-- .../solve_stacked_linear_problem.m | 38 +-- .../solve_stacked_problem.m | 52 +-- matlab/pruned_state_space_system.m | 62 ++-- matlab/set_all_parameters.m | 56 ++-- matlab/set_local_param_value.m | 9 +- matlab/solve_one_boundary.m | 60 ++-- matlab/solve_two_boundaries.m | 33 +- ...ll_variables_but_lagged_leaded_exogenous.m | 15 +- ...ndices_lagged_leaded_exogenous_variables.m | 11 +- matlab/utilities/general/log_variable.m | 8 +- 34 files changed, 826 insertions(+), 820 deletions(-) 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)