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

Further naming consistency improvements

kalman-mex
Johannes Pfeifer 2023-10-25 10:14:12 +02:00
parent 879d92fbf8
commit b3ce518433
71 changed files with 922 additions and 892 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

94
matlab/+bgp/write.m
View File

@ -1,5 +1,5 @@
function write(DynareModel) function write(M_)
% write(M_)
% Writes the nonlinear problem to be solved for computing the growth % Writes the nonlinear problem to be solved for computing the growth
% rates and levels along the Balanced Growth Path. Note that for % rates and levels along the Balanced Growth Path. Note that for
% the variables growing along the BGP, the identified levels are % the variables growing along the BGP, the identified levels are
@ -7,7 +7,7 @@ function write(DynareModel)
% sharing the same trend(s) are relevant. % sharing the same trend(s) are relevant.
% %
% INPUTS % INPUTS
% - DynareModel [struct] Dynare generated stucture describing the model (M_). % - M_ [struct] Dynare generated stucture describing the model (M_).
% %
% OUTPUTS % OUTPUTS
% None % None
@ -15,7 +15,7 @@ function write(DynareModel)
% REMARKS % REMARKS
% - The trends are assumed to be multiplicative. % - The trends are assumed to be multiplicative.
% Copyright © 2019 Dynare Team % Copyright © 2019-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -32,18 +32,18 @@ function write(DynareModel)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if DynareModel.maximum_lag && ~DynareModel.maximum_lead if M_.maximum_lag && ~M_.maximum_lead
i0 = transpose(DynareModel.lead_lag_incidence(1,:)); % Indices of the lagged variables. i0 = transpose(M_.lead_lag_incidence(1,:)); % Indices of the lagged variables.
i1 = transpose(DynareModel.lead_lag_incidence(2,:)); % Indices of the current variables. i1 = transpose(M_.lead_lag_incidence(2,:)); % Indices of the current variables.
i2 = []; % Indices of the leaded variables. i2 = []; % Indices of the leaded variables.
elseif DynareModel.maximum_lag && DynareModel.maximum_lead elseif M_.maximum_lag && M_.maximum_lead
i0 = transpose(DynareModel.lead_lag_incidence(1,:)); % Indices of the lagged variables. i0 = transpose(M_.lead_lag_incidence(1,:)); % Indices of the lagged variables.
i1 = transpose(DynareModel.lead_lag_incidence(2,:)); % Indices of the current variables. i1 = transpose(M_.lead_lag_incidence(2,:)); % Indices of the current variables.
i2 = transpose(DynareModel.lead_lag_incidence(3,:)); % Indices of the leaded variables. i2 = transpose(M_.lead_lag_incidence(3,:)); % Indices of the leaded variables.
elseif ~DynareModel.maximum_lag && DynareModel.maximum_lead elseif ~M_.maximum_lag && M_.maximum_lead
i0 = []; % Indices of the lagged variables. i0 = []; % Indices of the lagged variables.
i1 = transpose(DynareModel.lead_lag_incidence(1,:)); % Indices of the current variables. i1 = transpose(M_.lead_lag_incidence(1,:)); % Indices of the current variables.
i2 = transpose(DynareModel.lead_lag_incidence(2,:)); % Indices of the leaded variables. i2 = transpose(M_.lead_lag_incidence(2,:)); % Indices of the leaded variables.
else else
error('The model is static. The BGP is trivial.') error('The model is static. The BGP is trivial.')
end end
@ -71,7 +71,7 @@ else
end end
% Create function in mod namespace. % Create function in mod namespace.
fid = fopen(sprintf('+%s/bgpfun.m', DynareModel.fname), 'w'); fid = fopen(sprintf('+%s/bgpfun.m', M_.fname), 'w');
% Write header. % Write header.
fprintf(fid, 'function [F, JAC] = bgpfun(z)\n\n'); fprintf(fid, 'function [F, JAC] = bgpfun(z)\n\n');
@ -80,8 +80,8 @@ fprintf(fid, '%% This file has been generated by dynare (%s).\n\n', datestr(now)
% The function admits a unique vector as input argument. The first % The function admits a unique vector as input argument. The first
% half of the elements are for the levels of the endogenous % half of the elements are for the levels of the endogenous
% variables, the second half for the growth factors. % variables, the second half for the growth factors.
fprintf(fid, 'y = z(1:%u);\n\n', DynareModel.endo_nbr); fprintf(fid, 'y = z(1:%u);\n\n', M_.endo_nbr);
fprintf(fid, 'g = z(%u:%u);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr); fprintf(fid, 'g = z(%u:%u);\n', M_.endo_nbr+1, 2*M_.endo_nbr);
% Define the point where the dynamic model is to be evaluated. % Define the point where the dynamic model is to be evaluated.
fprintf(fid, 'Y = zeros(%u, 1);\n', 2*(n0+n1+n2)); fprintf(fid, 'Y = zeros(%u, 1);\n', 2*(n0+n1+n2));
@ -118,39 +118,39 @@ end
fprintf(fid, '\n'); fprintf(fid, '\n');
% Define the vector of parameters. % Define the vector of parameters.
fprintf(fid, 'p = zeros(%u, 1);\n', DynareModel.param_nbr); fprintf(fid, 'p = zeros(%u, 1);\n', M_.param_nbr);
for i = 1:DynareModel.param_nbr for i = 1:M_.param_nbr
fprintf(fid, 'p(%u) = %16.12f;\n', i, DynareModel.params(i)); fprintf(fid, 'p(%u) = %16.12f;\n', i, M_.params(i));
end end
fprintf(fid, '\n'); fprintf(fid, '\n');
% Initialize the vector holding the residuals over the two periods. % Initialize the vector holding the residuals over the two periods.
fprintf(fid, 'F = NaN(%u, 1);\n', 2*DynareModel.endo_nbr); fprintf(fid, 'F = NaN(%u, 1);\n', 2*M_.endo_nbr);
% Set vector of exogenous variables to 0. % Set vector of exogenous variables to 0.
fprintf(fid, 'x = zeros(1, %u);\n\n', DynareModel.exo_nbr); fprintf(fid, 'x = zeros(1, %u);\n\n', M_.exo_nbr);
% Evaluate the residuals and jacobian of the dynamic model in periods t and t+1. % Evaluate the residuals and jacobian of the dynamic model in periods t and t+1.
fprintf(fid, 'if nargout>1\n'); fprintf(fid, 'if nargout>1\n');
fprintf(fid, ' J = zeros(%u, %u);\n', 2*DynareModel.endo_nbr, n0+n1+n2+DynareModel.endo_nbr); fprintf(fid, ' J = zeros(%u, %u);\n', 2*M_.endo_nbr, n0+n1+n2+M_.endo_nbr);
fprintf(fid, ' [F(1:%u), tmp] = %s.dynamic(Y(1:%u), x, p, y, 1);\n', DynareModel.endo_nbr, DynareModel.fname, n0+n1+n2); fprintf(fid, ' [F(1:%u), tmp] = %s.dynamic(Y(1:%u), x, p, y, 1);\n', M_.endo_nbr, M_.fname, n0+n1+n2);
fprintf(fid, ' J(1:%u,1:%u) = tmp(:,1:%u);\n', DynareModel.endo_nbr, n0+n1+n2, n0+n1+n2); fprintf(fid, ' J(1:%u,1:%u) = tmp(:,1:%u);\n', M_.endo_nbr, n0+n1+n2, n0+n1+n2);
fprintf(fid, ' [F(%u:%u), tmp] = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr, DynareModel.fname, n0+n1+n2, 2*(n0+n1+n2)); fprintf(fid, ' [F(%u:%u), tmp] = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', M_.endo_nbr+1, 2*M_.endo_nbr, M_.fname, n0+n1+n2, 2*(n0+n1+n2));
fprintf(fid, ' J(%u:%u,1:%u) = tmp(:,1:%u);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr, n0+n1+n2, n0+n1+n2); fprintf(fid, ' J(%u:%u,1:%u) = tmp(:,1:%u);\n', M_.endo_nbr+1, 2*M_.endo_nbr, n0+n1+n2, n0+n1+n2);
fprintf(fid, 'else\n'); fprintf(fid, 'else\n');
fprintf(fid, ' F(1:%u) = %s.dynamic(Y(1:%u), x, p, y, 1);\n', DynareModel.endo_nbr, DynareModel.fname, n0+n1+n2); fprintf(fid, ' F(1:%u) = %s.dynamic(Y(1:%u), x, p, y, 1);\n', M_.endo_nbr, M_.fname, n0+n1+n2);
fprintf(fid, ' F(%u:%u) = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr, DynareModel.fname, n0+n1+n2, 2*(n0+n1+n2)); fprintf(fid, ' F(%u:%u) = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', M_.endo_nbr+1, 2*M_.endo_nbr, M_.fname, n0+n1+n2, 2*(n0+n1+n2));
fprintf(fid, 'end\n\n'); fprintf(fid, 'end\n\n');
% Compute the jacobian if required. % Compute the jacobian if required.
fprintf(fid, 'if nargout>1\n'); fprintf(fid, 'if nargout>1\n');
fprintf(fid, ' JAC = zeros(%u,%u);\n', 2*DynareModel.endo_nbr, 2*DynareModel.endo_nbr); fprintf(fid, ' JAC = zeros(%u,%u);\n', 2*M_.endo_nbr, 2*M_.endo_nbr);
% Compute the derivatives of the first block of equations (period t) % Compute the derivatives of the first block of equations (period t)
% with respect to the endogenous variables. % with respect to the endogenous variables.
if purely_backward_model || purely_forward_model if purely_backward_model || purely_forward_model
for i=1:DynareModel.eq_nbr for i=1:M_.eq_nbr
for j=1:DynareModel.endo_nbr for j=1:M_.endo_nbr
if I1(j) if I1(j)
if I0(j) if I0(j)
fprintf(fid, ' JAC(%u,%u) = J(%u,%u)+J(%u,%u)*g(%u);\n', i, j, i, I0(j), i, I1(j), j); fprintf(fid, ' JAC(%u,%u) = J(%u,%u)+J(%u,%u)*g(%u);\n', i, j, i, I0(j), i, I1(j), j);
@ -165,8 +165,8 @@ if purely_backward_model || purely_forward_model
end end
end end
else else
for i=1:DynareModel.eq_nbr for i=1:M_.eq_nbr
for j=1:DynareModel.endo_nbr for j=1:M_.endo_nbr
if I2(j) if I2(j)
if I1(j) if I1(j)
if I0(j) if I0(j)
@ -201,8 +201,8 @@ end
% Compute the derivatives of the second block of equations (period t+1) % Compute the derivatives of the second block of equations (period t+1)
% with respect to the endogenous variables. % with respect to the endogenous variables.
if purely_backward_model || purely_forward_model if purely_backward_model || purely_forward_model
for i=DynareModel.eq_nbr+1:2*DynareModel.eq_nbr for i=M_.eq_nbr+1:2*M_.eq_nbr
for j=1:DynareModel.endo_nbr for j=1:M_.endo_nbr
if I1(j) if I1(j)
if I0(j) if I0(j)
fprintf(fid, ' JAC(%u,%u) = J(%u,%u)*g(%u)+J(%u,%u)*g(%u)*g(%u);\n', i, j, i, I0(j), j, i, I1(j), j, j); fprintf(fid, ' JAC(%u,%u) = J(%u,%u)*g(%u)+J(%u,%u)*g(%u)*g(%u);\n', i, j, i, I0(j), j, i, I1(j), j, j);
@ -217,8 +217,8 @@ if purely_backward_model || purely_forward_model
end end
end end
else else
for i=DynareModel.eq_nbr+1:2*DynareModel.eq_nbr for i=M_.eq_nbr+1:2*M_.eq_nbr
for j=1:DynareModel.endo_nbr for j=1:M_.endo_nbr
if I2(j) if I2(j)
if I1(j) if I1(j)
if I0(j) if I0(j)
@ -253,16 +253,16 @@ end
% Compute the derivatives of the first block of equations (period t) % Compute the derivatives of the first block of equations (period t)
% with respect to the growth factors. % with respect to the growth factors.
if purely_backward_model || purely_forward_model if purely_backward_model || purely_forward_model
for i=1:DynareModel.eq_nbr for i=1:M_.eq_nbr
for j=1:DynareModel.endo_nbr for j=1:M_.endo_nbr
if I1(j) if I1(j)
fprintf(fid, ' J(%u,%u) = J(%u,%u)*y(%u);\n', i, n0+n1+n2+j, i, I1(j), j); fprintf(fid, ' J(%u,%u) = J(%u,%u)*y(%u);\n', i, n0+n1+n2+j, i, I1(j), j);
end end
end end
end end
else else
for i=1:DynareModel.eq_nbr for i=1:M_.eq_nbr
for j=1:DynareModel.endo_nbr for j=1:M_.endo_nbr
if I2(j) if I2(j)
if I1(j) if I1(j)
fprintf(fid, ' J(%u,%u) = J(%u,%u)*y(%u)+J(%u,%u)*2*g(%u)*y(%u);\n', i, n0+n1+n2+j, i, I1(j), j, i, I2(j), j, j); fprintf(fid, ' J(%u,%u) = J(%u,%u)*y(%u)+J(%u,%u)*2*g(%u)*y(%u);\n', i, n0+n1+n2+j, i, I1(j), j, i, I2(j), j, j);
@ -281,8 +281,8 @@ end
% Compute the derivatives of the second block of equations (period t+1) % Compute the derivatives of the second block of equations (period t+1)
% with respect to the endogenous variables. % with respect to the endogenous variables.
if purely_backward_model || purely_forward_model if purely_backward_model || purely_forward_model
for i=DynareModel.eq_nbr+1:2*DynareModel.eq_nbr for i=M_.eq_nbr+1:2*M_.eq_nbr
for j=1:DynareModel.endo_nbr for j=1:M_.endo_nbr
if I0(j) if I0(j)
fprintf(fid, ' J(%u,%u) = J(%u,%u)+J(%u,%u)*y(%u);\n', i, n0+n1+n2+j, i, n0+n1+n2+j, i, I0(j), j); fprintf(fid, ' J(%u,%u) = J(%u,%u)+J(%u,%u)*y(%u);\n', i, n0+n1+n2+j, i, n0+n1+n2+j, i, I0(j), j);
end end
@ -292,8 +292,8 @@ if purely_backward_model || purely_forward_model
end end
end end
else else
for i=DynareModel.eq_nbr+1:2*DynareModel.eq_nbr for i=M_.eq_nbr+1:2*M_.eq_nbr
for j=1:DynareModel.endo_nbr for j=1:M_.endo_nbr
if I2(j) if I2(j)
if I1(j) if I1(j)
if I0(j) if I0(j)
@ -325,7 +325,7 @@ else
end end
end end
fprintf(fid, ' JAC(:,%u:%u) = J(:,%u:%u);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr, n0+n1+n2+1, n0+n1+n2+DynareModel.endo_nbr); fprintf(fid, ' JAC(:,%u:%u) = J(:,%u:%u);\n', M_.endo_nbr+1, 2*M_.endo_nbr, n0+n1+n2+1, n0+n1+n2+M_.endo_nbr);
fprintf(fid,'end\n'); fprintf(fid,'end\n');
fclose(fid); fclose(fid);

2
matlab/+mom/run.m
View File

@ -554,7 +554,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
[xparam1, oo_, Woptflag] = mom.mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds); [xparam1, oo_, Woptflag] = mom.mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds);
% compute standard errors at mode % compute standard errors at mode
options_mom_.mom.vector_output = false; % make sure flag is reset options_mom_.mom.vector_output = false; % make sure flag is reset
M_ = set_all_parameters(xparam1,estim_params_,M_); % update M_ and DynareResults (in particular to get oo_.mom.model_moments) M_ = set_all_parameters(xparam1,estim_params_,M_); % update M_ and oo_ (in particular to get oo_.mom.model_moments)
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
options_mom_.mom.compute_derivs = true; % for GMM we compute derivatives analytically in the objective function with this flag options_mom_.mom.compute_derivs = true; % for GMM we compute derivatives analytically in the objective function with this flag
end end

62
matlab/+pac/+update/parameters.m
View File

@ -1,11 +1,11 @@
function DynareModel = parameters(pacname, DynareModel, DynareOutput, verbose) function M_ = parameters(pacname, M_, oo_, verbose)
% M_ = parameters(pacname, M_, oo_, verbose)
% Updates the parameters of a PAC equation. % Updates the parameters of a PAC equation.
% %
% INPUTS % INPUTS
% - pacname [string] Name of the pac equation. % - pacname [string] Name of the pac equation.
% - DynareModel [struct] M_ global structure (model properties) % - M_ [struct] M_ global structure (model properties)
% - DynareOutput [struct] oo_ global structure (model results) % - oo_ [struct] oo_ global structure (model results)
% %
% OUTPUTS % OUTPUTS
% - none % - none
@ -13,7 +13,7 @@ function DynareModel = parameters(pacname, DynareModel, DynareOutput, verbose)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2018-2021 Dynare Team % Copyright © 2018-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -40,27 +40,27 @@ if ~isrow(pacname)==1 || ~ischar(pacname)
end end
% Check the name of the PAC model. % Check the name of the PAC model.
if ~isfield(DynareModel.pac, pacname) if ~isfield(M_.pac, pacname)
error('PAC model %s is not defined in the model block!', pacname) error('PAC model %s is not defined in the model block!', pacname)
end end
% Get PAC model description % Get PAC model description
pacmodel = DynareModel.pac.(pacname); pacmodel = M_.pac.(pacname);
if pacmodel.model_consistent_expectations if pacmodel.model_consistent_expectations
error('This function cannot be used with Model Consistent Expectation. Try pac.mce.parameters instead.') error('This function cannot be used with Model Consistent Expectation. Try pac.mce.parameters instead.')
end end
% Get the name of the associated auxiliary model (VAR or TREND_COMPONENT) model and test its existence. % Get the name of the associated auxiliary model (VAR or TREND_COMPONENT) model and test its existence.
if ~isfield(DynareModel.(pacmodel.auxiliary_model_type), pacmodel.auxiliary_model_name) if ~isfield(M_.(pacmodel.auxiliary_model_type), pacmodel.auxiliary_model_name)
error('Unknown auxiliary model (%s) in PAC model (%s)!', pacmodel.auxiliary_model_name, pacname) error('Unknown auxiliary model (%s) in PAC model (%s)!', pacmodel.auxiliary_model_name, pacname)
end end
varmodel = DynareModel.(pacmodel.auxiliary_model_type).(pacmodel.auxiliary_model_name); varmodel = M_.(pacmodel.auxiliary_model_type).(pacmodel.auxiliary_model_name);
% Check that we have the values of the VAR or TREND_COMPONENT matrices. % Check that we have the values of the VAR or TREND_COMPONENT matrices.
if ~isfield(DynareOutput.(pacmodel.auxiliary_model_type), pacmodel.auxiliary_model_name) if ~isfield(oo_.(pacmodel.auxiliary_model_type), pacmodel.auxiliary_model_name)
error('Auxiliary model %s has to be estimated first!', pacmodel.auxiliary_model_name) error('Auxiliary model %s has to be estimated first!', pacmodel.auxiliary_model_name)
end end
varcalib = DynareOutput.(pacmodel.auxiliary_model_type).(pacmodel.auxiliary_model_name); varcalib = oo_.(pacmodel.auxiliary_model_type).(pacmodel.auxiliary_model_name);
if ~isfield(varcalib, 'CompanionMatrix') || any(isnan(varcalib.CompanionMatrix(:))) if ~isfield(varcalib, 'CompanionMatrix') || any(isnan(varcalib.CompanionMatrix(:)))
error('Auxiliary model %s has to be estimated first.', pacmodel.auxiliary_model_name) error('Auxiliary model %s has to be estimated first.', pacmodel.auxiliary_model_name)
end end
@ -79,13 +79,13 @@ else
end end
% Build the vector of PAC parameters (ECM parameter + autoregressive parameters). % Build the vector of PAC parameters (ECM parameter + autoregressive parameters).
pacvalues = DynareModel.params([pacmodel.ec.params; pacmodel.ar.params(1:pacmodel.max_lag)']); pacvalues = M_.params([pacmodel.ec.params; pacmodel.ar.params(1:pacmodel.max_lag)']);
% Get the indices for the stationary/nonstationary variables in the VAR system. % Get the indices for the stationary/nonstationary variables in the VAR system.
if numberofcomponents if numberofcomponents
id = cell(numberofcomponents, 1); id = cell(numberofcomponents, 1);
for i=1:numberofcomponents for i=1:numberofcomponents
id(i) = {find(strcmp(DynareModel.endo_names{pacmodel.components(i).endo_var}, varmodel.list_of_variables_in_companion_var))}; id(i) = {find(strcmp(M_.endo_names{pacmodel.components(i).endo_var}, varmodel.list_of_variables_in_companion_var))};
if isempty(id{i}) if isempty(id{i})
% Find the auxiliary variables if any % Find the auxiliary variables if any
ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false))); ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
@ -93,7 +93,7 @@ if numberofcomponents
error('Cannot find the trend variable in the Companion VAR/VECM model.') error('Cannot find the trend variable in the Companion VAR/VECM model.')
else else
for j=1:length(ad) for j=1:length(ad)
auxinfo = DynareModel.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(j)})); auxinfo = M_.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(j)}));
if isequal(auxinfo.endo_index, pacmodel.components(i).endo_var) if isequal(auxinfo.endo_index, pacmodel.components(i).endo_var)
id(i) = {ad(j)}; id(i) = {ad(j)};
break break
@ -110,7 +110,7 @@ if numberofcomponents
end end
end end
else else
id = {find(strcmp(DynareModel.endo_names{pacmodel.ec.vars(pacmodel.ec.istarget)}, varmodel.list_of_variables_in_companion_var))}; id = {find(strcmp(M_.endo_names{pacmodel.ec.vars(pacmodel.ec.istarget)}, varmodel.list_of_variables_in_companion_var))};
if isempty(id{1}) if isempty(id{1})
% Find the auxiliary variables if any % Find the auxiliary variables if any
ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false))); ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
@ -118,7 +118,7 @@ else
error('Cannot find the trend variable in the Companion VAR/VECM model.') error('Cannot find the trend variable in the Companion VAR/VECM model.')
else else
for i=1:length(ad) for i=1:length(ad)
auxinfo = DynareModel.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(i)})); auxinfo = M_.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(i)}));
if isequal(auxinfo.endo_index, pacmodel.ec.vars(pacmodel.ec.istarget)) if isequal(auxinfo.endo_index, pacmodel.ec.vars(pacmodel.ec.istarget))
id = {ad(i)}; id = {ad(i)};
break break
@ -180,7 +180,7 @@ else
end end
% Get the value of the discount factor. % Get the value of the discount factor.
beta = DynareModel.params(pacmodel.discount_index); beta = M_.params(pacmodel.discount_index);
% Is growth argument passed to pac_expectation? % Is growth argument passed to pac_expectation?
if isfield(pacmodel, 'growth_str') if isfield(pacmodel, 'growth_str')
@ -197,7 +197,7 @@ end
% Do we have rule of thumb agents? γ is the share of optimizing agents. % Do we have rule of thumb agents? γ is the share of optimizing agents.
if isfield(pacmodel, 'non_optimizing_behaviour') if isfield(pacmodel, 'non_optimizing_behaviour')
gamma = DynareModel.params(pacmodel.share_of_optimizing_agents_index); gamma = M_.params(pacmodel.share_of_optimizing_agents_index);
else else
gamma = 1.0; gamma = 1.0;
end end
@ -218,29 +218,29 @@ end
% Update M_.params with h % Update M_.params with h
if isequal(pacmodel.auxiliary_model_type, 'var') if isequal(pacmodel.auxiliary_model_type, 'var')
if DynareModel.var.(pacmodel.auxiliary_model_name).isconstant if M_.var.(pacmodel.auxiliary_model_name).isconstant
if isfield(pacmodel, 'h_param_indices') if isfield(pacmodel, 'h_param_indices')
% No decomposition % No decomposition
DynareModel.params(pacmodel.h_param_indices) = h{1}; M_.params(pacmodel.h_param_indices) = h{1};
else else
for i=1:numberofcomponents for i=1:numberofcomponents
DynareModel.params(pacmodel.components(i).h_param_indices) = h{i}; M_.params(pacmodel.components(i).h_param_indices) = h{i};
end end
end end
else else
if isfield(pacmodel, 'h_param_indices') if isfield(pacmodel, 'h_param_indices')
% No decomposition % No decomposition
DynareModel.params(pacmodel.h_param_indices(1)) = .0; M_.params(pacmodel.h_param_indices(1)) = .0;
DynareModel.params(pacmodel.h_param_indices(2:end)) = h{1}; M_.params(pacmodel.h_param_indices(2:end)) = h{1};
else else
for i=1:numberofcomponents for i=1:numberofcomponents
DynareModel.params(pacmodel.components(i).h_param_indices(1)) = .0; M_.params(pacmodel.components(i).h_param_indices(1)) = .0;
DynareModel.params(pacmodel.components(i).h_param_indices(2:end)) = h{i}; M_.params(pacmodel.components(i).h_param_indices(2:end)) = h{i};
end end
end end
end % If the auxiliary model (VAR) has no constant. end % If the auxiliary model (VAR) has no constant.
else else
DynareModel.params(pacmodel.h_param_indices) = h{1}; M_.params(pacmodel.h_param_indices) = h{1};
end % if auxiliary model is a VAR end % if auxiliary model is a VAR
% Update the parameter related to the growth neutrality correction. % Update the parameter related to the growth neutrality correction.
@ -263,7 +263,7 @@ if growth_flag
if isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i} if isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
tmp0 = tmp0 + pacmodel.optim_additive.scaling_factor(i); tmp0 = tmp0 + pacmodel.optim_additive.scaling_factor(i);
elseif ~isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i} elseif ~isnan(pacmodel.optim_additive.params(i)) && islogical(pacmodel.optim_additive.bgp{i}) && pacmodel.optim_additive.bgp{i}
tmp0 = tmp0 + DynareModel.params(pacmodel.optim_additive.params(i))*pacmodel.optim_additive.scaling_factor(i); tmp0 = tmp0 + M_.params(pacmodel.optim_additive.params(i))*pacmodel.optim_additive.scaling_factor(i);
elseif ~islogical(pacmodel.optim_additive.bgp{i}) elseif ~islogical(pacmodel.optim_additive.bgp{i})
error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.') error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
end end
@ -279,7 +279,7 @@ if growth_flag
if isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i} if isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i}
tmp0 = tmp0 + pacmodel.non_optimizing_behaviour.scaling_factor(i); tmp0 = tmp0 + pacmodel.non_optimizing_behaviour.scaling_factor(i);
elseif ~isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i} elseif ~isnan(pacmodel.non_optimizing_behaviour.params(i)) && islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && pacmodel.non_optimizing_behaviour.bgp{i}
tmp0 = tmp0 + DynareModel.params(pacmodel.non_optimizing_behaviour.params(i))*pacmodel.non_optimizing_behaviour.scaling_factor(i); tmp0 = tmp0 + M_.params(pacmodel.non_optimizing_behaviour.params(i))*pacmodel.non_optimizing_behaviour.scaling_factor(i);
elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && isnan(pacmodel.non_optimizing_behaviour.params(i)) elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && isnan(pacmodel.non_optimizing_behaviour.params(i))
tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.bgp{i}; tmp1 = tmp1 + pacmodel.non_optimizing_behaviour.scaling_factor(i)*pacmodel.non_optimizing_behaviour.bgp{i};
elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && ~isnan(pacmodel.non_optimizing_behaviour.params(i)) elseif ~islogical(pacmodel.non_optimizing_behaviour.bgp{i}) && isnumeric(pacmodel.non_optimizing_behaviour.bgp{i}) && ~isnan(pacmodel.non_optimizing_behaviour.params(i))
@ -298,7 +298,7 @@ if growth_flag
if isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i} if isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i}
tmp0 = tmp0 + pacmodel.additive.scaling_factor(i); tmp0 = tmp0 + pacmodel.additive.scaling_factor(i);
elseif ~isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i} elseif ~isnan(pacmodel.additive.params(i)) && islogical(pacmodel.additive.bgp{i}) && pacmodel.additive.bgp{i}
tmp0 = tmp0 + DynareModel.params(pacmodel.additive.params(i))*pacmodel.additive.scaling_factor(i); tmp0 = tmp0 + M_.params(pacmodel.additive.params(i))*pacmodel.additive.scaling_factor(i);
elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && isnan(pacmodel.additive.params(i)) elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && isnan(pacmodel.additive.params(i))
tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.bgp{i}; tmp1 = tmp1 + pacmodel.additive.scaling_factor(i)*pacmodel.additive.bgp{i};
elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && ~isnan(pacmodel.additive.params(i)) elseif ~islogical(pacmodel.additive.bgp{i}) && isnumeric(pacmodel.additive.bgp{i}) && ~isnan(pacmodel.additive.params(i))
@ -309,9 +309,9 @@ if growth_flag
ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation. ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
end end
if isfield(pacmodel, 'growth_neutrality_param_index') if isfield(pacmodel, 'growth_neutrality_param_index')
DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model. M_.params(pacmodel.growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
else else
DynareModel.params(pacmodel.components(j).growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model. M_.params(pacmodel.components(j).growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
end end
end end
end end

42
matlab/+var_expectation/update_parameters.m
View File

@ -1,19 +1,19 @@
function DynareModel = update_parameters(varexpectationmodelname, DynareModel, DynareOutput) function M_ = update_parameters(varexpectationmodelname, M_, oo_)
% function M_ = update_parameters(varexpectationmodelname, M_, oo_)
% Updates the VAR expectation reduced form parameters. % Updates the VAR expectation reduced form parameters.
% %
% INPUTS % INPUTS
% - varexpectationmodelname [string] Name of the pac equation. % - varexpectationmodelname [string] Name of the pac equation.
% - DynareModel [struct] M_ global structure (model properties) % - M_ [struct] global structure (model properties)
% - DynareOutput [struct] oo_ global structure (model results) % - oo_ [struct] oo_ global structure (model results)
% %
% OUTPUTS % OUTPUTS
% - DynareModel [struct] M_ global structure (with updated params field) % - M_ [struct] M_ global structure (with updated params field)
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2018-2021 Dynare Team % Copyright © 2018-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -36,26 +36,26 @@ if ~isrow(varexpectationmodelname)==1 || ~ischar(varexpectationmodelname)
end end
% Check that the model exists. % Check that the model exists.
if ~isfield(DynareModel.var_expectation, varexpectationmodelname) if ~isfield(M_.var_expectation, varexpectationmodelname)
error('VAR_EXPECTATION_MODEL %s is not defined!', varexpectationmodelname) error('VAR_EXPECTATION_MODEL %s is not defined!', varexpectationmodelname)
end end
% Get the VAR model description % Get the VAR model description
varexpectationmodel = DynareModel.var_expectation.(varexpectationmodelname); varexpectationmodel = M_.var_expectation.(varexpectationmodelname);
% Get the name of the associated VAR model and test its existence. % Get the name of the associated VAR model and test its existence.
if ~isfield(DynareModel.(varexpectationmodel.auxiliary_model_type), varexpectationmodel.auxiliary_model_name) if ~isfield(M_.(varexpectationmodel.auxiliary_model_type), varexpectationmodel.auxiliary_model_name)
error('Unknown VAR (%s) in VAR_EXPECTATION_MODEL (%s)!', varexpectationmodel.auxiliary_model_name, varexpectationmodelname) error('Unknown VAR (%s) in VAR_EXPECTATION_MODEL (%s)!', varexpectationmodel.auxiliary_model_name, varexpectationmodelname)
end end
auxmodel = DynareModel.(varexpectationmodel.auxiliary_model_type).(varexpectationmodel.auxiliary_model_name); auxmodel = M_.(varexpectationmodel.auxiliary_model_type).(varexpectationmodel.auxiliary_model_name);
% Check that we have the values of the VAR matrices. % Check that we have the values of the VAR matrices.
if ~isfield(DynareOutput.(varexpectationmodel.auxiliary_model_type), varexpectationmodel.auxiliary_model_name) if ~isfield(oo_.(varexpectationmodel.auxiliary_model_type), varexpectationmodel.auxiliary_model_name)
error('Auxiliary model %s has to be estimated or calibrated first!', varexpectationmodel.auxiliary_model_name) error('Auxiliary model %s has to be estimated or calibrated first!', varexpectationmodel.auxiliary_model_name)
end end
auxcalib = DynareOutput.(varexpectationmodel.auxiliary_model_type).(varexpectationmodel.auxiliary_model_name); auxcalib = oo_.(varexpectationmodel.auxiliary_model_type).(varexpectationmodel.auxiliary_model_name);
if ~isfield(auxcalib, 'CompanionMatrix') || any(isnan(auxcalib.CompanionMatrix(:))) if ~isfield(auxcalib, 'CompanionMatrix') || any(isnan(auxcalib.CompanionMatrix(:)))
message = sprintf('Auxiliary model %s has to be estimated first.', varexpectationmodel.auxiliary_model_name); message = sprintf('Auxiliary model %s has to be estimated first.', varexpectationmodel.auxiliary_model_name);
@ -68,7 +68,7 @@ if isfield(varexpectationmodel, 'discount_value')
discountfactor = varexpectationmodel.discount_value; discountfactor = varexpectationmodel.discount_value;
else else
if isfield(varexpectationmodel, 'discount_index') if isfield(varexpectationmodel, 'discount_index')
discountfactor = DynareModel.params(varexpectationmodel.discount_index); discountfactor = M_.params(varexpectationmodel.discount_index);
else else
error('This is most likely a bug. Pleasse conntact the Dynare Team.') error('This is most likely a bug. Pleasse conntact the Dynare Team.')
end end
@ -100,11 +100,11 @@ end
m = length(varexpectationmodel.expr.vars); m = length(varexpectationmodel.expr.vars);
variables_id_in_var = NaN(m,1); variables_id_in_var = NaN(m,1);
for i = 1:m for i = 1:m
j = find(strcmp(auxmodel.list_of_variables_in_companion_var, DynareModel.endo_names{varexpectationmodel.expr.vars(i)})); j = find(strcmp(auxmodel.list_of_variables_in_companion_var, M_.endo_names{varexpectationmodel.expr.vars(i)}));
if isempty(j) if isempty(j)
error('Cannot find variable %s in the companion VAR', DynareModel.endo_names{varexpectationmodel.expr.vars(i)}) error('Cannot find variable %s in the companion VAR', M_.endo_names{varexpectationmodel.expr.vars(i)})
else else
variables_id_in_var(i) = find(strcmp(auxmodel.list_of_variables_in_companion_var, DynareModel.endo_names{varexpectationmodel.expr.vars(i)})); variables_id_in_var(i) = find(strcmp(auxmodel.list_of_variables_in_companion_var, M_.endo_names{varexpectationmodel.expr.vars(i)}));
end end
end end
@ -210,12 +210,12 @@ end
% Update reduced form parameters in M_.params. % Update reduced form parameters in M_.params.
if isequal(varexpectationmodel.auxiliary_model_type, 'var') if isequal(varexpectationmodel.auxiliary_model_type, 'var')
if DynareModel.var.(varexpectationmodel.auxiliary_model_name).isconstant if M_.var.(varexpectationmodel.auxiliary_model_name).isconstant
DynareModel.params(varexpectationmodel.param_indices) = parameters; M_.params(varexpectationmodel.param_indices) = parameters;
else else
DynareModel.params(varexpectationmodel.param_indices(1)) = .0; M_.params(varexpectationmodel.param_indices(1)) = .0;
DynareModel.params(varexpectationmodel.param_indices(2:end)) = parameters; M_.params(varexpectationmodel.param_indices(2:end)) = parameters;
end end
else else
DynareModel.params(varexpectationmodel.param_indices) = parameters; M_.params(varexpectationmodel.param_indices) = parameters;
end end

64
matlab/@dprior/dprior.m
View File

@ -25,11 +25,11 @@ classdef dprior
methods methods
function o = dprior(BayesInfo, PriorTrunc, Uniform) function o = dprior(bayestopt_, PriorTrunc, Uniform)
% Class constructor. % Class constructor.
% %
% INPUTS % INPUTS
% - BayesInfo [struct] Informations about the prior distribution, aka bayestopt_. % - bayestopt_ [struct] Informations about the prior distribution, aka bayestopt_.
% - PriorTrunc [double] scalar, probability mass to be excluded, aka options_.prior_trunc % - PriorTrunc [double] scalar, probability mass to be excluded, aka options_.prior_trunc
% - Uniform [logical] scalar, produce uniform random deviates on the prior support. % - Uniform [logical] scalar, produce uniform random deviates on the prior support.
% %
@ -38,17 +38,17 @@ classdef dprior
% %
% REQUIREMENTS % REQUIREMENTS
% None. % None.
o.p6 = BayesInfo.p6; o.p6 = bayestopt_.p6;
o.p7 = BayesInfo.p7; o.p7 = bayestopt_.p7;
o.p3 = BayesInfo.p3; o.p3 = bayestopt_.p3;
o.p4 = BayesInfo.p4; o.p4 = bayestopt_.p4;
bounds = prior_bounds(BayesInfo, PriorTrunc); bounds = prior_bounds(bayestopt_, PriorTrunc);
o.lb = bounds.lb; o.lb = bounds.lb;
o.ub = bounds.ub; o.ub = bounds.ub;
if nargin>2 && Uniform if nargin>2 && Uniform
prior_shape = repmat(5, length(o.p6), 1); prior_shape = repmat(5, length(o.p6), 1);
else else
prior_shape = BayesInfo.pshape; prior_shape = bayestopt_.pshape;
end end
o.beta_index = find(prior_shape==1); o.beta_index = find(prior_shape==1);
if ~isempty(o.beta_index) if ~isempty(o.beta_index)
@ -246,23 +246,23 @@ end % classdef --*-- Unit tests --*--
%$ end %$ end
%$ end %$ end
%$ %$
%$ BayesInfo.pshape = p0; %$ bayestopt_.pshape = p0;
%$ BayesInfo.p1 = p1; %$ bayestopt_.p1 = p1;
%$ BayesInfo.p2 = p2; %$ bayestopt_.p2 = p2;
%$ BayesInfo.p3 = p3; %$ bayestopt_.p3 = p3;
%$ BayesInfo.p4 = p4; %$ bayestopt_.p4 = p4;
%$ BayesInfo.p5 = p5; %$ bayestopt_.p5 = p5;
%$ BayesInfo.p6 = p6; %$ bayestopt_.p6 = p6;
%$ BayesInfo.p7 = p7; %$ bayestopt_.p7 = p7;
%$ %$
%$ ndraws = 1e5; %$ ndraws = 1e5;
%$ m0 = BayesInfo.p1; %zeros(14,1); %$ m0 = bayestopt_.p1; %zeros(14,1);
%$ v0 = diag(BayesInfo.p2.^2); %zeros(14); %$ v0 = diag(bayestopt_.p2.^2); %zeros(14);
%$ %$
%$ % Call the tested routine %$ % Call the tested routine
%$ try %$ try
%$ % Instantiate dprior object %$ % Instantiate dprior object
%$ o = dprior(BayesInfo, prior_trunc, false); %$ o = dprior(bayestopt_, prior_trunc, false);
%$ % Do simulations in a loop and estimate recursively the mean and the variance. %$ % Do simulations in a loop and estimate recursively the mean and the variance.
%$ for i = 1:ndraws %$ for i = 1:ndraws
%$ draw = o.draw(); %$ draw = o.draw();
@ -277,8 +277,8 @@ end % classdef --*-- Unit tests --*--
%$ end %$ end
%$ %$
%$ if t(1) %$ if t(1)
%$ t(2) = all(abs(m0-BayesInfo.p1)<3e-3); %$ t(2) = all(abs(m0-bayestopt_.p1)<3e-3);
%$ t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3)); %$ t(3) = all(all(abs(v0-diag(bayestopt_.p2.^2))<5e-3));
%$ end %$ end
%$ T = all(t); %$ T = all(t);
%@eof:1 %@eof:1
@ -345,21 +345,21 @@ end % classdef --*-- Unit tests --*--
%$ end %$ end
%$ end %$ end
%$ %$
%$ BayesInfo.pshape = p0; %$ bayestopt_.pshape = p0;
%$ BayesInfo.p1 = p1; %$ bayestopt_.p1 = p1;
%$ BayesInfo.p2 = p2; %$ bayestopt_.p2 = p2;
%$ BayesInfo.p3 = p3; %$ bayestopt_.p3 = p3;
%$ BayesInfo.p4 = p4; %$ bayestopt_.p4 = p4;
%$ BayesInfo.p5 = p5; %$ bayestopt_.p5 = p5;
%$ BayesInfo.p6 = p6; %$ bayestopt_.p6 = p6;
%$ BayesInfo.p7 = p7; %$ bayestopt_.p7 = p7;
%$ %$
%$ ndraws = 1e5; %$ ndraws = 1e5;
%$ %$
%$ % Call the tested routine %$ % Call the tested routine
%$ try %$ try
%$ % Instantiate dprior object. %$ % Instantiate dprior object.
%$ o = dprior(BayesInfo, prior_trunc, false); %$ o = dprior(bayestopt_, prior_trunc, false);
%$ X = o.draws(ndraws); %$ X = o.draws(ndraws);
%$ m = mean(X, 2); %$ m = mean(X, 2);
%$ v = var(X, 0, 2); %$ v = var(X, 0, 2);
@ -369,8 +369,8 @@ end % classdef --*-- Unit tests --*--
%$ end %$ end
%$ %$
%$ if t(1) %$ if t(1)
%$ t(2) = all(abs(m-BayesInfo.p1)<3e-3); %$ t(2) = all(abs(m-bayestopt_.p1)<3e-3);
%$ t(3) = all(all(abs(v-BayesInfo.p2.^2)<5e-3)); %$ t(3) = all(all(abs(v-bayestopt_.p2.^2)<5e-3));
%$ end %$ end
%$ T = all(t); %$ T = all(t);
%@eof:2 %@eof:2

10
matlab/backward/backward_model_forecast.m
View File

@ -63,7 +63,7 @@ for i=1:M_.exo_nbr
end end
% Set up initial conditions % Set up initial conditions
[initialcondition, periods, innovations, DynareOptions, DynareModel, DynareOutput, endonames, exonames, dynamic_resid, dynamic_g1, y] = ... [initialcondition, periods, innovations, options_local, M_local, oo_local, endonames, exonames, dynamic_resid, dynamic_g1, y] = ...
simul_backward_model_init(initialcondition, periods, options_, M_, oo_, zeros(periods, M_.exo_nbr)); simul_backward_model_init(initialcondition, periods, options_, M_, oo_, zeros(periods, M_.exo_nbr));
% Get vector of indices for the selected endogenous variables. % Get vector of indices for the selected endogenous variables.
@ -90,9 +90,9 @@ end
% Compute forecast without shock % Compute forecast without shock
if options_.linear if options_.linear
[ysim__0, errorflag] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, dynamic_resid, dynamic_g1); [ysim__0, errorflag] = simul_backward_linear_model_(initialcondition, periods, options_local, M_local, oo_local, innovations, dynamic_resid, dynamic_g1);
else else
[ysim__0, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, dynamic_resid, dynamic_g1); [ysim__0, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, options_local, M_local, oo_local, innovations, dynamic_resid, dynamic_g1);
end end
if errorflag if errorflag
@ -110,9 +110,9 @@ if withuncertainty
for i=1:B for i=1:B
innovations = transpose(sigma*randn(M_.exo_nbr, periods)); innovations = transpose(sigma*randn(M_.exo_nbr, periods));
if options_.linear if options_.linear
[ysim__, xsim__, errorflag] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, dynamic_resid, dynamic_g1); [ysim__, xsim__, errorflag] = simul_backward_linear_model_(initialcondition, periods, options_local, M_local, oo_local, innovations, dynamic_resid, dynamic_g1);
else else
[ysim__, xsim__, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, dynamic_resid, dynamic_g1); [ysim__, xsim__, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, options_local, M_local, oo_local, innovations, dynamic_resid, dynamic_g1);
end end
if errorflag if errorflag
error('Simulation failed.') error('Simulation failed.')

16
matlab/backward/backward_model_irf.m
View File

@ -139,7 +139,7 @@ if ~isempty(innovationbaseline)
end end
% Set up initial conditions % Set up initial conditions
[initialcondition, periods, Innovations, DynareOptions, DynareModel, DynareOutput, endonames, exonames, dynamic_resid, dynamic_g1, y] = ... [initialcondition, periods, Innovations, options_local, M_local, oo_local, endonames, exonames, dynamic_resid, dynamic_g1, y] = ...
simul_backward_model_init(initialcondition, periods, options_, M_, oo_, Innovations); simul_backward_model_init(initialcondition, periods, options_, M_, oo_, Innovations);
% Get the covariance matrix of the shocks. % Get the covariance matrix of the shocks.
@ -161,9 +161,9 @@ irfs = struct();
% Baseline paths (get transition paths induced by the initial condition and % Baseline paths (get transition paths induced by the initial condition and
% baseline innovations). % baseline innovations).
if options_.linear if options_.linear
[ysim__0, errorflag] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, Innovations, dynamic_resid, dynamic_g1); [ysim__0, errorflag] = simul_backward_linear_model_(initialcondition, periods, options_local, M_local, oo_local, Innovations, dynamic_resid, dynamic_g1);
else else
[ysim__0, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, Innovations, dynamic_resid, dynamic_g1); [ysim__0, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, options_local, M_local, oo_local, Innovations, dynamic_resid, dynamic_g1);
end end
if errorflag if errorflag
@ -200,9 +200,9 @@ for i=1:length(listofshocks)
innovations(1,:) = innovations(1,:) + transpose(C(:,j)); innovations(1,:) = innovations(1,:) + transpose(C(:,j));
end end
if options_.linear if options_.linear
[ysim__1, errorflag] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, dynamic_resid, dynamic_g1); [ysim__1, errorflag] = simul_backward_linear_model_(initialcondition, periods, options_local, M_local, oo_local, innovations, dynamic_resid, dynamic_g1);
else else
[ysim__1, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, dynamic_resid, dynamic_g1); [ysim__1, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, options_local, M_local, oo_local, innovations, dynamic_resid, dynamic_g1);
end end
if errorflag if errorflag
warning('Simulation failed. Cannot compute IRF for %s.', listofshocks{i}) warning('Simulation failed. Cannot compute IRF for %s.', listofshocks{i})
@ -215,9 +215,9 @@ for i=1:length(listofshocks)
endo_simul__1 = feval(transform, ysim__1); endo_simul__1 = feval(transform, ysim__1);
end end
% Instantiate a dseries object (with all the endogenous variables) % Instantiate a dseries object (with all the endogenous variables)
alldeviations = dseries(transpose(endo_simul__1-endo_simul__0), initialcondition.init, endonames(1:M_.orig_endo_nbr), DynareModel.endo_names_tex(1:M_.orig_endo_nbr)); alldeviations = dseries(transpose(endo_simul__1-endo_simul__0), initialcondition.init, endonames(1:M_.orig_endo_nbr), M_local.endo_names_tex(1:M_.orig_endo_nbr));
if nargout>2 if nargout>2
allirfs = dseries(transpose(endo_simul__1), initialcondition.init, endonames(1:M_.orig_endo_nbr), DynareModel.endo_names_tex(1:M_.orig_endo_nbr)); allirfs = dseries(transpose(endo_simul__1), initialcondition.init, endonames(1:M_.orig_endo_nbr), M_local.endo_names_tex(1:M_.orig_endo_nbr));
end end
% Extract a sub-dseries object % Extract a sub-dseries object
if deterministicshockflag if deterministicshockflag
@ -234,6 +234,6 @@ for i=1:length(listofshocks)
end end
if nargout>1 if nargout>1
baseline = dseries(transpose(endo_simul__0), initialcondition.init, endonames(1:M_.orig_endo_nbr), DynareModel.endo_names_tex(1:M_.orig_endo_nbr)); baseline = dseries(transpose(endo_simul__0), initialcondition.init, endonames(1:M_.orig_endo_nbr), M_local.endo_names_tex(1:M_.orig_endo_nbr));
baseline = merge(baseline, innovationbaseline); baseline = merge(baseline, innovationbaseline);
end end

16
matlab/backward/calibrateresiduals.m
View File

@ -1,5 +1,5 @@
function [residuals, info] = calibrateresiduals(dbase, info, DynareModel) function [residuals, info] = calibrateresiduals(dbase, info, M_)
% [residuals, info] = calibrateresiduals(dbase, info, M_)
% Compute residuals in a backward model. Residuals are unobserved exogenous % Compute residuals in a backward model. Residuals are unobserved exogenous
% variables appearing additively in equations and without lags. An equation % variables appearing additively in equations and without lags. An equation
% cannot have more than one residual, and a residual cannot appear in more % cannot have more than one residual, and a residual cannot appear in more
@ -8,7 +8,7 @@ function [residuals, info] = calibrateresiduals(dbase, info, DynareModel)
% INPUTS % INPUTS
% - dbase [dseries] Object containing all the endogenous and observed exogenous variables. % - dbase [dseries] Object containing all the endogenous and observed exogenous variables.
% - info [struct] Informations about the residuals. % - info [struct] Informations about the residuals.
% - DynareModel [struct] M_ as produced by the preprocessor. % - M_ [struct] M_ as produced by the preprocessor.
% %
% OUTPUTS % OUTPUTS
% - residuals [dseries] Object containing the identified residuals. % - residuals [dseries] Object containing the identified residuals.
@ -38,13 +38,13 @@ function [residuals, info] = calibrateresiduals(dbase, info, DynareModel)
displayresidualsequationmapping = false; displayresidualsequationmapping = false;
% Get function handle for the dynamic model % Get function handle for the dynamic model
dynamic_resid = str2func([DynareModel.fname,'.sparse.dynamic_resid']); dynamic_resid = str2func([M_.fname,'.sparse.dynamic_resid']);
% Get data for all the endogenous variables. % Get data for all the endogenous variables.
ydata = dbase{info.endonames{:}}.data; ydata = dbase{info.endonames{:}}.data;
% Define function to retrieve an equation name % Define function to retrieve an equation name
eqname = @(z) DynareModel.equations_tags{cellfun(@(x) x==z, DynareModel.equations_tags(:,1)) & cellfun(@(x) isequal(x, 'name'), DynareModel.equations_tags(:,2)),3}; eqname = @(z) M_.equations_tags{cellfun(@(x) x==z, M_.equations_tags(:,1)) & cellfun(@(x) isequal(x, 'name'), M_.equations_tags(:,2)),3};
% Get data for all the exogenous variables. Missing exogenous variables, to be solved for, have NaN values. % Get data for all the exogenous variables. Missing exogenous variables, to be solved for, have NaN values.
exogenousvariablesindbase = intersect(info.exonames, dbase.name); exogenousvariablesindbase = intersect(info.exonames, dbase.name);
@ -56,7 +56,7 @@ xdata = allexogenousvariables.data;
% Evaluate the dynamic equation % Evaluate the dynamic equation
n = size(ydata, 2); n = size(ydata, 2);
y = [ydata(1,:)'; ydata(2,:)'; NaN(n, 1)]; y = [ydata(1,:)'; ydata(2,:)'; NaN(n, 1)];
r = dynamic_resid(y, xdata(2,:), DynareModel.params, zeros(n, 1)); r = dynamic_resid(y, xdata(2,:), M_.params, zeros(n, 1));
% Check that the number of equations evaluating to NaN matches the number of residuals % Check that the number of equations evaluating to NaN matches the number of residuals
idr = find(isnan(r)); idr = find(isnan(r));
@ -102,7 +102,7 @@ for i = 1:residuals.vobs
info.residualindex(i) = {strmatch(residualname, allexogenousvariables.name, 'exact')}; info.residualindex(i) = {strmatch(residualname, allexogenousvariables.name, 'exact')};
tmpxdata = xdata; tmpxdata = xdata;
tmpxdata(2, info.residualindex{i}) = 0; tmpxdata(2, info.residualindex{i}) = 0;
r = dynamic_resid(y, tmpxdata(2,:), DynareModel.params, zeros(n, 1)); r = dynamic_resid(y, tmpxdata(2,:), M_.params, zeros(n, 1));
info.equations(i) = { idr(find(~isnan(r(idr))))}; info.equations(i) = { idr(find(~isnan(r(idr))))};
end end
@ -130,7 +130,7 @@ xdata(:,cell2mat(info.residualindex)) = 0;
rdata = NaN(residuals.nobs, residuals.vobs); rdata = NaN(residuals.nobs, residuals.vobs);
for t=2:size(xdata, 1) for t=2:size(xdata, 1)
y = [ydata(t-1,:)'; ydata(t,:)'; NaN(n, 1)]; y = [ydata(t-1,:)'; ydata(t,:)'; NaN(n, 1)];
r = dynamic_resid(y, xdata(t,:), DynareModel.params, zeros(n, 1)); r = dynamic_resid(y, xdata(t,:), M_.params, zeros(n, 1));
rdata(t,:) = transpose(r(cell2mat(info.equations))); rdata(t,:) = transpose(r(cell2mat(info.equations)));
end end
residuals = dseries(rdata, dbase.init, info.residuals); residuals = dseries(rdata, dbase.init, info.residuals);

36
matlab/backward/checkdatabase.m
View File

@ -1,12 +1,12 @@
function [dbase, info] = checkdatabase(dbase, DynareModel, inversionflag, simulationflag) function [dbase, info] = checkdatabase(dbase, M_, inversionflag, simulationflag)
% [dbase, info] = checkdatabase(dbase, M_, inversionflag, simulationflag)
% Check that dbase contains all the endogenous variables of the model, and % Check that dbase contains all the endogenous variables of the model, and
% reorder the endogenous variables as declared in the mod file. If Dynare % reorder the endogenous variables as declared in the mod file. If Dynare
% adds auxiliary variables, for lags greater than 1 on endogenous variables, % adds auxiliary variables, for lags greater than 1 on endogenous variables,
% endogenous variables in difference (which may be lagged), or lags on the % endogenous variables in difference (which may be lagged), or lags on the
% exogenous variables, then thee routine complete the database. % exogenous variables, then thee routine complete the database.
% Copyright © 2018-2021 Dynare Team % Copyright © 2018-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -23,42 +23,42 @@ function [dbase, info] = checkdatabase(dbase, DynareModel, inversionflag, simula
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% if DynareModel.maximum_endo_lead % if M_.maximum_endo_lead
% error('The model (%s) is assumed to be backward!', DynareModel.fname) % error('The model (%s) is assumed to be backward!', M_.fname)
% end % end
if nargin<3 if nargin<3
inversionflag = false; inversionflag = false;
end end
if exist(sprintf('+%s/dynamic_set_auxiliary_series.m', DynareModel.fname), 'file') if exist(sprintf('+%s/dynamic_set_auxiliary_series.m', M_.fname), 'file')
dbase = feval(sprintf('%s.dynamic_set_auxiliary_series', DynareModel.fname), dbase, DynareModel.params); dbase = feval(sprintf('%s.dynamic_set_auxiliary_series', M_.fname), dbase, M_.params);
end end
listoflaggedexogenousvariables = {}; listoflaggedexogenousvariables = {};
if ~isempty(DynareModel.aux_vars) if ~isempty(M_.aux_vars)
listoflaggedexogenousvariables = DynareModel.exo_names([DynareModel.aux_vars(find([DynareModel.aux_vars.type]==3)).orig_index]); listoflaggedexogenousvariables = M_.exo_names([M_.aux_vars(find([M_.aux_vars.type]==3)).orig_index]);
end end
listoflaggedendogenousvariables = {}; listoflaggedendogenousvariables = {};
laggedendogenousvariablesidx = find(DynareModel.lead_lag_incidence(1,1:DynareModel.orig_endo_nbr)); laggedendogenousvariablesidx = find(M_.lead_lag_incidence(1,1:M_.orig_endo_nbr));
if ~isempty(laggedendogenousvariablesidx) if ~isempty(laggedendogenousvariablesidx)
listoflaggedendogenousvariables = DynareModel.endo_names(laggedendogenousvariablesidx); listoflaggedendogenousvariables = M_.endo_names(laggedendogenousvariablesidx);
end end
if ~isempty(DynareModel.aux_vars) if ~isempty(M_.aux_vars)
laggedendogenousvariablesidx = find([DynareModel.aux_vars.type]==1); laggedendogenousvariablesidx = find([M_.aux_vars.type]==1);
if ~isempty(laggedendogenousvariablesidx) if ~isempty(laggedendogenousvariablesidx)
listoflaggedendogenousvariables = union(listoflaggedendogenousvariables, DynareModel.endo_names([DynareModel.aux_vars(laggedendogenousvariablesidx).orig_index])); listoflaggedendogenousvariables = union(listoflaggedendogenousvariables, M_.endo_names([M_.aux_vars(laggedendogenousvariablesidx).orig_index]));
end end
laggedendogenousvariablesidx = find([DynareModel.aux_vars.type]==8); laggedendogenousvariablesidx = find([M_.aux_vars.type]==8);
if ~isempty(laggedendogenousvariablesidx) if ~isempty(laggedendogenousvariablesidx)
listoflaggedendogenousvariables = union(listoflaggedendogenousvariables, DynareModel.endo_names([DynareModel.aux_vars(laggedendogenousvariablesidx).orig_index])); listoflaggedendogenousvariables = union(listoflaggedendogenousvariables, M_.endo_names([M_.aux_vars(laggedendogenousvariablesidx).orig_index]));
end end
end end
info = struct; info = struct;
info.endonames = DynareModel.endo_names; info.endonames = M_.endo_names;
info.exonames = DynareModel.exo_names; info.exonames = M_.exo_names;
info.computeresiduals = false; info.computeresiduals = false;
% Check that all the endogenous variables are defined in dbase. % Check that all the endogenous variables are defined in dbase.

8
matlab/backward/checkdatabaseforinversion.m
View File

@ -1,11 +1,11 @@
function [dbase, info] = checkdatabaseforinversion(dbase, DynareModel) function [dbase, info] = checkdatabaseforinversion(dbase, M_)
% [dbase, info] = checkdatabaseforinversion(dbase, M_)
% Check that dbase contains all the endogenous variables of the model, and % Check that dbase contains all the endogenous variables of the model, and
% reorder the endogenous variables as declared in the mod file. If Dynare % reorder the endogenous variables as declared in the mod file. If Dynare
% adds auxiliary variables, for lags greater than 1 on endogebnous variables % adds auxiliary variables, for lags greater than 1 on endogebnous variables
% or lags on the exogenous variables. % or lags on the exogenous variables.
% Copyright © 2017-2018 Dynare Team % Copyright © 2017-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -22,4 +22,4 @@ function [dbase, info] = checkdatabaseforinversion(dbase, DynareModel)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[dbase, info] = checkdatabase(dbase, DynareModel, true, false); [dbase, info] = checkdatabase(dbase, M_, true, false);

18
matlab/backward/get_lags_on_endogenous_variables.m
View File

@ -1,8 +1,8 @@
function l = get_lags_on_endogenous_variables(DynareModel) function l = get_lags_on_endogenous_variables(M_)
% l = get_lags_on_endogenous_variables(M_)
% Returns a vector with the max lag for each endogenous variable. % Returns a vector with the max lag for each endogenous variable.
% Copyright © 2017 Dynare Team % Copyright © 2017-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -19,13 +19,13 @@ function l = get_lags_on_endogenous_variables(DynareModel)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
l = zeros(DynareModel.orig_endo_nbr, 1); l = zeros(M_.orig_endo_nbr, 1);
l(find(DynareModel.lead_lag_incidence(1,1:DynareModel.orig_endo_nbr))) = -1; l(find(M_.lead_lag_incidence(1,1:M_.orig_endo_nbr))) = -1;
if ~isempty(DynareModel.aux_vars) if ~isempty(M_.aux_vars)
aux_var_for_lagged_endogenous = find([DynareModel.aux_vars(:).type]==1); aux_var_for_lagged_endogenous = find([M_.aux_vars(:).type]==1);
for i=1:length(aux_var_for_lagged_endogenous) for i=1:length(aux_var_for_lagged_endogenous)
l(DynareModel.aux_vars(aux_var_for_lagged_endogenous(i)).orig_index) = ... l(M_.aux_vars(aux_var_for_lagged_endogenous(i)).orig_index) = ...
DynareModel.aux_vars(aux_var_for_lagged_endogenous(i)).orig_lead_lag; M_.aux_vars(aux_var_for_lagged_endogenous(i)).orig_lead_lag;
end end
end end

16
matlab/backward/get_lags_on_exogenous_variables.m
View File

@ -1,8 +1,8 @@
function l = get_lags_on_exogenous_variables(DynareModel) function l = get_lags_on_exogenous_variables(M_)
% l = get_lags_on_exogenous_variables(M_)
% Returns a vector with the max lag for each exogenous variable. % Returns a vector with the max lag for each exogenous variable.
% Copyright © 2017 Dynare Team % Copyright © 2017-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -19,12 +19,12 @@ function l = get_lags_on_exogenous_variables(DynareModel)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
l = zeros(DynareModel.exo_nbr, 1); l = zeros(M_.exo_nbr, 1);
if ~isempty(DynareModel.aux_vars) if ~isempty(M_.aux_vars)
aux_var_for_lagged_exogenous = find([DynareModel.aux_vars(:).type]==3); aux_var_for_lagged_exogenous = find([M_.aux_vars(:).type]==3);
for i=1:length(aux_var_for_lagged_exogenous) for i=1:length(aux_var_for_lagged_exogenous)
l(DynareModel.aux_vars(aux_var_for_lagged_exogenous(i)).orig_index) = ... l(M_.aux_vars(aux_var_for_lagged_exogenous(i)).orig_index) = ...
DynareModel.aux_vars(aux_var_for_lagged_exogenous(i)).orig_lead_lag-1; M_.aux_vars(aux_var_for_lagged_exogenous(i)).orig_lead_lag-1;
end end
end end

4
matlab/backward/simul_backward_model.m
View File

@ -12,8 +12,8 @@ function [simulation, errorflag] = simul_backward_model(initialconditions, sampl
% - errorflag [logical] scalar, equal to false iff the simulation did not fail. % - errorflag [logical] scalar, equal to false iff the simulation did not fail.
% %
% REMARKS % REMARKS
% [1] The innovations used for the simulation are saved in DynareOutput.exo_simul, and the resulting paths for the endogenous % [1] The innovations used for the simulation are saved in oo_.exo_simul, and the resulting paths for the endogenous
% variables are saved in DynareOutput.endo_simul. % variables are saved in oo_.endo_simul.
% [2] The last input argument is not mandatory. If absent we use random draws and rescale them with the informations provided % [2] The last input argument is not mandatory. If absent we use random draws and rescale them with the informations provided
% through the shocks block. % through the shocks block.
% [3] If the first input argument is empty, the endogenous variables are initialized with 0, or if available with the informations % [3] If the first input argument is empty, the endogenous variables are initialized with 0, or if available with the informations

14
matlab/check_dsge_var_model.m
View File

@ -1,4 +1,4 @@
function check_dsge_var_model(Model, EstimatedParameters, BayesInfo) function check_dsge_var_model(M_, estim_params_, bayestopt_)
% Check if the dsge model can be estimated with the DSGE-VAR approach. % Check if the dsge model can be estimated with the DSGE-VAR approach.
@ -19,26 +19,26 @@ function check_dsge_var_model(Model, EstimatedParameters, BayesInfo)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if EstimatedParameters.nvn if estim_params_.nvn
error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!') error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!')
end end
if EstimatedParameters.ncn if estim_params_.ncn
error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!') error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!')
end end
if any(vec(Model.H)) if any(vec(M_.H))
error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!') error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!')
end end
if EstimatedParameters.ncx if estim_params_.ncx
error('Estimation::DsgeVarLikelihood: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.') error('Estimation::DsgeVarLikelihood: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.')
end end
if Model.exo_nbr>1 && any(vec(tril(Model.Sigma_e,-1))) if M_.exo_nbr>1 && any(vec(tril(M_.Sigma_e,-1)))
error('Estimation::DsgeVarLikelihood: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.') error('Estimation::DsgeVarLikelihood: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.')
end end
if isequal(BayesInfo.with_trend,1) if isequal(bayestopt_.with_trend,1)
error('Estimation::DsgeVarLikelihood: Linear trend is not yet implemented!') error('Estimation::DsgeVarLikelihood: Linear trend is not yet implemented!')
end end

38
matlab/cli/prior.m
View File

@ -61,7 +61,7 @@ if (size(estim_params_.var_endo,1) || size(estim_params_.corrn,1))
end end
% Fill or update bayestopt_ structure % Fill or update bayestopt_ structure
[xparam1, estim_params_, BayesOptions, lb, ub, Model] = set_prior(estim_params_, M_, options_); [xparam1, estim_params_, BayesOptions, lb, ub, M_local] = set_prior(estim_params_, M_, options_);
% Set restricted state space % Set restricted state space
options_plot_priors_old=options_.plot_priors; options_plot_priors_old=options_.plot_priors;
options_.plot_priors=0; options_.plot_priors=0;
@ -77,27 +77,27 @@ if isempty(options_.qz_criterium)
changed_qz_criterium_flag = 1; changed_qz_criterium_flag = 1;
end end
Model.dname = Model.fname; M_local.dname = M_local.fname;
% Temporarly set options_.order equal to one % Temporarly set options_.order equal to one
order = options_.order; order = options_.order;
options_.order = 1; options_.order = 1;
if ismember('plot', varargin) if ismember('plot', varargin)
plot_priors(BayesOptions, Model, estim_params_, options_) plot_priors(BayesOptions, M_local, estim_params_, options_)
donesomething = true; donesomething = true;
end end
if ismember('table', varargin) if ismember('table', varargin)
print_table_prior(lb, ub, options_, Model, BayesOptions, estim_params_); print_table_prior(lb, ub, options_, M_local, BayesOptions, estim_params_);
donesomething = true; donesomething = true;
end end
if ismember('simulate', varargin) % Prior simulations (BK). if ismember('simulate', varargin) % Prior simulations (BK).
if ismember('moments(distribution)', varargin) if ismember('moments(distribution)', varargin)
results = prior_sampler(1, Model, BayesOptions, options_, oo_, estim_params_); results = prior_sampler(1, M_local, BayesOptions, options_, oo_, estim_params_);
else else
results = prior_sampler(0, Model, BayesOptions, options_, oo_, estim_params_); results = prior_sampler(0, M_local, BayesOptions, options_, oo_, estim_params_);
end end
% Display prior mass info % Display prior mass info
skipline(2) skipline(2)
@ -119,7 +119,7 @@ if ismember('simulate', varargin) % Prior simulations (BK).
end end
if ismember('optimize', varargin) % Prior optimization. if ismember('optimize', varargin) % Prior optimization.
optimize_prior(options_, Model, oo_, BayesOptions, estim_params_); optimize_prior(options_, M_local, oo_, BayesOptions, estim_params_);
donesomething = true; donesomething = true;
end end
@ -130,16 +130,16 @@ if ismember('moments', varargin) % Prior simulations (2nd order moments).
k = find(isnan(xparam1)); k = find(isnan(xparam1));
xparam1(k) = BayesOptions.p1(k); xparam1(k) = BayesOptions.p1(k);
% Update vector of parameters and covariance matrices % Update vector of parameters and covariance matrices
Model = set_all_parameters(xparam1, estim_params_, Model); M_local = set_all_parameters(xparam1, estim_params_, M_local);
% Check model. % Check model.
check_model(Model); check_model(M_local);
% Compute state space representation of the model. % Compute state space representation of the model.
oo__ = oo_; oo__ = oo_;
oo__.dr = set_state_space(oo__.dr, Model); oo__.dr = set_state_space(oo__.dr, M_local);
% Solve model % Solve model
[T,R,~,info,oo__.dr, Model.params] = dynare_resolve(Model , options_ , oo__.dr, oo__.steady_state, oo__.exo_steady_state, oo__.exo_det_steady_state,'restrict'); [T,R,~,info,oo__.dr, M_local.params] = dynare_resolve(M_local , options_ , oo__.dr, oo__.steady_state, oo__.exo_steady_state, oo__.exo_det_steady_state,'restrict');
if ~info(1) if ~info(1)
info=endogenous_prior_restrictions(T,R,Model , options__ , oo__.dr,oo__.steady_state,oo__.exo_steady_state,oo__.exo_det_steady_state); info=endogenous_prior_restrictions(T,R,M_local , options__ , oo__.dr,oo__.steady_state,oo__.exo_steady_state,oo__.exo_det_steady_state);
end end
if info if info
skipline() skipline()
@ -149,19 +149,19 @@ if ismember('moments', varargin) % Prior simulations (2nd order moments).
return return
end end
% Compute and display second order moments % Compute and display second order moments
oo__ = disp_th_moments(oo__.dr, [], Model, options__, oo__); oo__ = disp_th_moments(oo__.dr, [], M_local, options__, oo__);
skipline(2) skipline(2)
donesomething = true; donesomething = true;
end end
if ismember('moments(distribution)', varargin) % Prior simulations (BK). if ismember('moments(distribution)', varargin) % Prior simulations (BK).
if ~ismember('simulate', varargin) if ~ismember('simulate', varargin)
results = prior_sampler(1, Model, BayesOptions, options_, oo_, estim_params_); results = prior_sampler(1, M_local, BayesOptions, options_, oo_, estim_params_);
end end
priorpath = [Model.dname filesep() 'prior' filesep() 'draws' filesep()]; priorpath = [M_local.dname filesep() 'prior' filesep() 'draws' filesep()];
list_of_files = dir([priorpath 'prior_draws*']); list_of_files = dir([priorpath 'prior_draws*']);
FirstOrderMoments = NaN(Model.orig_endo_nbr, options_.prior_mc); FirstOrderMoments = NaN(M_local.orig_endo_nbr, options_.prior_mc);
SecondOrderMoments = NaN(Model.orig_endo_nbr, Model.orig_endo_nbr, options_.prior_mc); SecondOrderMoments = NaN(M_local.orig_endo_nbr, M_local.orig_endo_nbr, options_.prior_mc);
iter = 1; iter = 1;
noprint = options_.noprint; noprint = options_.noprint;
options_.noprint = 1; options_.noprint = 1;
@ -172,8 +172,8 @@ if ismember('moments(distribution)', varargin) % Prior simulations (BK).
dr = tmp.pdraws{j,3}; dr = tmp.pdraws{j,3};
oo__ = oo_; oo__ = oo_;
oo__.dr = dr; oo__.dr = dr;
Model=set_parameters_locally(Model,tmp.pdraws{j,1});% Needed to update the covariance matrix of the state innovations. M_local=set_parameters_locally(M_local,tmp.pdraws{j,1});% Needed to update the covariance matrix of the state innovations.
oo__ = disp_th_moments(oo__.dr, [], Model, options_, oo__); oo__ = disp_th_moments(oo__.dr, [], M_local, options_, oo__);
FirstOrderMoments(:,iter) = oo__.mean; FirstOrderMoments(:,iter) = oo__.mean;
SecondOrderMoments(:,:,iter) = oo__.var; SecondOrderMoments(:,:,iter) = oo__.var;
iter = iter+1; iter = iter+1;

20
matlab/disp_steady_state.m
View File

@ -1,11 +1,11 @@
function disp_steady_state(M,oo,options) function disp_steady_state(M_,oo_,options_)
% function disp_steady_state(M,oo,options) % function disp_steady_state(M_,oo_,options_)
% computes and prints the steady state calculations % computes and prints the steady state calculations
% %
% INPUTS % INPUTS
% M structure of parameters % M_ structure of parameters
% oo structure of results % oo_ structure of results
% options structure of options % options_ structure of options
% %
% OUTPUTS % OUTPUTS
% none % none
@ -13,7 +13,7 @@ function disp_steady_state(M,oo,options)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2001-2020 Dynare Team % Copyright © 2001-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -31,15 +31,15 @@ function disp_steady_state(M,oo,options)
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
skipline() skipline()
if options.loglinear if options_.loglinear
disp('STEADY-STATE RESULTS FOR THE UNLOGGED VARIABLES:') disp('STEADY-STATE RESULTS FOR THE UNLOGGED VARIABLES:')
else else
disp('STEADY-STATE RESULTS:') disp('STEADY-STATE RESULTS:')
end end
skipline() skipline()
endo_names = char(M.endo_names); endo_names = char(M_.endo_names);
steady_state = oo.steady_state; steady_state = oo_.steady_state;
for i = 1:M.orig_endo_nbr for i = 1:M_.orig_endo_nbr
fprintf('%s \t\t %g\n', endo_names(i,:), steady_state(i)); fprintf('%s \t\t %g\n', endo_names(i,:), steady_state(i));
end end

150
matlab/dyn_forecast.m
View File

@ -1,15 +1,15 @@
function forecast = dyn_forecast(var_list,M,options,oo,task,dataset_info) function forecast = dyn_forecast(var_list,M_,options_,oo_,task,dataset_info)
% function forecast = dyn_forecast(var_list,M,options,oo,task,dataset_info) % function forecast = dyn_forecast(var_list,M_,options_,oo_,task,dataset_info)
% computes mean forecast for a given value of the parameters % computes mean forecast for a given value of the parameters
% computes also confidence bands for the forecast % computes also confidence bands for the forecast
% %
% INPUTS % INPUTS
% var_list: list of variables (character matrix) % var_list: list of variables (character matrix)
% M: Dynare model structure % M_: Dynare model structure
% options: Dynare options structure % options_: Dynare options structure
% oo: Dynare results structure % oo_: Dynare results structure
% task: indicates how to initialize the forecast % task: indicates how to initialize the forecast
% either 'simul' or 'smoother' % either 'simul' or 'smoother'
% dataset_info: Various informations about the dataset (descriptive statistics and missing observations). % dataset_info: Various informations about the dataset (descriptive statistics and missing observations).
% OUTPUTS % OUTPUTS
@ -27,7 +27,7 @@ function forecast = dyn_forecast(var_list,M,options,oo,task,dataset_info)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2003-2022 Dynare Team % Copyright © 2003-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -44,23 +44,23 @@ function forecast = dyn_forecast(var_list,M,options,oo,task,dataset_info)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if ~isfield(oo,'dr') || isempty(oo.dr) if ~isfield(oo_,'dr') || isempty(oo_.dr)
error('dyn_forecast: the decision rules have not been computed. Did you forget a stoch_simul-command?') error('dyn_forecast: the decision rules have not been computed. Did you forget a stoch_simul-command?')
end end
if nargin<6 && options.prefilter if nargin<6 && options_.prefilter
error('The prefiltering option is not allowed without providing a dataset') error('The prefiltering option is not allowed without providing a dataset')
elseif nargin==6 elseif nargin==6
mean_varobs=dataset_info.descriptive.mean'; mean_varobs=dataset_info.descriptive.mean';
end end
oo=make_ex_(M,options,oo); oo_=make_ex_(M_,options_,oo_);
maximum_lag = M.maximum_lag; maximum_lag = M_.maximum_lag;
endo_names = M.endo_names; endo_names = M_.endo_names;
if isempty(var_list) if isempty(var_list)
var_list = endo_names(1:M.orig_endo_nbr); var_list = endo_names(1:M_.orig_endo_nbr);
end end
i_var = []; i_var = [];
for i = 1:length(var_list) for i = 1:length(var_list)
@ -76,108 +76,108 @@ n_var = length(i_var);
trend = 0; trend = 0;
switch task switch task
case 'simul' case 'simul'
horizon = options.periods; horizon = options_.periods;
if horizon == 0 if horizon == 0
horizon = 5; horizon = 5;
end end
if isempty(M.endo_histval) if isempty(M_.endo_histval)
if options.loglinear && ~options.logged_steady_state if options_.loglinear && ~options_.logged_steady_state
y0 = repmat(log(oo.dr.ys),1,maximum_lag); y0 = repmat(log(oo_.dr.ys),1,maximum_lag);
else else
y0 = repmat(oo.dr.ys,1,maximum_lag); y0 = repmat(oo_.dr.ys,1,maximum_lag);
end end
else else
if options.loglinear if options_.loglinear
y0 = log_variable(1:M.endo_nbr,M.endo_histval,M); y0 = log_variable(1:M_.endo_nbr,M_.endo_histval,M_);
else else
y0 = M.endo_histval; y0 = M_.endo_histval;
end end
end end
case 'smoother' case 'smoother'
horizon = options.forecast; horizon = options_.forecast;
if isnan(options.first_obs) if isnan(options_.first_obs)
first_obs=1; first_obs=1;
else else
first_obs=options.first_obs; first_obs=options_.first_obs;
end end
if isfield(oo.SmoothedVariables,'Mean') if isfield(oo_.SmoothedVariables,'Mean')
y_smoothed = oo.SmoothedVariables.Mean; y_smoothed = oo_.SmoothedVariables.Mean;
else else
y_smoothed = oo.SmoothedVariables; y_smoothed = oo_.SmoothedVariables;
end end
y0 = zeros(M.endo_nbr,maximum_lag); y0 = zeros(M_.endo_nbr,maximum_lag);
for i = 1:M.endo_nbr for i = 1:M_.endo_nbr
v_name = M.endo_names{i}; v_name = M_.endo_names{i};
y0(i,:) = y_smoothed.(v_name)(end-maximum_lag+1:end); %includes steady state or mean, but simult_ will subtract only steady state y0(i,:) = y_smoothed.(v_name)(end-maximum_lag+1:end); %includes steady state or mean, but simult_ will subtract only steady state
% 2. Subtract mean/steady state and add steady state; takes care of prefiltering % 2. Subtract mean/steady state and add steady state; takes care of prefiltering
if isfield(oo.Smoother,'Constant') && isfield(oo.Smoother.Constant,v_name) if isfield(oo_.Smoother,'Constant') && isfield(oo_.Smoother.Constant,v_name)
y0(i,:)=y0(i,:)-oo.Smoother.Constant.(v_name)(end-maximum_lag+1:end); %subtract mean or steady state y0(i,:)=y0(i,:)-oo_.Smoother.Constant.(v_name)(end-maximum_lag+1:end); %subtract mean or steady state
if options.loglinear if options_.loglinear
y0(i,:)=y0(i,:)+log_variable(i,oo.dr.ys,M); y0(i,:)=y0(i,:)+log_variable(i,oo_.dr.ys,M_);
else else
y0(i,:)=y0(i,:)+oo.dr.ys(strmatch(v_name, M.endo_names, 'exact')); y0(i,:)=y0(i,:)+oo_.dr.ys(strmatch(v_name, M_.endo_names, 'exact'));
end end
end end
% 2. Subtract trend % 2. Subtract trend
if isfield(oo.Smoother,'Trend') && isfield(oo.Smoother.Trend,v_name) if isfield(oo_.Smoother,'Trend') && isfield(oo_.Smoother.Trend,v_name)
y0(i,:)=y0(i,:)-oo.Smoother.Trend.(v_name)(end-maximum_lag+1:end); %subtract trend, which is not subtracted by simult_ y0(i,:)=y0(i,:)-oo_.Smoother.Trend.(v_name)(end-maximum_lag+1:end); %subtract trend, which is not subtracted by simult_
end end
end end
gend = options.nobs; gend = options_.nobs;
if isfield(oo.Smoother,'TrendCoeffs') if isfield(oo_.Smoother,'TrendCoeffs')
var_obs = options.varobs; var_obs = options_.varobs;
endo_names = M.endo_names; endo_names = M_.endo_names;
i_var_obs = []; i_var_obs = [];
trend_coeffs = []; trend_coeffs = [];
for i=1:length(var_obs) for i=1:length(var_obs)
tmp = strmatch(var_obs{i}, endo_names(i_var), 'exact'); tmp = strmatch(var_obs{i}, endo_names(i_var), 'exact');
trend_var_index = strmatch(var_obs{i}, M.endo_names, 'exact'); trend_var_index = strmatch(var_obs{i}, M_.endo_names, 'exact');
if ~isempty(tmp) if ~isempty(tmp)
i_var_obs = [ i_var_obs; tmp]; i_var_obs = [ i_var_obs; tmp];
trend_coeffs = [trend_coeffs; oo.Smoother.TrendCoeffs(trend_var_index)]; trend_coeffs = [trend_coeffs; oo_.Smoother.TrendCoeffs(trend_var_index)];
end end
end end
if ~isempty(trend_coeffs) if ~isempty(trend_coeffs)
trend = trend_coeffs*(first_obs+gend-1+(1-M.maximum_lag:horizon)); trend = trend_coeffs*(first_obs+gend-1+(1-M_.maximum_lag:horizon));
if options.prefilter if options_.prefilter
trend = trend - repmat(mean(trend_coeffs*[first_obs:first_obs+gend-1],2),1,horizon+1); %subtract mean trend trend = trend - repmat(mean(trend_coeffs*[first_obs:first_obs+gend-1],2),1,horizon+1); %subtract mean trend
end end
end end
else else
trend_coeffs=zeros(length(options.varobs),1); trend_coeffs=zeros(length(options_.varobs),1);
end end
otherwise otherwise
error('Wrong flag value') error('Wrong flag value')
end end
if M.exo_det_nbr == 0 if M_.exo_det_nbr == 0
if isequal(M.H,0) if isequal(M_.H,0)
[yf,int_width] = forcst(oo.dr,y0,horizon,var_list,M,oo,options); [yf,int_width] = forcst(oo_.dr,y0,horizon,var_list,M_,oo_,options_);
else else
[yf,int_width,int_width_ME] = forcst(oo.dr,y0,horizon,var_list,M,oo,options); [yf,int_width,int_width_ME] = forcst(oo_.dr,y0,horizon,var_list,M_,oo_,options_);
end end
else else
exo_det_length = size(oo.exo_det_simul,1)-M.maximum_lag; exo_det_length = size(oo_.exo_det_simul,1)-M_.maximum_lag;
if horizon > exo_det_length if horizon > exo_det_length
ex = zeros(horizon,M.exo_nbr); ex = zeros(horizon,M_.exo_nbr);
oo.exo_det_simul = [ oo.exo_det_simul;... oo_.exo_det_simul = [ oo_.exo_det_simul;...
repmat(oo.exo_det_steady_state',... repmat(oo_.exo_det_steady_state',...
horizon- ... horizon- ...
exo_det_length,1)]; exo_det_length,1)];
elseif horizon <= exo_det_length elseif horizon <= exo_det_length
ex = zeros(exo_det_length,M.exo_nbr); ex = zeros(exo_det_length,M_.exo_nbr);
end end
if options.linear if options_.linear
iorder = 1; iorder = 1;
else else
iorder = options.order; iorder = options_.order;
end end
if isequal(M.H,0) if isequal(M_.H,0)
[yf,int_width] = simultxdet(y0,ex,oo.exo_det_simul,... [yf,int_width] = simultxdet(y0,ex,oo_.exo_det_simul,...
iorder,var_list,M,oo,options); iorder,var_list,M_,oo_,options_);
else else
[yf,int_width,int_width_ME] = simultxdet(y0,ex,oo.exo_det_simul,... [yf,int_width,int_width_ME] = simultxdet(y0,ex,oo_.exo_det_simul,...
iorder,var_list,M,oo,options); iorder,var_list,M_,oo_,options_);
end end
end end
@ -185,13 +185,13 @@ if ~isscalar(trend) %add trend back to forecast
yf(i_var_obs,:) = yf(i_var_obs,:) + trend; yf(i_var_obs,:) = yf(i_var_obs,:) + trend;
end end
if options.loglinear if options_.loglinear
if options.prefilter == 1 %subtract steady state and add mean for observables if options_.prefilter == 1 %subtract steady state and add mean for observables
yf(i_var_obs,:)=yf(i_var_obs,:)-repmat(log(oo.dr.ys(i_var_obs)),1,horizon+M.maximum_lag)+ repmat(mean_varobs,1,horizon+M.maximum_lag); yf(i_var_obs,:)=yf(i_var_obs,:)-repmat(log(oo_.dr.ys(i_var_obs)),1,horizon+M_.maximum_lag)+ repmat(mean_varobs,1,horizon+M_.maximum_lag);
end end
else else
if options.prefilter == 1 %subtract steady state and add mean for observables if options_.prefilter == 1 %subtract steady state and add mean for observables
yf(i_var_obs,:)=yf(i_var_obs,:)-repmat(oo.dr.ys(i_var_obs),1,horizon+M.maximum_lag)+ repmat(mean_varobs,1,horizon+M.maximum_lag); yf(i_var_obs,:)=yf(i_var_obs,:)-repmat(oo_.dr.ys(i_var_obs),1,horizon+M_.maximum_lag)+ repmat(mean_varobs,1,horizon+M_.maximum_lag);
end end
end end
@ -200,17 +200,17 @@ for i=1:n_var
forecast.Mean.(vname) = yf(i,maximum_lag+(1:horizon))'; forecast.Mean.(vname) = yf(i,maximum_lag+(1:horizon))';
forecast.HPDinf.(vname)= yf(i,maximum_lag+(1:horizon))' - int_width(1:horizon,i); forecast.HPDinf.(vname)= yf(i,maximum_lag+(1:horizon))' - int_width(1:horizon,i);
forecast.HPDsup.(vname) = yf(i,maximum_lag+(1:horizon))' + int_width(1:horizon,i); forecast.HPDsup.(vname) = yf(i,maximum_lag+(1:horizon))' + int_width(1:horizon,i);
if ~isequal(M.H,0) && ismember(var_list{i},options.varobs) if ~isequal(M_.H,0) && ismember(var_list{i},options_.varobs)
forecast.HPDinf_ME.(vname)= yf(i,maximum_lag+(1:horizon))' - int_width_ME(1:horizon,i); forecast.HPDinf_ME.(vname)= yf(i,maximum_lag+(1:horizon))' - int_width_ME(1:horizon,i);
forecast.HPDsup_ME.(vname) = yf(i,maximum_lag+(1:horizon))' + int_width_ME(1:horizon,i); forecast.HPDsup_ME.(vname) = yf(i,maximum_lag+(1:horizon))' + int_width_ME(1:horizon,i);
end end
end end
for i=1:M.exo_det_nbr for i=1:M_.exo_det_nbr
forecast.Exogenous.(M.exo_det_names{i}) = oo.exo_det_simul(maximum_lag+(1:horizon),i); forecast.Exogenous.(M_.exo_det_names{i}) = oo_.exo_det_simul(maximum_lag+(1:horizon),i);
end end
if ~options.nograph if ~options_.nograph
oo.forecast = forecast; oo_.forecast = forecast;
forecast_graphs(var_list, M, oo, options) forecast_graphs(var_list, M_, oo_, options_)
end end

6
matlab/endogenous_prior.m
View File

@ -1,11 +1,11 @@
function [lnpriormom] = endogenous_prior(data,dataset_info, Pstar,BayesInfo,H) function [lnpriormom] = endogenous_prior(data,dataset_info, Pstar,bayestopt_,H)
% Computes the endogenous log prior addition to the initial prior % Computes the endogenous log prior addition to the initial prior
% %
% INPUTS % INPUTS
% data [double] n*T vector of data observations % data [double] n*T vector of data observations
% dataset_info [structure] various information about the dataset % dataset_info [structure] various information about the dataset
% Pstar [double] k*k matrix of % Pstar [double] k*k matrix of
% BayesInfo [structure] % bayestopt_ [structure]
% %
% OUTPUTS % OUTPUTS
% lnpriormom [double] scalar of log endogenous prior value % lnpriormom [double] scalar of log endogenous prior value
@ -80,7 +80,7 @@ Shat=C0 +(1-1/(2+1))*(C1+C1')...
+(1-2/(2+1))*(C2+C2'); +(1-2/(2+1))*(C2+C2');
% Model variances below: % Model variances below:
mf=BayesInfo.mf1; mf=bayestopt_.mf1;
II=eye(size(Pstar,2)); II=eye(size(Pstar,2));
Z=II(mf,:); Z=II(mf,:);
% This is Ftheta, variance of model variables, given param vector theta: % This is Ftheta, variance of model variables, given param vector theta:

36
matlab/ep/ep_accuracy_check.m
View File

@ -1,15 +1,33 @@
function e = ep_accuracy_check(M,options,oo) function e = ep_accuracy_check(M_,options_,oo_)
% e = ep_accuracy_check(M_,options_,oo_)
endo_simul = oo.endo_simul; % Copyright © 2016-2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
endo_simul = oo_.endo_simul;
n = size(endo_simul,2); n = size(endo_simul,2);
[initialconditions, innovations, pfm, ep, verbosity, options, oo] = ... [initialconditions, innovations, pfm, ep, verbosity, options_, oo_] = ...
extended_path_initialization([], n-1, [], options, M, oo); extended_path_initialization([], n-1, [], options_, M_, oo_);
options.ep.accuracy.stochastic.order = options.ep.stochastic.order; options_.ep.accuracy.stochastic.order = options_.ep.stochastic.order;
[nodes,weights,nnodes] = setup_integration_nodes(options.ep.accuracy,pfm); [nodes,weights,nnodes] = setup_integration_nodes(options_.ep.accuracy,pfm);
e = zeros(M.endo_nbr,n); e = zeros(M_.endo_nbr,n);
for i=1:n for i=1:n
e(:,i) = euler_equation_error(endo_simul(:,i),oo.exo_simul, ... e(:,i) = euler_equation_error(endo_simul(:,i),oo_.exo_simul, ...
innovations, M, options,oo,pfm,nodes,weights); innovations, M_, options_,oo_,pfm,nodes,weights);
end end

6
matlab/ep/ep_notes.org
View File

@ -55,9 +55,9 @@ solve_algo:
stack_solve_algo: stack_solve_algo:
ut: (unscented free parameter) ut: (unscented free parameter)
pfm.stochastic_order = DynareOptions.ep.stochastic.order; pfm.stochastic_order = options_.ep.stochastic.order;
pfm.periods = DynareOptions.ep.periods; pfm.periods = options_.ep.periods;
pfm.verbose = DynareOptions.ep.verbosity; pfm.verbose = options_.ep.verbosity;
* in extended_path_core, one passes options.ep and individual options * in extended_path_core, one passes options.ep and individual options

35
matlab/ep/euler_equation_error.m
View File

@ -1,6 +1,7 @@
function e = euler_equation_error(y0,x,innovations,M,options,oo,pfm,nodes,weights) function e = euler_equation_error(y0,x,innovations,M_,options_,oo_,pfm,nodes,weights)
% e = euler_equation_error(y0,x,innovations,M_,options_,oo_,pfm,nodes,weights)
% Copyright © 2016-2020 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -17,34 +18,34 @@ function e = euler_equation_error(y0,x,innovations,M,options,oo,pfm,nodes,weight
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
dynamic_model = str2func([M.fname '.dynamic']); dynamic_model = str2func([M_.fname '.dynamic']);
ep = options.ep; ep = options_.ep;
[y1, info_convergence, endogenousvariablespaths] = extended_path_core(ep.periods, ... [y1, info_convergence, endogenousvariablespaths] = extended_path_core(ep.periods, ...
M.endo_nbr, M.exo_nbr, ... M_.endo_nbr, M_.exo_nbr, ...
innovations.positive_var_indx, ... innovations.positive_var_indx, ...
x, ep.init, y0, oo.steady_state, ... x, ep.init, y0, oo_.steady_state, ...
0, ... 0, ...
ep.stochastic.order, M, ... ep.stochastic.order, M_, ...
pfm, ep.stochastic.algo, ... pfm, ep.stochastic.algo, ...
ep.solve_algo, ... ep.solve_algo, ...
ep.stack_solve_algo, ... ep.stack_solve_algo, ...
options.lmmcp, options, oo, ... options_.lmmcp, options_, oo_, ...
[]); []);
i_pred = find(M.lead_lag_incidence(1,:)); i_pred = find(M_.lead_lag_incidence(1,:));
i_fwrd = find(M.lead_lag_incidence(3,:)); i_fwrd = find(M_.lead_lag_incidence(3,:));
x1 = [x(2:end,:); zeros(1,M.exo_nbr)]; x1 = [x(2:end,:); zeros(1,M_.exo_nbr)];
for i=1:length(nodes) for i=1:length(nodes)
x2 = x1; x2 = x1;
x2(2,:) = x2(2,:) + nodes(i,:); x2(2,:) = x2(2,:) + nodes(i,:);
[y2, info_convergence, endogenousvariablespaths] = ... [y2, info_convergence, endogenousvariablespaths] = ...
extended_path_core(ep.periods, M.endo_nbr, M.exo_nbr, ... extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, ...
innovations.positive_var_indx, x2, ep.init, ... innovations.positive_var_indx, x2, ep.init, ...
y1, oo.steady_state, 0, ... y1, oo_.steady_state, 0, ...
ep.stochastic.order, M, pfm, ep.stochastic.algo, ... ep.stochastic.order, M_, pfm, ep.stochastic.algo, ...
ep.solve_algo, ep.stack_solve_algo, options.lmmcp, ... ep.solve_algo, ep.stack_solve_algo, options_.lmmcp, ...
options, oo, []); options_, oo_, []);
z = [y0(i_pred); y1; y2(i_fwrd)]; z = [y0(i_pred); y1; y2(i_fwrd)];
res(:,i) = dynamic_model(z,x,M.params,oo.steady_state,2); res(:,i) = dynamic_model(z,x,M_.params,oo_.steady_state,2);
end end
e = res*weights; e = res*weights;

48
matlab/ep/extended_path.m
View File

@ -1,5 +1,5 @@
function [ts, DynareResults] = extended_path(initialconditions, samplesize, exogenousvariables, DynareOptions, DynareModel, DynareResults) function [ts, oo_] = extended_path(initialconditions, samplesize, exogenousvariables, options_, M_, oo_)
% [ts, oo_] = extended_path(initialconditions, samplesize, exogenousvariables, options_, M_, oo_)
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time % Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
% series of size T is obtained by solving T perfect foresight models. % series of size T is obtained by solving T perfect foresight models.
% %
@ -7,9 +7,9 @@ function [ts, DynareResults] = extended_path(initialconditions, samplesize, exog
% o initialconditions [double] m*1 array, where m is the number of endogenous variables in the model. % o initialconditions [double] m*1 array, where m is the number of endogenous variables in the model.
% o samplesize [integer] scalar, size of the sample to be simulated. % o samplesize [integer] scalar, size of the sample to be simulated.
% o exogenousvariables [double] T*n array, values for the structural innovations. % o exogenousvariables [double] T*n array, values for the structural innovations.
% o DynareOptions [struct] options_ % o options_ [struct] options_
% o DynareModel [struct] M_ % o M_ [struct] Dynare's model structure
% o DynareResults [struct] oo_ % o oo_ [struct] Dynare's results structure
% %
% OUTPUTS % OUTPUTS
% o ts [dseries] m*samplesize array, the simulations. % o ts [dseries] m*samplesize array, the simulations.
@ -36,13 +36,13 @@ function [ts, DynareResults] = extended_path(initialconditions, samplesize, exog
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[initialconditions, innovations, pfm, ep, verbosity, DynareOptions, DynareResults] = ... [initialconditions, innovations, pfm, ep, verbosity, options_, oo_] = ...
extended_path_initialization(initialconditions, samplesize, exogenousvariables, DynareOptions, DynareModel, DynareResults); extended_path_initialization(initialconditions, samplesize, exogenousvariables, options_, M_, oo_);
[shocks, spfm_exo_simul, innovations, DynareResults] = extended_path_shocks(innovations, ep, exogenousvariables, samplesize,DynareModel,DynareOptions,DynareResults); [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, samplesize,M_,options_,oo_);
% Initialize the matrix for the paths of the endogenous variables. % Initialize the matrix for the paths of the endogenous variables.
endogenous_variables_paths = NaN(DynareModel.endo_nbr,samplesize+1); endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
endogenous_variables_paths(:,1) = initialconditions; endogenous_variables_paths(:,1) = initialconditions;
% Set waitbar (graphic or text mode) % Set waitbar (graphic or text mode)
@ -63,18 +63,18 @@ while (t <= samplesize)
if t>2 if t>2
% Set initial guess for the solver (using the solution of the % Set initial guess for the solver (using the solution of the
% previous period problem). % previous period problem).
initialguess = [endogenousvariablespaths(:, 2:end), DynareResults.steady_state]; initialguess = [endogenousvariablespaths(:, 2:end), oo_.steady_state];
else else
initialguess = []; initialguess = [];
end end
[endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(ep.periods, DynareModel.endo_nbr, DynareModel.exo_nbr, innovations.positive_var_indx, ... [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations.positive_var_indx, ...
spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ... spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
DynareResults.steady_state, ... oo_.steady_state, ...
verbosity, ep.stochastic.order, ... verbosity, ep.stochastic.order, ...
DynareModel, pfm, ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ... M_, pfm, ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
DynareOptions.lmmcp, ... options_.lmmcp, ...
DynareOptions, ... options_, ...
DynareResults, initialguess); oo_, initialguess);
if ~info_convergence if ~info_convergence
msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s)!', int2str(t)); msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s)!', int2str(t));
warning(msg) warning(msg)
@ -86,13 +86,13 @@ end % (while) loop over t
dyn_waitbar_close(hh_fig); dyn_waitbar_close(hh_fig);
% Set the initial period. % Set the initial period.
if isdates(DynareOptions.initial_period) if isdates(options_.initial_period)
if ischar(DynareOptions.initial_period) if ischar(options_.initial_period)
initial_period = dates(DynareOptions.initial_period); initial_period = dates(options_.initial_period);
else else
initial_period = DynareOptions.initial_period; initial_period = options_.initial_period;
end end
elseif isnan(DynareOptions.initial_period) elseif isnan(options_.initial_period)
initial_period = dates(1,1); initial_period = dates(1,1);
else else
error('Type of option initial_period is wrong.') error('Type of option initial_period is wrong.')
@ -104,11 +104,11 @@ if any(isnan(endogenous_variables_paths(:)))
nn = size(endogenous_variables_paths, 1); nn = size(endogenous_variables_paths, 1);
endogenous_variables_paths = reshape(endogenous_variables_paths(sl), nn, length(sl)/nn); endogenous_variables_paths = reshape(endogenous_variables_paths(sl), nn, length(sl)/nn);
end end
ts = dseries(transpose(endogenous_variables_paths), initial_period, DynareModel.endo_names); ts = dseries(transpose(endogenous_variables_paths), initial_period, M_.endo_names);
DynareResults.endo_simul = transpose(ts.data); oo_.endo_simul = transpose(ts.data);
assignin('base', 'Simulated_time_series', ts); assignin('base', 'Simulated_time_series', ts);
if ~nargout || nargout<2 if ~nargout || nargout<2
assignin('base', 'oo_', DynareResults); assignin('base', 'oo_', oo_);
end end

38
matlab/ep/extended_path_core.m
View File

@ -1,8 +1,8 @@
function [y, info_convergence, endogenousvariablespaths] = extended_path_core(periods,endo_nbr,exo_nbr,positive_var_indx, ... function [y, info_convergence, endogenousvariablespaths] = extended_path_core(periods,endo_nbr,exo_nbr,positive_var_indx, ...
exo_simul,init,initial_conditions,... exo_simul,init,initial_conditions,...
steady_state, ... steady_state, ...
debug,order,M,pfm,algo,solve_algo,stack_solve_algo,... debug,order,M_,pfm,algo,solve_algo,stack_solve_algo,...
olmmcp,options,oo,initialguess) olmmcp,options_,oo_,initialguess)
% Copyright © 2016-2023 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
@ -21,10 +21,10 @@ function [y, info_convergence, endogenousvariablespaths] = extended_path_core(pe
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
ep = options.ep; ep = options_.ep;
if init% Compute first order solution (Perturbation)... if init% Compute first order solution (Perturbation)...
endo_simul = simult_(M,options,initial_conditions,oo.dr,exo_simul(2:end,:),1); endo_simul = simult_(M_,options_,initial_conditions,oo_.dr,exo_simul(2:end,:),1);
else else
if nargin==19 && ~isempty(initialguess) if nargin==19 && ~isempty(initialguess)
% Note that the first column of initialguess should be equal to initial_conditions. % Note that the first column of initialguess should be equal to initial_conditions.
@ -34,29 +34,29 @@ else
end end
end end
oo.endo_simul = endo_simul; oo_.endo_simul = endo_simul;
if debug if debug
save ep_test_1.mat endo_simul exo_simul save ep_test_1.mat endo_simul exo_simul
end end
if options.bytecode && order > 0 if options_.bytecode && order > 0
error('Option order > 0 of extended_path command is not compatible with bytecode option.') error('Option order > 0 of extended_path command is not compatible with bytecode option.')
end end
if options.block && order > 0 if options_.block && order > 0
error('Option order > 0 of extended_path command is not compatible with block option.') error('Option order > 0 of extended_path command is not compatible with block option.')
end end
if order == 0 if order == 0
options.periods = periods; options_.periods = periods;
options.block = pfm.block; options_.block = pfm.block;
oo.endo_simul = endo_simul; oo_.endo_simul = endo_simul;
oo.exo_simul = exo_simul; oo_.exo_simul = exo_simul;
oo.steady_state = steady_state; oo_.steady_state = steady_state;
options.lmmcp = olmmcp; options_.lmmcp = olmmcp;
options.solve_algo = solve_algo; options_.solve_algo = solve_algo;
options.stack_solve_algo = stack_solve_algo; options_.stack_solve_algo = stack_solve_algo;
[endogenousvariablespaths, info_convergence] = perfect_foresight_solver_core(oo.endo_simul, oo.exo_simul, oo.steady_state, oo.exo_steady_state, M, options); [endogenousvariablespaths, info_convergence] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
else else
switch(algo) switch(algo)
case 0 case 0
@ -64,13 +64,13 @@ else
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order); solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1 case 1
[flag, endogenousvariablespaths] = ... [flag, endogenousvariablespaths] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order); solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);
end end
info_convergence = ~flag; info_convergence = ~flag;
end end
if ~info_convergence && ~options.no_homotopy if ~info_convergence && ~options_.no_homotopy
[info_convergence, endogenousvariablespaths] = extended_path_homotopy(endo_simul, exo_simul, M, options, oo, pfm, ep, order, algo, 2, debug); [info_convergence, endogenousvariablespaths] = extended_path_homotopy(endo_simul, exo_simul, M_, options_, oo_, pfm, ep, order, algo, 2, debug);
end end
if info_convergence if info_convergence

20
matlab/ep/extended_path_homotopy.m
View File

@ -1,4 +1,4 @@
function [info_convergence, endo_simul] = extended_path_homotopy(endo_simul, exo_simul, M, options, oo, pfm, ep, order, algo, method, debug) function [info_convergence, endo_simul] = extended_path_homotopy(endo_simul, exo_simul, M_, options_, oo_, pfm, ep, order, algo, method, debug)
% Copyright © 2016-2023 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
@ -31,11 +31,11 @@ if ismember(method, [1, 2])
oldweight = weight; oldweight = weight;
while noconvergence while noconvergence
iteration = iteration + 1; iteration = iteration + 1;
oo.endo_simul = endo_simul; oo_.endo_simul = endo_simul;
oo.endo_simul(:,1) = oo.steady_state + weight*(endo_simul0(:,1) - oo.steady_state); oo_.endo_simul(:,1) = oo_.steady_state + weight*(endo_simul0(:,1) - oo_.steady_state);
oo.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo.exo_steady_state)); oo_.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo_.exo_steady_state));
if order==0 if order==0
[endo_simul_new, success] = perfect_foresight_solver_core(oo.endo_simul, oo.exo_simul, oo.steady_state, oo.exo_steady_state, M, options); [endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
else else
switch(algo) switch(algo)
case 0 case 0
@ -43,7 +43,7 @@ if ismember(method, [1, 2])
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order); solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1 case 1
[flag, endo_simul_new] = ... [flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order); solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);
end end
end end
if isequal(order, 0) if isequal(order, 0)
@ -99,10 +99,10 @@ if isequal(method, 3) || (isequal(method, 2) && noconvergence)
nweights = length(weights); nweights = length(weights);
while noconvergence while noconvergence
weight = weights(index); weight = weights(index);
oo.endo_simul = endo_simul; oo_.endo_simul = endo_simul;
oo.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo.exo_steady_state)); oo_.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo_.exo_steady_state));
if order==0 if order==0
[endo_simul_new, success] = perfect_foresight_solver_core(oo.endo_simul, oo.exo_simul, oo.steady_state, oo.exo_steady_state, M, options); [endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
else else
switch(algo) switch(algo)
case 0 case 0
@ -110,7 +110,7 @@ if isequal(method, 3) || (isequal(method, 2) && noconvergence)
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order); solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1 case 1
[flag, endo_simul_new] = ... [flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order); solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);
end end
end end
if isequal(order, 0) if isequal(order, 0)

70
matlab/ep/extended_path_initialization.m
View File

@ -1,14 +1,14 @@
function [initial_conditions, innovations, pfm, ep, verbosity, DynareOptions, DynareResults] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, DynareOptions, DynareModel, DynareResults) function [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, options_, M_, oo_)
% [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, options_, M_, oo_)
% Initialization of the extended path routines. % Initialization of the extended path routines.
% %
% INPUTS % INPUTS
% o initial_conditions [double] m*1 array, where m is the number of endogenous variables in the model. % o initial_conditions [double] m*1 array, where m is the number of endogenous variables in the model.
% o sample_size [integer] scalar, size of the sample to be simulated. % o sample_size [integer] scalar, size of the sample to be simulated.
% o exogenousvariables [double] T*n array, values for the structural innovations. % o exogenousvariables [double] T*n array, values for the structural innovations.
% o DynareOptions [struct] options_ % o options_ [struct] Dynare's options structure
% o DynareModel [struct] M_ % o M_ [struct] Dynare's model structure
% o DynareResults [struct] oo_ % o oo_ [struct] Dynare's result structure
% %
% OUTPUTS % OUTPUTS
% %
@ -16,7 +16,7 @@ function [initial_conditions, innovations, pfm, ep, verbosity, DynareOptions, Dy
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% Copyright © 2016-2020 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -33,29 +33,29 @@ function [initial_conditions, innovations, pfm, ep, verbosity, DynareOptions, Dy
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
ep = DynareOptions.ep; ep = options_.ep;
% Set verbosity levels. % Set verbosity levels.
DynareOptions.verbosity = ep.verbosity; options_.verbosity = ep.verbosity;
verbosity = ep.verbosity+ep.debug; verbosity = ep.verbosity+ep.debug;
% Set maximum number of iterations for the deterministic solver. % Set maximum number of iterations for the deterministic solver.
DynareOptions.simul.maxit = ep.maxit; options_.simul.maxit = ep.maxit;
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver % Prepare a structure needed by the matlab implementation of the perfect foresight model solver
pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel, DynareOptions, DynareResults); pfm = setup_stochastic_perfect_foresight_model_solver(M_, options_, oo_);
% Check that the user did not use varexo_det % Check that the user did not use varexo_det
if DynareModel.exo_det_nbr~=0 if M_.exo_det_nbr~=0
error('Extended path does not support varexo_det.') error('Extended path does not support varexo_det.')
end end
% Set default initial conditions. % Set default initial conditions.
if isempty(initial_conditions) if isempty(initial_conditions)
if isempty(DynareModel.endo_histval) if isempty(M_.endo_histval)
initial_conditions = DynareResults.steady_state; initial_conditions = oo_.steady_state;
else else
initial_conditions = DynareModel.endo_histval; initial_conditions = M_.endo_histval;
end end
end end
@ -64,10 +64,10 @@ pfm.periods = ep.periods;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny); pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
pfm.block = DynareOptions.block; pfm.block = options_.block;
% Set the algorithm for the perfect foresight solver % Set the algorithm for the perfect foresight solver
DynareOptions.stack_solve_algo = ep.stack_solve_algo; options_.stack_solve_algo = ep.stack_solve_algo;
% Compute the first order reduced form if needed. % Compute the first order reduced form if needed.
% %
@ -76,20 +76,20 @@ DynareOptions.stack_solve_algo = ep.stack_solve_algo;
dr = struct(); dr = struct();
if ep.init if ep.init
DynareOptions.order = 1; options_.order = 1;
DynareResults.dr=set_state_space(dr,DynareModel); oo_.dr=set_state_space(dr,M_);
[DynareResults.dr,Info,DynareModel.params] = resol(0,DynareModel,DynareOptions,DynareResults.dr,DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state); [oo_.dr,Info,M_.params] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
end end
% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks) % Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
DynareOptions.minimal_solving_period = DynareOptions.ep.periods; options_.minimal_solving_period = options_.ep.periods;
% Set the covariance matrix of the structural innovations. % Set the covariance matrix of the structural innovations.
if isempty(exogenousvariables) if isempty(exogenousvariables)
innovations = struct(); innovations = struct();
innovations.positive_var_indx = find(diag(DynareModel.Sigma_e)>0); innovations.positive_var_indx = find(diag(M_.Sigma_e)>0);
innovations.effective_number_of_shocks = length(innovations.positive_var_indx); innovations.effective_number_of_shocks = length(innovations.positive_var_indx);
innovations.covariance_matrix = DynareModel.Sigma_e(innovations.positive_var_indx,innovations.positive_var_indx); innovations.covariance_matrix = M_.Sigma_e(innovations.positive_var_indx,innovations.positive_var_indx);
innovations.covariance_matrix_upper_cholesky = chol(innovations.covariance_matrix); innovations.covariance_matrix_upper_cholesky = chol(innovations.covariance_matrix);
else else
innovations = struct(); innovations = struct();
@ -97,26 +97,26 @@ end
% Set seed. % Set seed.
if ep.set_dynare_seed_to_default if ep.set_dynare_seed_to_default
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
end end
% hybrid correction % hybrid correction
pfm.hybrid_order = ep.stochastic.hybrid_order; pfm.hybrid_order = ep.stochastic.hybrid_order;
if pfm.hybrid_order if pfm.hybrid_order
DynareResults.dr = set_state_space(DynareResults.dr, DynareModel); oo_.dr = set_state_space(oo_.dr, M_);
options = DynareOptions; options = options_;
options.order = pfm.hybrid_order; options.order = pfm.hybrid_order;
[pfm.dr, DynareModel.params] = resol(0, DynareModel, options, DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state); [pfm.dr, M_.params] = resol(0, M_, options, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
else else
pfm.dr = []; pfm.dr = [];
end end
% number of nonzero derivatives % number of nonzero derivatives
pfm.nnzA = DynareModel.NNZDerivatives(1); pfm.nnzA = M_.NNZDerivatives(1);
% setting up integration nodes if order > 0 % setting up integration nodes if order > 0
if ep.stochastic.order > 0 if ep.stochastic.order > 0
[nodes,weights,nnodes] = setup_integration_nodes(DynareOptions.ep,pfm); [nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm);
pfm.nodes = nodes; pfm.nodes = nodes;
pfm.weights = weights; pfm.weights = weights;
pfm.nnodes = nnodes; pfm.nnodes = nnodes;
@ -127,12 +127,12 @@ else
end end
% set boundaries if mcp % set boundaries if mcp
[lb,ub,pfm.eq_index] = get_complementarity_conditions(DynareModel, DynareOptions.ramsey_policy); [lb,ub,pfm.eq_index] = get_complementarity_conditions(M_, options_.ramsey_policy);
if DynareOptions.ep.solve_algo == 10 if options_.ep.solve_algo == 10
DynareOptions.lmmcp.lb = repmat(lb,block_nbr,1); options_.lmmcp.lb = repmat(lb,block_nbr,1);
DynareOptions.lmmcp.ub = repmat(ub,block_nbr,1); options_.lmmcp.ub = repmat(ub,block_nbr,1);
elseif DynareOptions.ep.solve_algo == 11 elseif options_.ep.solve_algo == 11
DynareOptions.mcppath.lb = repmat(lb,block_nbr,1); options_.mcppath.lb = repmat(lb,block_nbr,1);
DynareOptions.mcppath.ub = repmat(ub,block_nbr,1); options_.mcppath.ub = repmat(ub,block_nbr,1);
end end
pfm.block_nbr = block_nbr; pfm.block_nbr = block_nbr;

44
matlab/ep/extended_path_mc.m
View File

@ -1,15 +1,15 @@
function Simulations = extended_path_mc(initialconditions, samplesize, replic, exogenousvariables, DynareOptions, DynareModel, DynareResults) function Simulations = extended_path_mc(initialconditions, samplesize, replic, exogenousvariables, options_, M_, oo_)
% Simulations = extended_path_mc(initialconditions, samplesize, replic, exogenousvariables, options_, M_, oo_)
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time % Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
% series of size T is obtained by solving T perfect foresight models. % series of size T is obtained by solving T perfect foresight models.
% %
% INPUTS % INPUTS
% o initialconditions [double] m*1 array, where m is the number of endogenous variables in the model. % o initialconditions [double] m*1 array, where m is the number of endogenous variables in the model.
% o samplesize [integer] scalar, size of the sample to be simulated. % o samplesize [integer] scalar, size of the sample to be simulated.
% o exogenousvariables [double] T*n array, values for the structural innovations. % o exogenousvariables [double] T*n array, values for the structural innovations.
% o DynareOptions [struct] options_ % o options_ [struct] Dynare's options structure
% o DynareModel [struct] M_ % o M_ [struct] Dynare's model structure
% o DynareResults [struct] oo_ % o oo_ [struct] Dynare's results structure
% %
% OUTPUTS % OUTPUTS
% o ts [dseries] m*samplesize array, the simulations. % o ts [dseries] m*samplesize array, the simulations.
@ -19,7 +19,7 @@ function Simulations = extended_path_mc(initialconditions, samplesize, replic, e
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% Copyright © 2016-2020 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -36,8 +36,8 @@ function Simulations = extended_path_mc(initialconditions, samplesize, replic, e
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[initialconditions, innovations, pfm, ep, verbosity, DynareOptions, DynareResults] = ... [initialconditions, innovations, pfm, ep, verbosity, options_, oo_] = ...
extended_path_initialization(initialconditions, samplesize, exogenousvariables, DynareOptions, DynareModel, DynareResults); extended_path_initialization(initialconditions, samplesize, exogenousvariables, options_, M_, oo_);
% Check the dimension of the first input argument % Check the dimension of the first input argument
if isequal(size(initialconditions, 2), 1) if isequal(size(initialconditions, 2), 1)
@ -68,9 +68,9 @@ if ep.parallel
% Use the Parallel toolbox. % Use the Parallel toolbox.
parfor i=1:replic parfor i=1:replic
innovations_ = innovations; innovations_ = innovations;
DynareResults_ = DynareResults; oo__ = oo_;
[shocks, spfm_exo_simul, innovations_, DynareResults_] = extended_path_shocks(innovations_, ep, exogenousvariables(:,:,i), samplesize, DynareModel, DynareOptions, DynareResults_); [shocks, spfm_exo_simul, innovations_, oo__] = extended_path_shocks(innovations_, ep, exogenousvariables(:,:,i), samplesize, M_, options_, oo__);
endogenous_variables_paths = NaN(DynareModel.endo_nbr,samplesize+1); endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
endogenous_variables_paths(:,1) = initialconditions(:,1); endogenous_variables_paths(:,1) = initialconditions(:,1);
exogenous_variables_paths = NaN(innovations_.effective_number_of_shocks,samplesize+1); exogenous_variables_paths = NaN(innovations_.effective_number_of_shocks,samplesize+1);
exogenous_variables_paths(:,1) = 0; exogenous_variables_paths(:,1) = 0;
@ -80,12 +80,12 @@ if ep.parallel
t = t+1; t = t+1;
spfm_exo_simul(2,:) = shocks(t-1,:); spfm_exo_simul(2,:) = shocks(t-1,:);
exogenous_variables_paths(:,t) = shocks(t-1,:); exogenous_variables_paths(:,t) = shocks(t-1,:);
[endogenous_variables_paths(:,t), info_convergence] = extended_path_core(ep.periods, DynareModel.endo_nbr, DynareModel.exo_nbr, innovations_.positive_var_indx, ... [endogenous_variables_paths(:,t), info_convergence] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations_.positive_var_indx, ...
spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ... spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
DynareResults_.steady_state, ... oo__.steady_state, ...
ep.verbosity, ep.stochastic.order, ... ep.verbosity, ep.stochastic.order, ...
DynareModel, pfm,ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ... M_, pfm,ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
DynareOptions.lmmcp, DynareOptions, DynareResults_); options_.lmmcp, options_, oo__);
if ~info_convergence if ~info_convergence
msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s, iteration %s)!', int2str(t), int2str(i)); msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s, iteration %s)!', int2str(t), int2str(i));
warning(msg) warning(msg)
@ -99,8 +99,8 @@ if ep.parallel
else else
% Sequential approach. % Sequential approach.
for i=1:replic for i=1:replic
[shocks, spfm_exo_simul, innovations, DynareResults] = extended_path_shocks(innovations, ep, exogenousvariables(:,:,i), samplesize, DynareModel, DynareOptions, DynareResults); [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables(:,:,i), samplesize, M_, options_, oo_);
endogenous_variables_paths = NaN(DynareModel.endo_nbr,samplesize+1); endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
endogenous_variables_paths(:,1) = initialconditions(:,1); endogenous_variables_paths(:,1) = initialconditions(:,1);
exogenous_variables_paths = NaN(innovations.effective_number_of_shocks,samplesize+1); exogenous_variables_paths = NaN(innovations.effective_number_of_shocks,samplesize+1);
exogenous_variables_paths(:,1) = 0; exogenous_variables_paths(:,1) = 0;
@ -109,12 +109,12 @@ else
t = t+1; t = t+1;
spfm_exo_simul(2,:) = shocks(t-1,:); spfm_exo_simul(2,:) = shocks(t-1,:);
exogenous_variables_paths(:,t) = shocks(t-1,:); exogenous_variables_paths(:,t) = shocks(t-1,:);
[endogenous_variables_paths(:,t), info_convergence] = extended_path_core(ep.periods, DynareModel.endo_nbr, DynareModel.exo_nbr, innovations.positive_var_indx, ... [endogenous_variables_paths(:,t), info_convergence] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations.positive_var_indx, ...
spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ... spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
DynareResults.steady_state, ... oo_.steady_state, ...
ep.verbosity, ep.stochastic.order, ... ep.verbosity, ep.stochastic.order, ...
DynareModel, pfm,ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ... M_, pfm,ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
DynareOptions.lmmcp, DynareOptions, DynareResults); options_.lmmcp, options_, oo_);
if ~info_convergence if ~info_convergence
msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s, iteration %s)!', int2str(t), int2str(i)); msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s, iteration %s)!', int2str(t), int2str(i));
warning(msg) warning(msg)

18
matlab/ep/extended_path_shocks.m
View File

@ -1,6 +1,6 @@
function [shocks, spfm_exo_simul, innovations, DynareResults] = extended_path_shocks(innovations, ep, exogenousvariables, sample_size,DynareModel,DynareOptions, DynareResults) function [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, sample_size,M_,options_, oo_)
% [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, sample_size,M_,options_, oo_)
% Copyright © 2016-2017 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -21,13 +21,13 @@ function [shocks, spfm_exo_simul, innovations, DynareResults] = extended_path_sh
if isempty(exogenousvariables) if isempty(exogenousvariables)
switch ep.innovation_distribution switch ep.innovation_distribution
case 'gaussian' case 'gaussian'
shocks = zeros(sample_size, DynareModel.exo_nbr); shocks = zeros(sample_size, M_.exo_nbr);
shocks(:,innovations.positive_var_indx) = transpose(transpose(innovations.covariance_matrix_upper_cholesky)*randn(innovations.effective_number_of_shocks,sample_size)); shocks(:,innovations.positive_var_indx) = transpose(transpose(innovations.covariance_matrix_upper_cholesky)*randn(innovations.effective_number_of_shocks,sample_size));
case 'calibrated' case 'calibrated'
options = DynareOptions; options = options_;
options.periods = options.ep.periods; options.periods = options.ep.periods;
oo = make_ex_(DynareModel,options,DynareResults); oo_local = make_ex_(M_,options,oo_);
shocks = oo.exo_simul(2:end,:); shocks = oo_local.exo_simul(2:end,:);
otherwise otherwise
error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!']) error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end end
@ -37,5 +37,5 @@ else
end end
% Copy the shocks in exo_simul % Copy the shocks in exo_simul
DynareResults.exo_simul = shocks; oo_.exo_simul = shocks;
spfm_exo_simul = repmat(DynareResults.exo_steady_state',ep.periods+2,1); spfm_exo_simul = repmat(oo_.exo_steady_state',ep.periods+2,1);

39
matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
View File

@ -1,6 +1,11 @@
function pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,DynareOptions,DynareOutput) function pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_)
% pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_)
% INPUTS
% o M_ [struct] Dynare's model structure
% o options_ [struct] Dynare's options structure
% o oo_ [struct] Dynare's results structure
% Copyright © 2013-2020 Dynare Team % Copyright © 2013-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -17,15 +22,15 @@ function pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,Dynar
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
pfm.lead_lag_incidence = DynareModel.lead_lag_incidence; pfm.lead_lag_incidence = M_.lead_lag_incidence;
pfm.ny = DynareModel.endo_nbr; pfm.ny = M_.endo_nbr;
pfm.Sigma = DynareModel.Sigma_e; pfm.Sigma = M_.Sigma_e;
if det(pfm.Sigma) > 0 if det(pfm.Sigma) > 0
pfm.Omega = chol(pfm.Sigma,'upper'); % Sigma = Omega'*Omega pfm.Omega = chol(pfm.Sigma,'upper'); % Sigma = Omega'*Omega
end end
pfm.number_of_shocks = length(pfm.Sigma); pfm.number_of_shocks = length(pfm.Sigma);
pfm.stochastic_order = DynareOptions.ep.stochastic.order; pfm.stochastic_order = options_.ep.stochastic.order;
pfm.max_lag = DynareModel.maximum_endo_lag; pfm.max_lag = M_.maximum_endo_lag;
if pfm.max_lag > 0 if pfm.max_lag > 0
pfm.nyp = nnz(pfm.lead_lag_incidence(1,:)); pfm.nyp = nnz(pfm.lead_lag_incidence(1,:));
pfm.iyp = find(pfm.lead_lag_incidence(1,:)>0); pfm.iyp = find(pfm.lead_lag_incidence(1,:)>0);
@ -35,7 +40,7 @@ else
end end
pfm.ny0 = nnz(pfm.lead_lag_incidence(pfm.max_lag+1,:)); pfm.ny0 = nnz(pfm.lead_lag_incidence(pfm.max_lag+1,:));
pfm.iy0 = find(pfm.lead_lag_incidence(pfm.max_lag+1,:)>0); pfm.iy0 = find(pfm.lead_lag_incidence(pfm.max_lag+1,:)>0);
if DynareModel.maximum_endo_lead if M_.maximum_endo_lead
pfm.nyf = nnz(pfm.lead_lag_incidence(pfm.max_lag+2,:)); pfm.nyf = nnz(pfm.lead_lag_incidence(pfm.max_lag+2,:));
pfm.iyf = find(pfm.lead_lag_incidence(pfm.max_lag+2,:)>0); pfm.iyf = find(pfm.lead_lag_incidence(pfm.max_lag+2,:)>0);
else else
@ -49,10 +54,10 @@ pfm.is = [pfm.nyp+1:pfm.ny+pfm.nyp];
pfm.isf = pfm.iyf+pfm.nyp; pfm.isf = pfm.iyf+pfm.nyp;
pfm.isf1 = [pfm.nyp+pfm.ny+1:pfm.nyf+pfm.nyp+pfm.ny+1]; pfm.isf1 = [pfm.nyp+pfm.ny+1:pfm.nyf+pfm.nyp+pfm.ny+1];
pfm.iz = [1:pfm.ny+pfm.nyp+pfm.nyf]; pfm.iz = [1:pfm.ny+pfm.nyp+pfm.nyf];
pfm.periods = DynareOptions.ep.periods; pfm.periods = options_.ep.periods;
pfm.steady_state = DynareOutput.steady_state; pfm.steady_state = oo_.steady_state;
pfm.params = DynareModel.params; pfm.params = M_.params;
if DynareModel.maximum_endo_lead if M_.maximum_endo_lead
pfm.i_cols_1 = nonzeros(pfm.lead_lag_incidence(pfm.max_lag+(1:2),:)'); pfm.i_cols_1 = nonzeros(pfm.lead_lag_incidence(pfm.max_lag+(1:2),:)');
pfm.i_cols_A1 = find(pfm.lead_lag_incidence(pfm.max_lag+(1:2),:)'); pfm.i_cols_A1 = find(pfm.lead_lag_incidence(pfm.max_lag+(1:2),:)');
else else
@ -66,9 +71,9 @@ else
end end
pfm.i_cols_j = 1:pfm.nd; pfm.i_cols_j = 1:pfm.nd;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny); pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
if ~DynareOptions.bytecode if ~options_.bytecode
pfm.dynamic_model = str2func([DynareModel.fname,'.dynamic']); pfm.dynamic_model = str2func([M_.fname,'.dynamic']);
end end
pfm.verbose = DynareOptions.ep.verbosity; pfm.verbose = options_.ep.verbosity;
pfm.maxit_ = DynareOptions.simul.maxit; pfm.maxit_ = options_.simul.maxit;
pfm.tolerance = DynareOptions.dynatol.f; pfm.tolerance = options_.dynatol.f;

6
matlab/get_equation_number_by_tag.m
View File

@ -1,15 +1,15 @@
function eqnumber = get_equation_number_by_tag(eqname, M_) function eqnumber = get_equation_number_by_tag(eqname, M_)
% eqnumber = get_equation_number_by_tag(eqname, M_)
% Translates an equation name into an equation number. % Translates an equation name into an equation number.
% %
% INPUTS % INPUTS
% - eqname [char] 1×n array, name of the equation. % - eqname [char] 1×n array, name of the equation.
% - DynareModel [struct] Structure describing the model, aka M_. % - M_ [struct] Structure describing the model
% %
% OUTPUTS % OUTPUTS
% - eqnumber [integer] Equation number. % - eqnumber [integer] Equation number.
% Copyright © 2018-2022 Dynare Team % Copyright © 2018-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %

12
matlab/get_lhs_and_rhs.m
View File

@ -1,10 +1,10 @@
function [lhs, rhs, json] = get_lhs_and_rhs(eqname, DynareModel, original, json) function [lhs, rhs, json] = get_lhs_and_rhs(eqname, M_, original, json)
% [lhs, rhs, json] = get_lhs_and_rhs(eqname, M_, original, json)
% Returns the left and right handsides of an equation. % Returns the left and right handsides of an equation.
% %
% INPUTS % INPUTS
% - eqname [char] Name of the equation. % - eqname [char] Name of the equation.
% - DynareModel [struct] Structure describing the current model (M_). % - M_ [struct] Structure describing the current model.
% - original [logical] fetch equation in modfile-original.json or modfile.json % - original [logical] fetch equation in modfile-original.json or modfile.json
% - json [char] content of the JSON file % - json [char] content of the JSON file
% %
@ -23,7 +23,7 @@ function [lhs, rhs, json] = get_lhs_and_rhs(eqname, DynareModel, original, json)
% [name='Phillips curve'] % [name='Phillips curve']
% pi = beta*pi(1) + slope*y + lam; % pi = beta*pi(1) + slope*y + lam;
% Copyright © 2018-2020 Dynare Team % Copyright © 2018-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -47,9 +47,9 @@ end
% Load JSON file if nargin<4 % Load JSON file if nargin<4
if nargin<4 if nargin<4
if original if original
json = loadjson_([DynareModel.fname filesep() 'model' filesep() 'json' filesep() 'modfile-original.json']); json = loadjson_([M_.fname filesep() 'model' filesep() 'json' filesep() 'modfile-original.json']);
else else
json = loadjson_([DynareModel.fname filesep() 'model' filesep() 'json' filesep() 'modfile.json']); json = loadjson_([M_.fname filesep() 'model' filesep() 'json' filesep() 'modfile.json']);
end end
end end

18
matlab/get_variables_and_parameters_in_equation.m
View File

@ -1,11 +1,11 @@
function [pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, DynareModel) function [pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, M_)
% [pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, M_)
% Returns the lists of parameters, endogenous variables and exogenous variables in an equation. % Returns the lists of parameters, endogenous variables and exogenous variables in an equation.
% %
% INPUTS % INPUTS
% - lhs [string] Left hand side of an equation. % - lhs [string] Left hand side of an equation.
% - rhs [string] Right hand side of an equation. % - rhs [string] Right hand side of an equation.
% - DynareModel [struct] Structure describing the current model (M_). % - M_ [struct] Structure describing the current model.
% %
% OUTPUTS % OUTPUTS
% - pnames [cell] Cell of row char arrays (p elements), names of the parameters. % - pnames [cell] Cell of row char arrays (p elements), names of the parameters.
@ -39,33 +39,33 @@ rhs_ = get_variables_and_parameters_in_expression(rhs);
lhs_ = get_variables_and_parameters_in_expression(lhs); lhs_ = get_variables_and_parameters_in_expression(lhs);
% Get list of parameters. % Get list of parameters.
pnames = DynareModel.param_names; pnames = M_.param_names;
pnames = intersect([rhs_, lhs_], pnames); pnames = intersect([rhs_, lhs_], pnames);
if nargout>1 if nargout>1
% Get list of endogenous variables. % Get list of endogenous variables.
enames = DynareModel.endo_names; enames = M_.endo_names;
enames = intersect([rhs_, lhs_], enames); enames = intersect([rhs_, lhs_], enames);
if nargout>2 if nargout>2
% Get list of exogenous variables % Get list of exogenous variables
xnames = DynareModel.exo_names; xnames = M_.exo_names;
xnames = intersect([rhs_,lhs_], xnames); xnames = intersect([rhs_,lhs_], xnames);
if nargout>3 if nargout>3
% Returns vector of indices for parameters endogenous and exogenous variables if required. % Returns vector of indices for parameters endogenous and exogenous variables if required.
p = length(pnames); p = length(pnames);
pid = zeros(p, 1); pid = zeros(p, 1);
for i = 1:p for i = 1:p
pid(i) = find(strcmp(pnames{i}, DynareModel.param_names)); pid(i) = find(strcmp(pnames{i}, M_.param_names));
end end
p = length(enames); p = length(enames);
eid = zeros(p, 1); eid = zeros(p, 1);
for i = 1:p for i = 1:p
eid(i) = find(strcmp(enames{i}, DynareModel.endo_names)); eid(i) = find(strcmp(enames{i}, M_.endo_names));
end end
p = length(xnames); p = length(xnames);
xid = zeros(p, 1); xid = zeros(p, 1);
for i = 1:p for i = 1:p
xid(i) = find(strcmp(xnames{i}, DynareModel.exo_names)); xid(i) = find(strcmp(xnames{i}, M_.exo_names));
end end
end end
end end

107
matlab/gsa/map_calibration.m
View File

@ -1,5 +1,6 @@
function map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResults, EstimatedParameters, BayesInfo) function map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_, bayestopt_)
% map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResults, EstimatedParameters, BayesInfo) % map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_, bayestopt_)
% Written by Marco Ratto % Written by Marco Ratto
% Joint Research Centre, The European Commission, % Joint Research Centre, The European Commission,
@ -23,33 +24,33 @@ function map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResult
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
fname_ = Model.fname; fname_ = M_.fname;
np = EstimatedParameters.np; np = estim_params_.np;
nshock = EstimatedParameters.nvx + EstimatedParameters.nvn + EstimatedParameters.ncx + EstimatedParameters.ncn; nshock = estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
pnames=cell(np,1); pnames=cell(np,1);
pnames_tex=cell(np,1); pnames_tex=cell(np,1);
for jj=1:np for jj=1:np
if DynareOptions.TeX if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj, DynareOptions.TeX, Model, EstimatedParameters, DynareOptions); [param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj, options_.TeX, M_, estim_params_, options_);
pnames_tex{jj,1} = strrep(param_name_tex_temp,'$',''); pnames_tex{jj,1} = strrep(param_name_tex_temp,'$','');
pnames{jj,1} = param_name_temp; pnames{jj,1} = param_name_temp;
else else
param_name_temp = get_the_name(nshock+jj, DynareOptions.TeX, Model, EstimatedParameters, DynareOptions); param_name_temp = get_the_name(nshock+jj, options_.TeX, M_, estim_params_, options_);
pnames{jj,1} = param_name_temp; pnames{jj,1} = param_name_temp;
end end
end end
pvalue_ks = DynareOptions.opt_gsa.pvalue_ks; pvalue_ks = options_.opt_gsa.pvalue_ks;
indx_irf = []; indx_irf = [];
indx_moment = []; indx_moment = [];
init = ~DynareOptions.opt_gsa.load_stab; init = ~options_.opt_gsa.load_stab;
options_mcf.pvalue_ks = DynareOptions.opt_gsa.pvalue_ks; options_mcf.pvalue_ks = options_.opt_gsa.pvalue_ks;
options_mcf.pvalue_corr = DynareOptions.opt_gsa.pvalue_corr; options_mcf.pvalue_corr = options_.opt_gsa.pvalue_corr;
options_mcf.alpha2 = DynareOptions.opt_gsa.alpha2_stab; options_mcf.alpha2 = options_.opt_gsa.alpha2_stab;
options_mcf.param_names = pnames; options_mcf.param_names = pnames;
if DynareOptions.TeX if options_.TeX
options_mcf.param_names_tex = pnames_tex; options_mcf.param_names_tex = pnames_tex;
end end
options_mcf.fname_ = fname_; options_mcf.fname_ = fname_;
@ -58,17 +59,17 @@ options_mcf.OutputDirectoryName = OutputDirectoryName;
skipline() skipline()
disp('Sensitivity analysis for calibration criteria') disp('Sensitivity analysis for calibration criteria')
if DynareOptions.opt_gsa.ppost if options_.opt_gsa.ppost
filetoload=dir([Model.dname filesep 'metropolis' filesep fname_ '_param_irf*.mat']); filetoload=dir([M_.dname filesep 'metropolis' filesep fname_ '_param_irf*.mat']);
lpmat=[]; lpmat=[];
for j=1:length(filetoload) for j=1:length(filetoload)
load([Model.dname filesep 'metropolis' filesep fname_ '_param_irf',int2str(j),'.mat']) load([M_.dname filesep 'metropolis' filesep fname_ '_param_irf',int2str(j),'.mat'])
lpmat = [lpmat; stock]; lpmat = [lpmat; stock];
clear stock clear stock
end end
type = 'post'; type = 'post';
else else
if DynareOptions.opt_gsa.pprior if options_.opt_gsa.pprior
filetoload=[OutputDirectoryName '/' fname_ '_prior']; filetoload=[OutputDirectoryName '/' fname_ '_prior'];
load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox') load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox')
lpmat = [lpmat0 lpmat]; lpmat = [lpmat0 lpmat];
@ -84,31 +85,31 @@ end
npar = size(pnames,1); npar = size(pnames,1);
nshock = np - npar; nshock = np - npar;
nbr_irf_restrictions = size(DynareOptions.endogenous_prior_restrictions.irf,1); nbr_irf_restrictions = size(options_.endogenous_prior_restrictions.irf,1);
nbr_moment_restrictions = size(DynareOptions.endogenous_prior_restrictions.moment,1); nbr_moment_restrictions = size(options_.endogenous_prior_restrictions.moment,1);
if init if init
mat_irf=cell(nbr_irf_restrictions,1); mat_irf=cell(nbr_irf_restrictions,1);
for ij=1:nbr_irf_restrictions for ij=1:nbr_irf_restrictions
mat_irf{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.irf{ij,3})); mat_irf{ij}=NaN(Nsam,length(options_.endogenous_prior_restrictions.irf{ij,3}));
end end
mat_moment=cell(nbr_moment_restrictions,1); mat_moment=cell(nbr_moment_restrictions,1);
for ij=1:nbr_moment_restrictions for ij=1:nbr_moment_restrictions
mat_moment{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.moment{ij,3})); mat_moment{ij}=NaN(Nsam,length(options_.endogenous_prior_restrictions.moment{ij,3}));
end end
irestrictions = [1:Nsam]; irestrictions = [1:Nsam];
h = dyn_waitbar(0,'Please wait...'); h = dyn_waitbar(0,'Please wait...');
for j=1:Nsam for j=1:Nsam
Model = set_all_parameters(lpmat(j,:)',EstimatedParameters,Model); M_ = set_all_parameters(lpmat(j,:)',estim_params_,M_);
if nbr_moment_restrictions if nbr_moment_restrictions
[Tt,Rr,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state); [Tt,Rr,SteadyState,info,oo_.dr, M_.params] = dynare_resolve(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
else else
[Tt,Rr,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict'); [Tt,Rr,SteadyState,info,oo_.dr, M_.params] = dynare_resolve(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state,'restrict');
end end
if info(1)==0 if info(1)==0
[info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults.dr,DynareResults.steady_state,DynareResults.exo_steady_state,DynareResults.exo_det_steady_state); [info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if ~isempty(info_irf) if ~isempty(info_irf)
for ij=1:nbr_irf_restrictions for ij=1:nbr_irf_restrictions
mat_irf{ij}(j,:)=data_irf{ij}(:,2)'; mat_irf{ij}(j,:)=data_irf{ij}(:,2)';
@ -131,7 +132,7 @@ if init
irestrictions=irestrictions(find(irestrictions)); irestrictions=irestrictions(find(irestrictions));
xmat=lpmat(irestrictions,:); xmat=lpmat(irestrictions,:);
skipline() skipline()
endo_prior_restrictions=DynareOptions.endogenous_prior_restrictions; endo_prior_restrictions=options_.endogenous_prior_restrictions;
save([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions'); save([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
else else
load([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions'); load([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
@ -190,8 +191,8 @@ if ~isempty(indx_irf)
iplot_indx = ones(size(plot_indx)); iplot_indx = ones(size(plot_indx));
indx_irf = indx_irf(irestrictions,:); indx_irf = indx_irf(irestrictions,:);
if ~DynareOptions.nograph if ~options_.nograph
h1=dyn_figure(DynareOptions.nodisplay,'name',[type ' evaluation of irf restrictions']); h1=dyn_figure(options_.nodisplay,'name',[type ' evaluation of irf restrictions']);
nrow=ceil(sqrt(nbr_irf_couples)); nrow=ceil(sqrt(nbr_irf_couples));
ncol=nrow; ncol=nrow;
if nrow*(nrow-1)>nbr_irf_couples if nrow*(nrow-1)>nbr_irf_couples
@ -204,7 +205,7 @@ if ~isempty(indx_irf)
indx_irf_matrix(:,plot_indx(ij)) = indx_irf_matrix(:,plot_indx(ij)) + indx_irf(:,ij); indx_irf_matrix(:,plot_indx(ij)) = indx_irf_matrix(:,plot_indx(ij)) + indx_irf(:,ij);
for ik=1:size(mat_irf{ij},2) for ik=1:size(mat_irf{ij},2)
[Mean,Median,Var,HPD,Distrib] = ... [Mean,Median,Var,HPD,Distrib] = ...
posterior_moments(mat_irf{ij}(:,ik),0,DynareOptions.mh_conf_sig); posterior_moments(mat_irf{ij}(:,ik),0,options_.mh_conf_sig);
irf_mean{plot_indx(ij)} = [irf_mean{plot_indx(ij)}; Mean]; irf_mean{plot_indx(ij)} = [irf_mean{plot_indx(ij)}; Mean];
irf_median{plot_indx(ij)} = [irf_median{plot_indx(ij)}; Median]; irf_median{plot_indx(ij)} = [irf_median{plot_indx(ij)}; Median];
irf_var{plot_indx(ij)} = [irf_var{plot_indx(ij)}; Var]; irf_var{plot_indx(ij)} = [irf_var{plot_indx(ij)}; Var];
@ -218,7 +219,7 @@ if ~isempty(indx_irf)
aleg = [aleg,'-' ,num2str(endo_prior_restrictions.irf{ij,3}(end))]; aleg = [aleg,'-' ,num2str(endo_prior_restrictions.irf{ij,3}(end))];
iplot_indx(ij)=0; iplot_indx(ij)=0;
end end
if ~DynareOptions.nograph && length(time_matrix{plot_indx(ij)})==1 if ~options_.nograph && length(time_matrix{plot_indx(ij)})==1
set(0,'currentfigure',h1), set(0,'currentfigure',h1),
subplot(nrow,ncol, plot_indx(ij)), subplot(nrow,ncol, plot_indx(ij)),
hc = cumplot(mat_irf{ij}(:,ik)); hc = cumplot(mat_irf{ij}(:,ik));
@ -258,7 +259,7 @@ if ~isempty(indx_irf)
options_mcf.nobeha_title = 'NO IRF restriction'; options_mcf.nobeha_title = 'NO IRF restriction';
options_mcf.title = atitle0; options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2) if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, DynareOptions); mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, options_);
end end
% [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0); % [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
@ -269,7 +270,7 @@ if ~isempty(indx_irf)
end end
for ij=1:nbr_irf_couples for ij=1:nbr_irf_couples
if length(time_matrix{ij})>1 if length(time_matrix{ij})>1
if ~DynareOptions.nograph if ~options_.nograph
set(0,'currentfigure',h1); set(0,'currentfigure',h1);
subplot(nrow,ncol, ij) subplot(nrow,ncol, ij)
itmp = (find(plot_indx==ij)); itmp = (find(plot_indx==ij));
@ -319,14 +320,14 @@ if ~isempty(indx_irf)
options_mcf.nobeha_title = 'NO IRF restriction'; options_mcf.nobeha_title = 'NO IRF restriction';
options_mcf.title = atitle0; options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2) if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, DynareOptions); mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, options_);
end end
end end
end end
end end
if ~DynareOptions.nograph if ~options_.nograph
dyn_saveas(h1,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],DynareOptions.nodisplay,DynareOptions.graph_format); dyn_saveas(h1,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],options_.nodisplay,options_.graph_format);
create_TeX_loader(DynareOptions,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],[type ' evaluation of irf restrictions'],'irf_restrictions',type,DynareOptions.figures.textwidth*min(ij/ncol,1)) create_TeX_loader(options_,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],[type ' evaluation of irf restrictions'],'irf_restrictions',type,options_.figures.textwidth*min(ij/ncol,1))
end end
skipline() skipline()
end end
@ -362,24 +363,24 @@ if ~isempty(indx_moment)
skipline() skipline()
%get parameter names including standard deviations %get parameter names including standard deviations
np=size(BayesInfo.name,1); np=size(bayestopt_.name,1);
name=cell(np,1); name=cell(np,1);
name_tex=cell(np,1); name_tex=cell(np,1);
for jj=1:np for jj=1:np
if DynareOptions.TeX if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(jj,DynareOptions.TeX,Model,EstimatedParameters,DynareOptions); [param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_);
name_tex{jj,1} = strrep(param_name_tex_temp,'$',''); name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
name{jj,1} = param_name_temp; name{jj,1} = param_name_temp;
else else
param_name_temp = get_the_name(jj,DynareOptions.TeX,Model,EstimatedParameters,DynareOptions); param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_);
name{jj,1} = param_name_temp; name{jj,1} = param_name_temp;
end end
end end
options_mcf.param_names = name; options_mcf.param_names = name;
if DynareOptions.TeX if options_.TeX
options_mcf.param_names_tex = name_tex; options_mcf.param_names_tex = name_tex;
end end
options_mcf.param_names = BayesInfo.name; options_mcf.param_names = bayestopt_.name;
all_moment_couples = cellstr([char(endo_prior_restrictions.moment(:,1)) char(endo_prior_restrictions.moment(:,2))]); all_moment_couples = cellstr([char(endo_prior_restrictions.moment(:,1)) char(endo_prior_restrictions.moment(:,2))]);
moment_couples = unique(all_moment_couples); moment_couples = unique(all_moment_couples);
nbr_moment_couples = size(moment_couples,1); nbr_moment_couples = size(moment_couples,1);
@ -405,8 +406,8 @@ if ~isempty(indx_moment)
iplot_indx = ones(size(plot_indx)); iplot_indx = ones(size(plot_indx));
indx_moment = indx_moment(irestrictions,:); indx_moment = indx_moment(irestrictions,:);
if ~DynareOptions.nograph if ~options_.nograph
h2=dyn_figure(DynareOptions.nodisplay,'name',[type ' evaluation of moment restrictions']); h2=dyn_figure(options_.nodisplay,'name',[type ' evaluation of moment restrictions']);
nrow=ceil(sqrt(nbr_moment_couples)); nrow=ceil(sqrt(nbr_moment_couples));
ncol=nrow; ncol=nrow;
if nrow*(nrow-1)>nbr_moment_couples if nrow*(nrow-1)>nbr_moment_couples
@ -420,7 +421,7 @@ if ~isempty(indx_moment)
indx_moment_matrix(:,plot_indx(ij)) = indx_moment_matrix(:,plot_indx(ij)) + indx_moment(:,ij); indx_moment_matrix(:,plot_indx(ij)) = indx_moment_matrix(:,plot_indx(ij)) + indx_moment(:,ij);
for ik=1:size(mat_moment{ij},2) for ik=1:size(mat_moment{ij},2)
[Mean,Median,Var,HPD,Distrib] = ... [Mean,Median,Var,HPD,Distrib] = ...
posterior_moments(mat_moment{ij}(:,ik),0,DynareOptions.mh_conf_sig); posterior_moments(mat_moment{ij}(:,ik),0,options_.mh_conf_sig);
moment_mean{plot_indx(ij)} = [moment_mean{plot_indx(ij)}; Mean]; moment_mean{plot_indx(ij)} = [moment_mean{plot_indx(ij)}; Mean];
moment_median{plot_indx(ij)} = [moment_median{plot_indx(ij)}; Median]; moment_median{plot_indx(ij)} = [moment_median{plot_indx(ij)}; Median];
moment_var{plot_indx(ij)} = [moment_var{plot_indx(ij)}; Var]; moment_var{plot_indx(ij)} = [moment_var{plot_indx(ij)}; Var];
@ -434,7 +435,7 @@ if ~isempty(indx_moment)
aleg = [aleg,'_' ,num2str(endo_prior_restrictions.moment{ij,3}(end))]; aleg = [aleg,'_' ,num2str(endo_prior_restrictions.moment{ij,3}(end))];
iplot_indx(ij)=0; iplot_indx(ij)=0;
end end
if ~DynareOptions.nograph && length(time_matrix{plot_indx(ij)})==1 if ~options_.nograph && length(time_matrix{plot_indx(ij)})==1
set(0,'currentfigure',h2); set(0,'currentfigure',h2);
subplot(nrow,ncol,plot_indx(ij)), subplot(nrow,ncol,plot_indx(ij)),
hc = cumplot(mat_moment{ij}(:,ik)); hc = cumplot(mat_moment{ij}(:,ik));
@ -469,7 +470,7 @@ if ~isempty(indx_moment)
options_mcf.nobeha_title = 'NO moment restriction'; options_mcf.nobeha_title = 'NO moment restriction';
options_mcf.title = atitle0; options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2) if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat, indx1, indx2, options_mcf, DynareOptions); mcf_analysis(xmat, indx1, indx2, options_mcf, options_);
end end
% [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0); % [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
@ -481,7 +482,7 @@ if ~isempty(indx_moment)
for ij=1:nbr_moment_couples for ij=1:nbr_moment_couples
time_matrix{ij} = sort(time_matrix{ij}); time_matrix{ij} = sort(time_matrix{ij});
if length(time_matrix{ij})>1 if length(time_matrix{ij})>1
if ~DynareOptions.nograph if ~options_.nograph
itmp = (find(plot_indx==ij)); itmp = (find(plot_indx==ij));
set(0,'currentfigure',h2); set(0,'currentfigure',h2);
subplot(nrow,ncol, ij) subplot(nrow,ncol, ij)
@ -531,14 +532,14 @@ if ~isempty(indx_moment)
options_mcf.nobeha_title = 'NO moment restriction'; options_mcf.nobeha_title = 'NO moment restriction';
options_mcf.title = atitle0; options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2) if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat, indx1, indx2, options_mcf, DynareOptions); mcf_analysis(xmat, indx1, indx2, options_mcf, options_);
end end
end end
end end
end end
if ~DynareOptions.nograph if ~options_.nograph
dyn_saveas(h2,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],DynareOptions.nodisplay,DynareOptions.graph_format); dyn_saveas(h2,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],options_.nodisplay,options_.graph_format);
create_TeX_loader(DynareOptions,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],[type ' evaluation of moment restrictions'],'moment_restrictions',type,DynareOptions.figures.textwidth*min(ij/ncol,1)) create_TeX_loader(options_,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],[type ' evaluation of moment restrictions'],'moment_restrictions',type,options_.figures.textwidth*min(ij/ncol,1))
end end
skipline() skipline()

15
matlab/gsa/mcf_analysis.m
View File

@ -1,12 +1,13 @@
function indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, DynareOptions) function indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, options_)
% % indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, options_)
% Written by Marco Ratto % Written by Marco Ratto
% Joint Research Centre, The European Commission, % Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu % marco.ratto@ec.europa.eu
% %
% Copyright © 2014 European Commission % Copyright © 2014 European Commission
% Copyright © 2016-2018 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -28,7 +29,7 @@ pvalue_corr = options_mcf.pvalue_corr;
alpha2 = options_mcf.alpha2; alpha2 = options_mcf.alpha2;
param_names = options_mcf.param_names; param_names = options_mcf.param_names;
if DynareOptions.TeX if options_.TeX
if ~isfield(options_mcf,'param_names_tex') if ~isfield(options_mcf,'param_names_tex')
param_names_tex = options_mcf.param_names; param_names_tex = options_mcf.param_names;
else else
@ -60,7 +61,7 @@ if ~isempty(indmcf)
data_mat=[dproba(indmcf)' proba(indmcf)']; data_mat=[dproba(indmcf)' proba(indmcf)'];
options_temp.noprint=0; options_temp.noprint=0;
dyntable(options_temp,['Smirnov statistics in driving ', title],headers,labels,data_mat,size(labels,2)+2,16,3); dyntable(options_temp,['Smirnov statistics in driving ', title],headers,labels,data_mat,size(labels,2)+2,16,3);
if DynareOptions.TeX if options_.TeX
labels_TeX=param_names_tex(indmcf); labels_TeX=param_names_tex(indmcf);
M_temp.dname=OutputDirectoryName ; M_temp.dname=OutputDirectoryName ;
M_temp.fname=fname_; M_temp.fname=fname_;
@ -76,11 +77,11 @@ if length(ibeha)>10 && length(inobeha)>10
indcorr = indcorr(~ismember(indcorr(:),indmcf)); indcorr = indcorr(~ismember(indcorr(:),indmcf));
indmcf = [indmcf(:); indcorr(:)]; indmcf = [indmcf(:); indcorr(:)];
end end
if ~isempty(indmcf) && ~DynareOptions.nograph if ~isempty(indmcf) && ~options_.nograph
skipline() skipline()
xx=[]; xx=[];
if ~ isempty(xparam1), xx=xparam1(indmcf); end if ~ isempty(xparam1), xx=xparam1(indmcf); end
scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ... scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
'.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, DynareOptions, ... '.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ...
beha_title, nobeha_title) beha_title, nobeha_title)
end end

8
matlab/gsa/scatter_analysis.m
View File

@ -1,4 +1,4 @@
function indmcf = scatter_analysis(lpmat, xdata, options_scatter, DynareOptions) function indmcf = scatter_analysis(lpmat, xdata, options_scatter, options_)
% %
% Written by Marco Ratto % Written by Marco Ratto
% Joint Research Centre, The European Commission, % Joint Research Centre, The European Commission,
@ -25,7 +25,7 @@ function indmcf = scatter_analysis(lpmat, xdata, options_scatter, DynareOptions)
param_names = options_scatter.param_names; param_names = options_scatter.param_names;
if DynareOptions.TeX if options_.TeX
if ~isfield(options_scatter,'param_names_tex') if ~isfield(options_scatter,'param_names_tex')
param_names_tex = options_scatter.param_names; param_names_tex = options_scatter.param_names;
else else
@ -42,11 +42,11 @@ if isfield(options_scatter,'xparam1')
end end
OutputDirectoryName = options_scatter.OutputDirectoryName; OutputDirectoryName = options_scatter.OutputDirectoryName;
if ~DynareOptions.nograph if ~options_.nograph
skipline() skipline()
xx=[]; xx=[];
if ~isempty(xparam1) if ~isempty(xparam1)
xx=xparam1; xx=xparam1;
end end
scatter_plots(lpmat, xdata, param_names, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, DynareOptions) scatter_plots(lpmat, xdata, param_names, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, options_)
end end

9
matlab/gsa/scatter_mcf.m
View File

@ -1,4 +1,5 @@
function scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, DynareOptions, beha_name, non_beha_name) function scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, options_, beha_name, non_beha_name)
% scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, options_, beha_name, non_beha_name)
% %
% Written by Marco Ratto % Written by Marco Ratto
% Joint Research Centre, The European Commission, % Joint Research Centre, The European Commission,
@ -84,7 +85,7 @@ figtitle_tex=strrep(figtitle,'_','\_');
fig_nam_=[fnam]; fig_nam_=[fnam];
if ~nograph if ~nograph
hh_fig=dyn_figure(DynareOptions.nodisplay,'name',figtitle); hh_fig=dyn_figure(options_.nodisplay,'name',figtitle);
end end
bf = 0.1; bf = 0.1;
@ -166,8 +167,8 @@ if ~isoctave
end end
if ~nograph if ~nograph
dyn_saveas(hh_fig,[dirname,filesep,fig_nam_],DynareOptions.nodisplay,DynareOptions.graph_format); dyn_saveas(hh_fig,[dirname,filesep,fig_nam_],options_.nodisplay,options_.graph_format);
if DynareOptions.TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format))) if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([dirname,'/',fig_nam_ '.tex'],'w'); fidTeX = fopen([dirname,'/',fig_nam_ '.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by scatter_mcf.m (Dynare).\n'); fprintf(fidTeX,'%% TeX eps-loader file generated by scatter_mcf.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']); fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);

10
matlab/gsa/scatter_plots.m
View File

@ -1,4 +1,4 @@
function scatter_plots(X,xp,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, DynareOptions) function scatter_plots(X,xp,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, options_)
% %
% Written by Marco Ratto % Written by Marco Ratto
% Joint Research Centre, The European Commission, % Joint Research Centre, The European Commission,
@ -58,7 +58,7 @@ end
if nargin<6 || isempty(dirname) if nargin<6 || isempty(dirname)
dirname=''; dirname='';
nograph=1; nograph=1;
DynareOptions.nodisplay=0; options_.nodisplay=0;
else else
nograph=0; nograph=0;
end end
@ -73,7 +73,7 @@ figtitle_tex=strrep(figtitle,'_','\_');
fig_nam_=[fnam]; fig_nam_=[fnam];
hh_fig=dyn_figure(DynareOptions.nodisplay,'name',figtitle); hh_fig=dyn_figure(options_.nodisplay,'name',figtitle);
set(hh_fig,'userdata',{X,xp}) set(hh_fig,'userdata',{X,xp})
bf = 0.1; bf = 0.1;
@ -172,8 +172,8 @@ end
% end % end
if ~nograph if ~nograph
dyn_saveas(hh_fig,[dirname,filesep,fig_nam_],DynareOptions.nodisplay,DynareOptions.graph_format); dyn_saveas(hh_fig,[dirname,filesep,fig_nam_],options_.nodisplay,options_.graph_format);
if DynareOptions.TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format))) if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([dirname,'/',fig_nam_ '.tex'],'w'); fidTeX = fopen([dirname,'/',fig_nam_ '.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by scatter_plots.m (Dynare).\n'); fprintf(fidTeX,'%% TeX eps-loader file generated by scatter_plots.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']); fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);

62
matlab/histvalf_initvalf.m
View File

@ -1,10 +1,10 @@
function [series, p] = histvalf_initvalf(caller, DynareModel, options) function [series, p] = histvalf_initvalf(caller, M_, options)
% [series, p] = histvalf_initvalf(caller, M_, options)
% Handles options for histval_file and initval_file % Handles options for histval_file and initval_file
% %
% INPUTS % INPUTS
% - caller [char] row array, name of calling function % - caller [char] row array, name of calling function
% - DynareModel [struct] model description, a.k.a M_ % - M_ [struct] model description, a.k.a
% - options [struct] options specific to histval_file and initval_file % - options [struct] options specific to histval_file and initval_file
% %
% OUTPUTS % OUTPUTS
@ -79,23 +79,23 @@ end
% checking that all variable are present % checking that all variable are present
error_flag = false; error_flag = false;
for i = 1:DynareModel.orig_endo_nbr for i = 1:M_.orig_endo_nbr
if ~series.exist(DynareModel.endo_names{i}) if ~series.exist(M_.endo_names{i})
dprintf('%s_FILE: endogenous variable %s is missing', caller, DynareModel.endo_names{i}) dprintf('%s_FILE: endogenous variable %s is missing', caller, M_.endo_names{i})
error_flag = true; error_flag = true;
end end
end end
for i = 1:DynareModel.exo_nbr for i = 1:M_.exo_nbr
if ~series.exist(DynareModel.exo_names{i}) if ~series.exist(M_.exo_names{i})
dprintf('%s_FILE: exogenous variable %s is missing', caller, DynareModel.exo_names{i}) dprintf('%s_FILE: exogenous variable %s is missing', caller, M_.exo_names{i})
error_flag = true; error_flag = true;
end end
end end
for i = 1:DynareModel.exo_det_nbr for i = 1:M_.exo_det_nbr
if ~series.exist(DynareModel.exo_det_names{i}) if ~series.exist(M_.exo_det_names{i})
dprintf('%s_FILE: exo_det variable %s is missing', caller, DynareModel.exo_det_names{i}) dprintf('%s_FILE: exo_det variable %s is missing', caller, M_.exo_det_names{i})
error_flag = true; error_flag = true;
end end
end end
@ -104,8 +104,8 @@ if error_flag
error('%s_FILE: some variables are missing', caller) error('%s_FILE: some variables are missing', caller)
end end
if exist(sprintf('+%s/dynamic_set_auxiliary_series.m', DynareModel.fname), 'file') if exist(sprintf('+%s/dynamic_set_auxiliary_series.m', M_.fname), 'file')
series = feval(sprintf('%s.dynamic_set_auxiliary_series', DynareModel.fname), series, DynareModel.params); series = feval(sprintf('%s.dynamic_set_auxiliary_series', M_.fname), series, M_.params);
end end
% selecting observations % selecting observations
@ -129,18 +129,18 @@ if isfield(options, 'first_obs') && ~isempty(options.first_obs)
error('%s_FILE: first_obs = %d is larger than the number of observations in the data file (%d)', ... error('%s_FILE: first_obs = %d is larger than the number of observations in the data file (%d)', ...
caller, options.first_obs, nobs0) caller, options.first_obs, nobs0)
elseif isfield(options, 'first_simulation_period') elseif isfield(options, 'first_simulation_period')
if options.first_obs == options.first_simulation_period - DynareModel.orig_maximum_lag if options.first_obs == options.first_simulation_period - M_.orig_maximum_lag
first_obs = periods(options.first_obs); first_obs = periods(options.first_obs);
else else
error('%s_FILE: first_obs = %d and first_simulation_period = %d have values inconsistent with a maximum lag of %d periods', ... error('%s_FILE: first_obs = %d and first_simulation_period = %d have values inconsistent with a maximum lag of %d periods', ...
caller, options.first_obs, options.first_simulation_period, DynareModel.orig_maximum_lag) caller, options.first_obs, options.first_simulation_period, M_.orig_maximum_lag)
end end
elseif isfield(options, 'firstsimulationperiod') elseif isfield(options, 'firstsimulationperiod')
if periods(options.first_obs) == options.firstsimulationperiod - DynareModel.orig_maximum_lag if periods(options.first_obs) == options.firstsimulationperiod - M_.orig_maximum_lag
first_obs = periods(options.first_obs); first_obs = periods(options.first_obs);
else else
error('%s_FILE: first_obs = %d and first_simulation_period = %s have values inconsistent with a maximum lag of %d periods', ... error('%s_FILE: first_obs = %d and first_simulation_period = %s have values inconsistent with a maximum lag of %d periods', ...
caller, options.first_obs, options.firstsimulationperiod, DynareModel.orig_maximum_lag) caller, options.first_obs, options.firstsimulationperiod, M_.orig_maximum_lag)
end end
else else
first_obs = periods(options.first_obs); first_obs = periods(options.first_obs);
@ -150,18 +150,18 @@ end
if isfield(options, 'firstobs') && ~isempty(options.firstobs) if isfield(options, 'firstobs') && ~isempty(options.firstobs)
if isfield(options, 'first_simulation_period') if isfield(options, 'first_simulation_period')
if options.firstobs == periods(options.first_simulation_period) - DynareModel.orig_maximum_lag if options.firstobs == periods(options.first_simulation_period) - M_.orig_maximum_lag
first_obs = options.firstobs; first_obs = options.firstobs;
else else
error('%s_FILE: first_obs = %s and first_simulation_period = %d have values inconsistent with a maximum lag of %d periods', ... error('%s_FILE: first_obs = %s and first_simulation_period = %d have values inconsistent with a maximum lag of %d periods', ...
caller, options.firstobs, options.first_simulation_period, DynareModel.orig_maximum_lag) caller, options.firstobs, options.first_simulation_period, M_.orig_maximum_lag)
end end
elseif isfield(options, 'firstsimulationperiod') elseif isfield(options, 'firstsimulationperiod')
if options.firstobs == options.firstsimulationperiod - DynareModel.orig_maximum_lag if options.firstobs == options.firstsimulationperiod - M_.orig_maximum_lag
first_obs = options.firstobs; first_obs = options.firstobs;
else else
error('%s_FILE: firstobs = %s and first_simulation_period = %s have values inconsistent with a maximum lag of %d periods', ... error('%s_FILE: firstobs = %s and first_simulation_period = %s have values inconsistent with a maximum lag of %d periods', ...
caller, options.firstobs, options.firstsimulationperiod, DynareModel.orig_maximum_lag) caller, options.firstobs, options.firstsimulationperiod, M_.orig_maximum_lag)
end end
else else
first_obs = options.firstobs; first_obs = options.firstobs;
@ -171,18 +171,18 @@ end
if ~first_obs_ispresent if ~first_obs_ispresent
if isfield(options, 'first_simulation_period') if isfield(options, 'first_simulation_period')
if options.first_simulation_period < DynareModel.orig_maximum_lag if options.first_simulation_period < M_.orig_maximum_lag
error('%s_FILE: first_simulation_period = %d must be larger than the maximum lag (%d)', ... error('%s_FILE: first_simulation_period = %d must be larger than the maximum lag (%d)', ...
caller, options.first_simulation_period, DynareModel.orig_maximum_lag) caller, options.first_simulation_period, M_.orig_maximum_lag)
elseif options.first_simulation_period > nobs0 elseif options.first_simulation_period > nobs0
error('%s_FILE: first_simulations_period = %d is larger than the number of observations in the data file (%d)', ... error('%s_FILE: first_simulations_period = %d is larger than the number of observations in the data file (%d)', ...
caller, options.first_obs, nobs0) caller, options.first_obs, nobs0)
else else
first_obs = periods(options.first_simulation_period) - DynareModel.orig_maximum_lag; first_obs = periods(options.first_simulation_period) - M_.orig_maximum_lag;
end end
first_obs_ispresent = true; first_obs_ispresent = true;
elseif isfield(options, 'firstsimulationperiod') elseif isfield(options, 'firstsimulationperiod')
first_obs = options.firstsimulationperiod - DynareModel.orig_maximum_lag; first_obs = options.firstsimulationperiod - M_.orig_maximum_lag;
first_obs_ispresent = true; first_obs_ispresent = true;
end end
end end
@ -237,8 +237,8 @@ end
if exist('lastsimulationperiod', 'var') if exist('lastsimulationperiod', 'var')
if lastsimulationperiod<=last_obs-DynareModel.orig_maximum_lead if lastsimulationperiod<=last_obs-M_.orig_maximum_lead
last_obs = lastsimulationperiod+DynareModel.orig_maximum_lead; last_obs = lastsimulationperiod+M_.orig_maximum_lead;
else else
error('%s_FILE: LAST_SIMULATION_PERIOD is too large compared to the available data.', caller) error('%s_FILE: LAST_SIMULATION_PERIOD is too large compared to the available data.', caller)
end end
@ -247,11 +247,11 @@ end
if exist('lastsimulationperiod', 'var') && exist('firstsimulationperiod', 'var') if exist('lastsimulationperiod', 'var') && exist('firstsimulationperiod', 'var')
p = lastsimulationperiod-firstsimulationperiod+1; p = lastsimulationperiod-firstsimulationperiod+1;
elseif exist('lastsimulationperiod', 'var') elseif exist('lastsimulationperiod', 'var')
p = lastsimulationperiod-(first_obs+DynareModel.orig_maximum_lag)+1; p = lastsimulationperiod-(first_obs+M_.orig_maximum_lag)+1;
elseif exist('firstsimulationperiod', 'var') elseif exist('firstsimulationperiod', 'var')
p = (last_obs-DynareModel.orig_maximum_lead)-firstsimulationperiod+1; p = (last_obs-M_.orig_maximum_lead)-firstsimulationperiod+1;
else else
p = (last_obs-DynareModel.orig_maximum_lead)-(first_obs+DynareModel.orig_maximum_lag)+1; p = (last_obs-M_.orig_maximum_lead)-(first_obs+M_.orig_maximum_lag)+1;
end end
series = series(first_obs:last_obs); series = series(first_obs:last_obs);

16
matlab/initial_condition_decomposition.m
View File

@ -99,10 +99,10 @@ if ~isfield(oo_,'initval_decomposition') || isequal(varlist,0)
with_epilogue = options_.initial_condition_decomp.with_epilogue; with_epilogue = options_.initial_condition_decomp.with_epilogue;
options_.selected_variables_only = 0; %make sure all variables are stored options_.selected_variables_only = 0; %make sure all variables are stored
options_.plot_priors=0; options_.plot_priors=0;
[oo,M,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_); [oo_local,M,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
% reduced form % reduced form
dr = oo.dr; dr = oo_local.dr;
% data reordering % data reordering
order_var = dr.order_var; order_var = dr.order_var;
@ -114,7 +114,7 @@ if ~isfield(oo_,'initval_decomposition') || isequal(varlist,0)
B = dr.ghu; B = dr.ghu;
% initialization % initialization
gend = length(oo.SmoothedShocks.(M_.exo_names{1})); %+options_.forecast; gend = length(oo_local.SmoothedShocks.(M_.exo_names{1})); %+options_.forecast;
z = zeros(endo_nbr,endo_nbr+2,gend); z = zeros(endo_nbr,endo_nbr+2,gend);
z(:,end,:) = Smoothed_Variables_deviation_from_mean; z(:,end,:) = Smoothed_Variables_deviation_from_mean;
@ -155,15 +155,15 @@ end
if ~isequal(varlist,0) if ~isequal(varlist,0)
% if ~options_.no_graph.shock_decomposition % if ~options_.no_graph.shock_decomposition
oo=oo_; oo_local=oo_;
oo.shock_decomposition = oo_.initval_decomposition; oo_local.shock_decomposition = oo_.initval_decomposition;
if ~isempty(init2shocks) if ~isempty(init2shocks)
init2shocks = M_.init2shocks.(init2shocks); init2shocks = M_.init2shocks.(init2shocks);
n=size(init2shocks,1); n=size(init2shocks,1);
for i=1:n for i=1:n
j=strmatch(init2shocks{i}{1},M_.endo_names,'exact'); j=strmatch(init2shocks{i}{1},M_.endo_names,'exact');
oo.shock_decomposition(:,end-1,:)=oo.shock_decomposition(:,j,:)+oo.shock_decomposition(:,end-1,:); oo_local.shock_decomposition(:,end-1,:)=oo_local.shock_decomposition(:,j,:)+oo_local.shock_decomposition(:,end-1,:);
oo.shock_decomposition(:,j,:)=0; oo_local.shock_decomposition(:,j,:)=0;
end end
end end
M_.exo_names = M_.endo_names; M_.exo_names = M_.endo_names;
@ -173,5 +173,5 @@ if ~isequal(varlist,0)
options_.plot_shock_decomp.use_shock_groups = ''; options_.plot_shock_decomp.use_shock_groups = '';
options_.plot_shock_decomp.init_cond_decomp = 1; % private flag to plotting utilities options_.plot_shock_decomp.init_cond_decomp = 1; % private flag to plotting utilities
plot_shock_decomposition(M_,oo,options_,varlist); plot_shock_decomposition(M_,oo_local,options_,varlist);
end end

22
matlab/kalman/get_Qvec_heteroskedastic_filter.m
View File

@ -1,14 +1,14 @@
function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,Model) function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,M_)
% function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,Model) % function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,M_)
% %
% INPUTS % INPUTS
% Q: baseline non-heteroskadastic covariance matrix of shocks % Q: baseline non-heteroskadastic covariance matrix of shocks
% smpl: scalar storing end of sample % smpl: scalar storing end of sample
% Model: structure storing the model information % M_: structure storing the model information
% Outputs: % Outputs:
% Qvec: [n_exo by n_exo by smpl] array of covariance matrices % Qvec: [n_exo by n_exo by smpl] array of covariance matrices
% Copyright © 2020-21 Dynare Team % Copyright © 2020-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -28,31 +28,31 @@ function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,Model)
isqdiag = all(all(abs(Q-diag(diag(Q)))<1e-14)); % ie, the covariance matrix is diagonal... isqdiag = all(all(abs(Q-diag(diag(Q)))<1e-14)); % ie, the covariance matrix is diagonal...
Qvec=repmat(Q,[1 1 smpl+1]); Qvec=repmat(Q,[1 1 smpl+1]);
for k=1:smpl for k=1:smpl
inx = ~isnan(Model.heteroskedastic_shocks.Qvalue(:,k)); inx = ~isnan(M_.heteroskedastic_shocks.Qvalue(:,k));
if any(inx) if any(inx)
if isqdiag if isqdiag
Qvec(inx,inx,k)=diag(Model.heteroskedastic_shocks.Qvalue(inx,k)); Qvec(inx,inx,k)=diag(M_.heteroskedastic_shocks.Qvalue(inx,k));
else else
inx = find(inx); inx = find(inx);
for s=1:length(inx) for s=1:length(inx)
if Q(inx(s),inx(s))>1.e-14 if Q(inx(s),inx(s))>1.e-14
tmpscale = sqrt(Model.heteroskedastic_shocks.Qvalue(inx(s),k)./Q(inx(s),inx(s))); tmpscale = sqrt(M_.heteroskedastic_shocks.Qvalue(inx(s),k)./Q(inx(s),inx(s)));
Qvec(inx(s),:,k) = Qvec(inx(s),:,k).*tmpscale; Qvec(inx(s),:,k) = Qvec(inx(s),:,k).*tmpscale;
Qvec(:,inx(s),k) = Qvec(:,inx(s),k).*tmpscale; Qvec(:,inx(s),k) = Qvec(:,inx(s),k).*tmpscale;
else else
Qvec(inx(s),inx(s),k)=Model.heteroskedastic_shocks.Qvalue(inx(s),k); Qvec(inx(s),inx(s),k)=M_.heteroskedastic_shocks.Qvalue(inx(s),k);
end end
end end
end end
end end
inx = ~isnan(Model.heteroskedastic_shocks.Qscale(:,k)); inx = ~isnan(M_.heteroskedastic_shocks.Qscale(:,k));
if any(inx) if any(inx)
if isqdiag if isqdiag
Qvec(inx,inx,k)=Qvec(inx,inx,k).*diag(Model.heteroskedastic_shocks.Qscale(inx,k)); Qvec(inx,inx,k)=Qvec(inx,inx,k).*diag(M_.heteroskedastic_shocks.Qscale(inx,k));
else else
inx = find(inx); inx = find(inx);
for s=1:length(inx) for s=1:length(inx)
tmpscale = sqrt(Model.heteroskedastic_shocks.Qscale(inx(s),k)); tmpscale = sqrt(M_.heteroskedastic_shocks.Qscale(inx(s),k));
Qvec(inx(s),:,k) = Qvec(inx(s),:,k).*tmpscale; Qvec(inx(s),:,k) = Qvec(inx(s),:,k).*tmpscale;
Qvec(:,inx(s),k) = Qvec(:,inx(s),k).*tmpscale; Qvec(:,inx(s),k) = Qvec(:,inx(s),k).*tmpscale;
end end

23
matlab/lmmcp/dyn_lmmcp.m
View File

@ -1,6 +1,7 @@
function [endo_simul,info] = dyn_lmmcp(M,options,oo) function [endo_simul,info] = dyn_lmmcp(M_,options,oo_)
% [endo_simul,info] = dyn_lmmcp(M_,options,oo_)
% Copyright © 2014 Dynare Team % Copyright © 2014-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -17,13 +18,13 @@ function [endo_simul,info] = dyn_lmmcp(M,options,oo)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[lb,ub,eq_index] = get_complementarity_conditions(M); [lb,ub,eq_index] = get_complementarity_conditions(M_);
lead_lag_incidence = M.lead_lag_incidence; lead_lag_incidence = M_.lead_lag_incidence;
ny = M.endo_nbr; ny = M_.endo_nbr;
max_lag = M.maximum_endo_lag; max_lag = M_.maximum_endo_lag;
nyp = nnz(lead_lag_incidence(1,:)) ; nyp = nnz(lead_lag_incidence(1,:)) ;
iyp = find(lead_lag_incidence(1,:)>0) ; iyp = find(lead_lag_incidence(1,:)>0) ;
@ -42,10 +43,10 @@ stop = 0 ;
iz = [1:ny+nyp+nyf]; iz = [1:ny+nyp+nyf];
periods = options.periods; periods = options.periods;
steady_state = oo.steady_state; steady_state = oo_.steady_state;
params = M.params; params = M_.params;
endo_simul = oo.endo_simul; endo_simul = oo_.endo_simul;
exo_simul = oo.exo_simul; exo_simul = oo_.exo_simul;
i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)'); i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)');
i_cols_A1 = find(lead_lag_incidence(2:3,:)'); i_cols_A1 = find(lead_lag_incidence(2:3,:)');
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)'); i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
@ -54,7 +55,7 @@ i_upd = ny+(1:periods*ny);
x = endo_simul(:); x = endo_simul(:);
model_dynamic = str2func([M.fname,'.dynamic']); model_dynamic = str2func([M_.fname,'.dynamic']);
z = x(find(lead_lag_incidence')); z = x(find(lead_lag_incidence'));
[res,A] = model_dynamic(z, exo_simul, params, steady_state,2); [res,A] = model_dynamic(z, exo_simul, params, steady_state,2);
nnzA = nnz(A); nnzA = nnz(A);

10
matlab/lnsrch1.m
View File

@ -56,18 +56,18 @@ if ~isfinite(summ)
end end
eq_number_string=[eq_number_string, num2str(j1(end))]; eq_number_string=[eq_number_string, num2str(j1(end))];
var_string=[]; var_string=[];
Model=evalin('base','M_'); M_=evalin('base','M_');
for ii=1:length(j2)-1 for ii=1:length(j2)-1
var_string=[var_string, Model.endo_names{j2(ii)}, ', ']; var_string=[var_string, M_.endo_names{j2(ii)}, ', '];
end end
var_string=[var_string, Model.endo_names{j2(end)}]; var_string=[var_string, M_.endo_names{j2(end)}];
fprintf('\nAn infinite element was encountered when trying to solve equation(s) %s \n',eq_number_string) fprintf('\nAn infinite element was encountered when trying to solve equation(s) %s \n',eq_number_string)
fprintf('with respect to the variable(s): %s.\n',var_string) fprintf('with respect to the variable(s): %s.\n',var_string)
fprintf('The values of the endogenous variables when the problem was encountered were:\n') fprintf('The values of the endogenous variables when the problem was encountered were:\n')
label_width=size(char(Model.endo_names),2)+2; label_width=size(char(M_.endo_names),2)+2;
label_string=sprintf('%%-%us %%8.4g \\n',label_width); label_string=sprintf('%%-%us %%8.4g \\n',label_width);
for ii=1:length(xold) for ii=1:length(xold)
fprintf(label_string, Model.endo_names{ii}, xold(ii)); fprintf(label_string, M_.endo_names{ii}, xold(ii));
end end
skipline(); skipline();
end end

8
matlab/mex/k_order_perturbation.m
View File

@ -1,10 +1,10 @@
% [dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,DynareModel,DynareOptions) % [dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,M_,options_)
% computes a k-th order perturbation solution % computes a k-th order perturbation solution
% %
% INPUTS % INPUTS
% dr: struct describing the reduced form solution of the model. % dr: struct describing the reduced form solution of the model.
% DynareModel: struct jobs's parameters % M_: struct jobs's parameters
% DynareOptions: struct job's options % options_: struct job's options
% %
% OUTPUTS % OUTPUTS
% dynpp_derivs struct Derivatives of the decision rule in Dynare++ format. % dynpp_derivs struct Derivatives of the decision rule in Dynare++ format.
@ -25,7 +25,7 @@
% dynare/mex/sources/k_order_perturbation.cc and it uses code provided by % dynare/mex/sources/k_order_perturbation.cc and it uses code provided by
% dynare++ % dynare++
% Copyright © 2013-2021 Dynare Team % Copyright © 2013-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %

14
matlab/model_comparison.m
View File

@ -1,17 +1,17 @@
function oo = model_comparison(ModelNames,ModelPriors,oo,options_,fname) function oo_ = model_comparison(ModelNames,ModelPriors,oo_,options_,fname)
% function oo = model_comparison(ModelNames,ModelPriors,oo,options_,fname) % function oo_ = model_comparison(ModelNames,ModelPriors,oo_,options_,fname)
% Conducts Bayesian model comparison. This function computes Odds ratios and % Conducts Bayesian model comparison. This function computes Odds ratios and
% estimates a posterior density over a collection of models. % estimates a posterior density over a collection of models.
% %
% INPUTS % INPUTS
% ModelNames [string] m*1 cell array of string. % ModelNames [string] m*1 cell array of string.
% ModelPriors [double] m*1 vector of prior probabilities % ModelPriors [double] m*1 vector of prior probabilities
% oo [struct] Dynare results structure % oo_ [struct] Dynare results structure
% options_ [struct] Dynare options structure % options_ [struct] Dynare options structure
% fname [string] name of the current mod-file % fname [string] name of the current mod-file
% %
% OUTPUTS % OUTPUTS
% oo [struct] Dynare results structure containing the % oo_ [struct] Dynare results structure containing the
% results in a field PosteriorOddsTable % results in a field PosteriorOddsTable
% %
% ALGORITHM % ALGORITHM
@ -20,7 +20,7 @@ function oo = model_comparison(ModelNames,ModelPriors,oo,options_,fname)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2007-2018 Dynare Team % Copyright © 2007-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -78,7 +78,7 @@ iname = strmatch(fname,ShortModelNames,'exact');
for i=1:NumberOfModels for i=1:NumberOfModels
if i==iname if i==iname
mstruct.oo_ = oo; mstruct.oo_ = oo_;
else else
if length(ModelNames{i})>3 && (strcmpi(ModelNames{i}(end-3:end),'.mod') || strcmpi(ModelNames{i}(end-3:end),'.dyn')) if length(ModelNames{i})>3 && (strcmpi(ModelNames{i}(end-3:end),'.mod') || strcmpi(ModelNames{i}(end-3:end),'.dyn'))
mstruct = load([ModelNames{i}(1:end-4) filesep 'Output' ModelNames{i}(1:end-4) '_results.mat' ],'oo_'); mstruct = load([ModelNames{i}(1:end-4) filesep 'Output' ModelNames{i}(1:end-4) '_results.mat' ],'oo_');
@ -127,7 +127,7 @@ end
for model_iter = 1:NumberOfModels for model_iter = 1:NumberOfModels
for var_iter = 1:length(labels) for var_iter = 1:length(labels)
oo.Model_Comparison.(headers{1+model_iter}).(field_labels{var_iter}) = values(var_iter, model_iter); oo_.Model_Comparison.(headers{1+model_iter}).(field_labels{var_iter}) = values(var_iter, model_iter);
end end
end end

4
matlab/nonlinear-filters/auxiliary_initialization.m
View File

@ -46,7 +46,7 @@ if isempty(init_flag)
init_flag = 1; init_flag = 1;
end end
order = DynareOptions.order; order = options_.order;
% Set local state space model (first order approximation). % Set local state space model (first order approximation).
ghx = ReducedForm.ghx; ghx = ReducedForm.ghx;
@ -83,7 +83,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
%end %end
% Set seed for randn(). % Set seed for randn().
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
% Initialization of the likelihood. % Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables; const_lik = log(2*pi)*number_of_observed_variables;

14
matlab/nonlinear-filters/auxiliary_particle_filter.m
View File

@ -1,9 +1,9 @@
function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model) function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, options_, M_)
% [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, options_, M_)
% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter % Evaluates the likelihood of a nonlinear model with the auxiliary particle filter
% allowing eventually resampling. % allowing eventually resampling.
% %
% Copyright © 2011-2022 Dynare Team % Copyright © 2011-2023 Dynare Team
% %
% This file is part of Dynare (particles module). % This file is part of Dynare (particles module).
% %
@ -25,7 +25,7 @@ if isempty(start)
start = 1; start = 1;
end end
% Get perturbation order % Get perturbation order
order = DynareOptions.order; order = options_.order;
% Set flag for prunning % Set flag for prunning
pruning = ParticleOptions.pruning; pruning = ParticleOptions.pruning;
@ -77,7 +77,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = chol(Q)'; Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn(). % Set seed for randn().
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
% Initialization of the likelihood. % Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables+log(det(H)); const_lik = log(2*pi)*number_of_observed_variables+log(det(H));
@ -119,7 +119,7 @@ for t=1:sample_size
end end
else else
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations,number_of_particles), dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations,number_of_particles), dr, M_, options_, udr);
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
@ -152,7 +152,7 @@ for t=1:sample_size
StateVectors_ = tmp_(mf0_,:); StateVectors_ = tmp_(mf0_,:);
else else
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);

12
matlab/nonlinear-filters/conditional_filter_proposal.m
View File

@ -1,5 +1,5 @@
function [ProposalStateVector, Weights, flag] = conditional_filter_proposal(ReducedForm, y, StateVectors, SampleWeights, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ... function [ProposalStateVector, Weights, flag] = conditional_filter_proposal(ReducedForm, y, StateVectors, SampleWeights, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ...
H, ParticleOptions, ThreadsOptions, DynareOptions, Model) H, ParticleOptions, ThreadsOptions, options_, M_)
% Computes the proposal for each past particle using Gaussian approximations % Computes the proposal for each past particle using Gaussian approximations
% for the state errors and the Kalman filter % for the state errors and the Kalman filter
@ -14,8 +14,8 @@ function [ProposalStateVector, Weights, flag] = conditional_filter_proposal(Redu
% - H % - H
% - ParticleOptions % - ParticleOptions
% - ThreadsOptions % - ThreadsOptions
% - DynareOptions % - options_
% - Model % - M_
% %
% OUTPUTS % OUTPUTS
% - ProposalStateVector % - ProposalStateVector
@ -41,7 +41,7 @@ function [ProposalStateVector, Weights, flag] = conditional_filter_proposal(Redu
flag = false; flag = false;
order = DynareOptions.order; order = options_.order;
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
dr = ReducedForm.dr; dr = ReducedForm.dr;
@ -93,7 +93,7 @@ epsilon = Q_lower_triangular_cholesky*nodes';
yhat = repmat(StateVectors-state_variables_steady_state, 1, size(epsilon, 2)); yhat = repmat(StateVectors-state_variables_steady_state, 1, size(epsilon, 2));
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
@ -159,6 +159,6 @@ if ParticleOptions.cpf_weights_method.murrayjonesparslow
end end
Prior = probability2(PredictedStateMean, PredictedStateVarianceSquareRoot, ProposalStateVector); Prior = probability2(PredictedStateMean, PredictedStateVarianceSquareRoot, ProposalStateVector);
Posterior = probability2(StateVectorMean, StateVectorVarianceSquareRoot, ProposalStateVector); Posterior = probability2(StateVectorMean, StateVectorVarianceSquareRoot, ProposalStateVector);
Likelihood = probability2(y, H_lower_triangular_cholesky, measurement_equations(ProposalStateVector, ReducedForm, ThreadsOptions, DynareOptions, Model)); Likelihood = probability2(y, H_lower_triangular_cholesky, measurement_equations(ProposalStateVector, ReducedForm, ThreadsOptions, options_, M_));
Weights = SampleWeights.*Likelihood.*(Prior./Posterior); Weights = SampleWeights.*Likelihood.*(Prior./Posterior);
end end

10
matlab/nonlinear-filters/conditional_particle_filter.m
View File

@ -1,4 +1,4 @@
function [LIK,lik] = conditional_particle_filter(ReducedForm, Y, s, ParticleOptions, ThreadsOptions, DynareOptions, Model) function [LIK,lik] = conditional_particle_filter(ReducedForm, Y, s, ParticleOptions, ThreadsOptions, options_, M_)
% Evaluates the likelihood of a non-linear model with a particle filter % Evaluates the likelihood of a non-linear model with a particle filter
% %
@ -8,8 +8,8 @@ function [LIK,lik] = conditional_particle_filter(ReducedForm, Y, s, ParticleOpti
% - s [integer] scalar, likelihood evaluation starts at s (has to be smaller than T, the sample length provided in Y). % - s [integer] scalar, likelihood evaluation starts at s (has to be smaller than T, the sample length provided in Y).
% - ParticlesOptions [struct] % - ParticlesOptions [struct]
% - ThreadsOptions [struct] % - ThreadsOptions [struct]
% - DynareOptions [struct] % - options_ [struct]
% - Model [struct] % - M_ [struct]
% %
% OUTPUTS % OUTPUTS
% - LIK [double] scalar, likelihood % - LIK [double] scalar, likelihood
@ -78,7 +78,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot, 2);
Q_lower_triangular_cholesky = chol(Q)'; Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn(). % Set seed for randn().
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
% Initialization of the likelihood. % Initialization of the likelihood.
lik = NaN(T, 1); lik = NaN(T, 1);
@ -90,7 +90,7 @@ for t=1:T
flags = false(n, 1); flags = false(n, 1);
for i=1:n for i=1:n
[StateParticles(:,i), SampleWeights(i), flags(i)] = ... [StateParticles(:,i), SampleWeights(i), flags(i)] = ...
conditional_filter_proposal(ReducedForm, Y(:,t), StateParticles(:,i), SampleWeights(i), Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ParticleOptions, ThreadsOptions, DynareOptions, Model); conditional_filter_proposal(ReducedForm, Y(:,t), StateParticles(:,i), SampleWeights(i), Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ParticleOptions, ThreadsOptions, options_, M_);
end end
if any(flags) if any(flags)
LIK = -Inf; LIK = -Inf;

8
matlab/nonlinear-filters/gaussian_densities.m
View File

@ -1,5 +1,5 @@
function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions,DynareOptions, Model) function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions,options_, M_)
% % IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions,options_, M_)
% Elements to calculate the importance sampling ratio % Elements to calculate the importance sampling ratio
% %
% INPUTS % INPUTS
@ -20,7 +20,7 @@ function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sq
% NOTES % NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood. % The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright © 2009-2019 Dynare Team % Copyright © 2009-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -44,7 +44,7 @@ proposal = probability2(mut_t, sqr_Pss_t_t, particles);
prior = probability2(st_t_1, sqr_Pss_t_t_1, particles); prior = probability2(st_t_1, sqr_Pss_t_t_1, particles);
% likelihood % likelihood
yt_t_1_i = measurement_equations(particles, ReducedForm, ThreadsOptions, DynareOptions, Model); yt_t_1_i = measurement_equations(particles, ReducedForm, ThreadsOptions, options_, M_);
likelihood = probability2(obs, sqrt(H), yt_t_1_i); likelihood = probability2(obs, sqrt(H), yt_t_1_i);
IncrementalWeights = likelihood.*prior./proposal; IncrementalWeights = likelihood.*prior./proposal;

12
matlab/nonlinear-filters/gaussian_filter.m
View File

@ -1,5 +1,5 @@
function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model) function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, options_, M_)
% [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, options_, M_)
% Evaluates the likelihood of a non-linear model approximating the % Evaluates the likelihood of a non-linear model approximating the
% predictive (prior) and filtered (posterior) densities for state variables % predictive (prior) and filtered (posterior) densities for state variables
% by gaussian distributions. % by gaussian distributions.
@ -74,7 +74,7 @@ else
end end
if ParticleOptions.distribution_approximation.montecarlo if ParticleOptions.distribution_approximation.montecarlo
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
end end
% Get covariance matrices % Get covariance matrices
@ -101,20 +101,20 @@ LIK = NaN;
for t=1:sample_size for t=1:sample_size
[PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ... [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm, Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ... gaussian_filter_bank(ReducedForm, Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ...
H, ParticleOptions, ThreadsOptions, DynareOptions, Model); H, ParticleOptions, ThreadsOptions, options_, M_);
if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented
StateParticles = bsxfun(@plus, StateVectorMean, StateVectorVarianceSquareRoot*nodes2'); StateParticles = bsxfun(@plus, StateVectorMean, StateVectorVarianceSquareRoot*nodes2');
IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ... IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ...
PredictedStateVarianceSquareRoot, StateParticles, H, const_lik, ... PredictedStateVarianceSquareRoot, StateParticles, H, const_lik, ...
weights2, weights_c2, ReducedForm, ThreadsOptions, ... weights2, weights_c2, ReducedForm, ThreadsOptions, ...
DynareOptions, Model); options_, M_);
SampleWeights = weights2.*IncrementalWeights; SampleWeights = weights2.*IncrementalWeights;
else else
StateParticles = bsxfun(@plus, StateVectorVarianceSquareRoot*randn(state_variance_rank, number_of_particles), StateVectorMean) ; StateParticles = bsxfun(@plus, StateVectorVarianceSquareRoot*randn(state_variance_rank, number_of_particles), StateVectorMean) ;
IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ... IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik, ... PredictedStateVarianceSquareRoot,StateParticles,H,const_lik, ...
1/number_of_particles,1/number_of_particles,ReducedForm,ThreadsOptions, ... 1/number_of_particles,1/number_of_particles,ReducedForm,ThreadsOptions, ...
DynareOptions, Model); options_, M_);
SampleWeights = IncrementalWeights/number_of_particles; SampleWeights = IncrementalWeights/number_of_particles;
end end
SampleWeights = SampleWeights + 1e-6*ones(size(SampleWeights, 1), 1); SampleWeights = SampleWeights + 1e-6*ones(size(SampleWeights, 1), 1);

6
matlab/nonlinear-filters/gaussian_filter_bank.m
View File

@ -1,6 +1,6 @@
function [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ... function [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm, obs, StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ... gaussian_filter_bank(ReducedForm, obs, StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ...
ParticleOptions, ThreadsOptions, DynareOptions, Model) ParticleOptions, ThreadsOptions, options_, M_)
% %
% Computes the proposal with a gaussian approximation for importance % Computes the proposal with a gaussian approximation for importance
% sampling % sampling
@ -39,7 +39,7 @@ function [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean,
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
order = DynareOptions.order; order = options_.order;
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
dr = ReducedForm.dr; dr = ReducedForm.dr;
@ -95,7 +95,7 @@ StateVectors = sigma_points(1:number_of_state_variables,:);
epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:); epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
yhat = bsxfun(@minus, StateVectors, state_variables_steady_state); yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);

6
matlab/nonlinear-filters/gaussian_mixture_densities.m
View File

@ -1,6 +1,6 @@
function IncrementalWeights = gaussian_mixture_densities(obs, StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ... function IncrementalWeights = gaussian_mixture_densities(obs, StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ... StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
ReducedForm, ThreadsOptions, DynareOptions, Model) ReducedForm, ThreadsOptions, options_, M_)
% Elements to calculate the importance sampling ratio % Elements to calculate the importance sampling ratio
% %
@ -22,7 +22,7 @@ function IncrementalWeights = gaussian_mixture_densities(obs, StateMuPrior, Sta
% NOTES % NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood. % The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright © 2009-2019 Dynare Team % Copyright © 2009-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -48,7 +48,7 @@ prior = prior';
proposal = proposal'; proposal = proposal';
% Compute the density of the current observation conditionally to each particle % Compute the density of the current observation conditionally to each particle
yt_t_1_i = measurement_equations(StateParticles, ReducedForm, ThreadsOptions, DynareOptions, Model); yt_t_1_i = measurement_equations(StateParticles, ReducedForm, ThreadsOptions, options_, M_);
% likelihood % likelihood
likelihood = probability2(obs, sqrt(H), yt_t_1_i); likelihood = probability2(obs, sqrt(H), yt_t_1_i);

10
matlab/nonlinear-filters/gaussian_mixture_filter.m
View File

@ -1,4 +1,4 @@
function [LIK, lik] = gaussian_mixture_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model) function [LIK, lik] = gaussian_mixture_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, options_, M_)
% Evaluates the likelihood of a non-linear model approximating the state % Evaluates the likelihood of a non-linear model approximating the state
% variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing % variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing
@ -79,7 +79,7 @@ else
end end
if ParticleOptions.distribution_approximation.montecarlo if ParticleOptions.distribution_approximation.montecarlo
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
end end
% Get covariance matrices % Get covariance matrices
@ -178,7 +178,7 @@ for t=1:sample_size
gaussian_mixture_filter_bank(ReducedForm,Y(:,t), StateMu(:,g), StateSqrtP(:,:,g), StateWeights(g),... gaussian_mixture_filter_bank(ReducedForm,Y(:,t), StateMu(:,g), StateSqrtP(:,:,g), StateWeights(g),...
StructuralShocksMu(:,i), StructuralShocksSqrtP(:,:,i), StructuralShocksWeights(i),... StructuralShocksMu(:,i), StructuralShocksSqrtP(:,:,i), StructuralShocksWeights(i),...
ObservationShocksWeights(j), H, H_lower_triangular_cholesky, const_lik, ... ObservationShocksWeights(j), H, H_lower_triangular_cholesky, const_lik, ...
ParticleOptions, ThreadsOptions, DynareOptions, Model); ParticleOptions, ThreadsOptions, options_, M_);
end end
end end
end end
@ -192,7 +192,7 @@ for t=1:sample_size
StateParticles = bsxfun(@plus, StateMuPost(:,i), StateSqrtPPost(:,:,i)*nodes'); StateParticles = bsxfun(@plus, StateMuPost(:,i), StateSqrtPPost(:,:,i)*nodes');
IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ... IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ... StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
ReducedForm, ThreadsOptions, DynareOptions, Model); ReducedForm, ThreadsOptions, options_, M_);
SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights); SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights);
end end
SumSampleWeights = sum(SampleWeights); SumSampleWeights = sum(SampleWeights);
@ -210,7 +210,7 @@ for t=1:sample_size
StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles); StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles);
IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ... IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ... StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
ReducedForm, ThreadsOptions, DynareOptions, Model); ReducedForm, ThreadsOptions, options_, M_);
SampleWeights = IncrementalWeights/number_of_particles; SampleWeights = IncrementalWeights/number_of_particles;
SumSampleWeights = sum(SampleWeights,1); SumSampleWeights = sum(SampleWeights,1);
SampleWeights = SampleWeights./SumSampleWeights; SampleWeights = SampleWeights./SumSampleWeights;

6
matlab/nonlinear-filters/gaussian_mixture_filter_bank.m
View File

@ -2,7 +2,7 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPP
gaussian_mixture_filter_bank(ReducedForm, obs, StateMu, StateSqrtP, StateWeights, ... gaussian_mixture_filter_bank(ReducedForm, obs, StateMu, StateSqrtP, StateWeights, ...
StructuralShocksMu, StructuralShocksSqrtP, StructuralShocksWeights, ... StructuralShocksMu, StructuralShocksSqrtP, StructuralShocksWeights, ...
ObservationShocksWeights, H, H_lower_triangular_cholesky, normfactO, ... ObservationShocksWeights, H, H_lower_triangular_cholesky, normfactO, ...
ParticleOptions, ThreadsOptions, DynareOptions, Model) ParticleOptions, ThreadsOptions, options_, M_)
% Computes the proposal with a gaussian approximation for importance % Computes the proposal with a gaussian approximation for importance
% sampling % sampling
@ -41,7 +41,7 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPP
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
order = DynareOptions.order; order = options_.order;
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
dr = ReducedForm.dr; dr = ReducedForm.dr;
@ -91,7 +91,7 @@ epsilon = bsxfun(@plus, StructuralShocksSqrtP*nodes3(:,number_of_state_variables
StateVectors = bsxfun(@plus, StateSqrtP*nodes3(:,1:number_of_state_variables)', StateMu); StateVectors = bsxfun(@plus, StateSqrtP*nodes3(:,1:number_of_state_variables)', StateMu);
yhat = bsxfun(@minus, StateVectors, state_variables_steady_state); yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);

6
matlab/nonlinear-filters/measurement_equations.m
View File

@ -1,4 +1,4 @@
function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions, DynareOptions, Model) function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions, options_, M_)
% Copyright © 2013-2022 Dynare Team % Copyright © 2013-2022 Dynare Team
% %
@ -17,7 +17,7 @@ function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
order = DynareOptions.order; order = options_.order;
mf1 = ReducedForm.mf1; mf1 = ReducedForm.mf1;
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
dr = ReducedForm.dr; dr = ReducedForm.dr;
@ -44,7 +44,7 @@ state_variables_steady_state = ReducedForm.state_variables_steady_state;
number_of_structural_innovations = length(ReducedForm.Q); number_of_structural_innovations = length(ReducedForm.Q);
yhat = bsxfun(@minus, StateVectors, state_variables_steady_state); yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, M_, options_, udr);
measure = tmp(mf1,:); measure = tmp(mf1,:);
else else
if order == 2 if order == 2

8
matlab/nonlinear-filters/nonlinear_kalman_filter.m
View File

@ -1,4 +1,4 @@
function [LIK,lik] = nonlinear_kalman_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model) function [LIK,lik] = nonlinear_kalman_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, options_, M_)
% Evaluates the likelihood of a non-linear model approximating the % Evaluates the likelihood of a non-linear model approximating the
% predictive (prior) and filtered (posterior) densities for state variables % predictive (prior) and filtered (posterior) densities for state variables
@ -54,7 +54,7 @@ if isempty(start)
start = 1; start = 1;
end end
order = DynareOptions.order; order = options_.order;
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
dr = ReducedForm.dr; dr = ReducedForm.dr;
@ -105,7 +105,7 @@ else
end end
if ParticleOptions.distribution_approximation.montecarlo if ParticleOptions.distribution_approximation.montecarlo
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
end end
% Get covariance matrices % Get covariance matrices
@ -130,7 +130,7 @@ for t=1:sample_size
epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:); epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);

98
matlab/nonlinear-filters/online_auxiliary_filter.m
View File

@ -1,16 +1,16 @@
function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, DynareResults) function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
% [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters. % Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
% %
% INPUTS % INPUTS
% - xparam1 [double] n×1 vector, Initial condition for the estimated parameters. % - xparam1 [double] n×1 vector, Initial condition for the estimated parameters.
% - DynareDataset [dseries] Sample used for estimation. % - dataset_ [dseries] Sample used for estimation.
% - dataset_info [struct] Description of the sample. % - dataset_info [struct] Description of the sample.
% - DynareOptions [struct] Option values (options_). % - options_ [struct] Option values.
% - Model [struct] Description of the model (M_). % - M_ [struct] Description of the model.
% - EstimatedParameters [struct] Description of the estimated parameters (estim_params_). % - estim_params_ [struct] Description of the estimated parameters.
% - BayesInfo [struct] Prior definition (bayestopt_). % - bayestopt_ [struct] Prior definition.
% - DynareResults [struct] Results (oo_). % - oo_ [struct] Results.
% %
% OUTPUTS % OUTPUTS
% - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters). % - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters).
@ -39,26 +39,26 @@ function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxili
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Set seed for randn(). % Set seed for randn().
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
pruning = DynareOptions.particle.pruning; pruning = options_.particle.pruning;
second_resample = DynareOptions.particle.resampling.status.systematic; second_resample = options_.particle.resampling.status.systematic;
variance_update = true; variance_update = true;
bounds = prior_bounds(BayesInfo, DynareOptions.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding bounds = prior_bounds(bayestopt_, options_.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
% initialization of state particles % initialization of state particles
[~, Model, DynareOptions, DynareResults, ReducedForm] = solve_model_for_online_filter(true, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults); [~, M_, options_, oo_, ReducedForm] = solve_model_for_online_filter(true, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
order = DynareOptions.order; order = options_.order;
mf0 = ReducedForm.mf0; mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1; mf1 = ReducedForm.mf1;
number_of_particles = DynareOptions.particle.number_of_particles; number_of_particles = options_.particle.number_of_particles;
number_of_parameters = size(xparam1,1); number_of_parameters = size(xparam1,1);
Y = DynareDataset.data; Y = dataset_.data;
sample_size = size(Y,1); sample_size = size(Y,1);
number_of_observed_variables = length(mf1); number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q); number_of_structural_innovations = length(ReducedForm.Q);
liu_west_delta = DynareOptions.particle.liu_west_delta; liu_west_delta = options_.particle.liu_west_delta;
% Get initial conditions for the state particles % Get initial conditions for the state particles
StateVectorMean = ReducedForm.StateVectorMean; StateVectorMean = ReducedForm.StateVectorMean;
@ -81,12 +81,12 @@ b_square = 1-small_a*small_a;
% Initialization of parameter particles % Initialization of parameter particles
xparam = zeros(number_of_parameters,number_of_particles); xparam = zeros(number_of_parameters,number_of_particles);
Prior = dprior(BayesInfo, DynareOptions.prior_trunc); Prior = dprior(bayestopt_, options_.prior_trunc);
for i=1:number_of_particles for i=1:number_of_particles
info = 12042009; info = 12042009;
while info while info
candidate = Prior.draw(); candidate = Prior.draw();
[info, Model, DynareOptions, DynareResults] = solve_model_for_online_filter(false, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults); [info, M_, options_, oo_] = solve_model_for_online_filter(false, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
if ~info if ~info
xparam(:,i) = candidate(:); xparam(:,i) = candidate(:);
end end
@ -124,8 +124,8 @@ for t=1:sample_size
tau_tilde = zeros(1,number_of_particles); tau_tilde = zeros(1,number_of_particles);
for i=1:number_of_particles for i=1:number_of_particles
% model resolution % model resolution
[info, Model, DynareOptions, DynareResults, ReducedForm] = ... [info, M_, options_, oo_, ReducedForm] = ...
solve_model_for_online_filter(false, fore_xparam(:,i), DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults); solve_model_for_online_filter(false, fore_xparam(:,i), dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
if ~info(1) if ~info(1)
steadystate = ReducedForm.steadystate; steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state; state_variables_steady_state = ReducedForm.state_variables_steady_state;
@ -167,22 +167,22 @@ for t=1:sample_size
% particle likelihood contribution % particle likelihood contribution
yhat = bsxfun(@minus, StateVectors(:,i), state_variables_steady_state); yhat = bsxfun(@minus, StateVectors(:,i), state_variables_steady_state);
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, M_, options_, udr);
else else
if pruning if pruning
yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state_); yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state_);
if order == 2 if order == 2
[tmp, ~] = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2); [tmp, ~] = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, options_.threads.local_state_space_iteration_2);
elseif order == 3 elseif order == 3
[tmp, tmp_] = local_state_space_iteration_3(yhat_, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning); [tmp, tmp_] = local_state_space_iteration_3(yhat_, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, options_.threads.local_state_space_iteration_3, pruning);
else else
error('Pruning is not available for orders > 3'); error('Pruning is not available for orders > 3');
end end
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, options_.threads.local_state_space_iteration_2);
elseif order == 3 elseif order == 3
tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning); tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, options_.threads.local_state_space_iteration_3, pruning);
else else
error('Order > 3: use_k_order_solver should be set to true'); error('Order > 3: use_k_order_solver should be set to true');
end end
@ -196,7 +196,7 @@ for t=1:sample_size
end end
% particles selection % particles selection
tau_tilde = tau_tilde/sum(tau_tilde); tau_tilde = tau_tilde/sum(tau_tilde);
indx = resample(0, tau_tilde', DynareOptions.particle); indx = resample(0, tau_tilde', options_.particle);
StateVectors = StateVectors(:,indx); StateVectors = StateVectors(:,indx);
xparam = fore_xparam(:,indx); xparam = fore_xparam(:,indx);
if pruning if pruning
@ -208,13 +208,13 @@ for t=1:sample_size
for i=1:number_of_particles for i=1:number_of_particles
info = 12042009; info = 12042009;
counter=0; counter=0;
while info(1) && counter <DynareOptions.particle.liu_west_max_resampling_tries while info(1) && counter <options_.particle.liu_west_max_resampling_tries
counter=counter+1; counter=counter+1;
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1); candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
if all(candidate>=bounds.lb) && all(candidate<=bounds.ub) if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
% model resolution for new parameters particles % model resolution for new parameters particles
[info, Model, DynareOptions, DynareResults, ReducedForm] = ... [info, M_, options_, oo_, ReducedForm] = ...
solve_model_for_online_filter(false, candidate, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults) ; solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_) ;
if ~info(1) if ~info(1)
xparam(:,i) = candidate ; xparam(:,i) = candidate ;
steadystate = ReducedForm.steadystate; steadystate = ReducedForm.steadystate;
@ -263,23 +263,23 @@ for t=1:sample_size
% compute particles likelihood contribution % compute particles likelihood contribution
yhat = bsxfun(@minus,StateVectors(:,i), state_variables_steady_state); yhat = bsxfun(@minus,StateVectors(:,i), state_variables_steady_state);
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
else else
if pruning if pruning
yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state_); yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state_);
if order == 2 if order == 2
[tmp, tmp_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2); [tmp, tmp_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, options_.threads.local_state_space_iteration_2);
elseif order == 3 elseif order == 3
[tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning); [tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, options_.threads.local_state_space_iteration_3, pruning);
else else
error('Pruning is not available for orders > 3'); error('Pruning is not available for orders > 3');
end end
StateVectors_(:,i) = tmp_(mf0_,:); StateVectors_(:,i) = tmp_(mf0_,:);
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, options_.threads.local_state_space_iteration_2);
elseif order == 3 elseif order == 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning); tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, options_.threads.local_state_space_iteration_3, pruning);
else else
error('Order > 3: use_k_order_solver should be set to true'); error('Order > 3: use_k_order_solver should be set to true');
end end
@ -290,8 +290,8 @@ for t=1:sample_size
wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1))); wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1)));
end end
end end
if counter==DynareOptions.particle.liu_west_max_resampling_tries if counter==options_.particle.liu_west_max_resampling_tries
fprintf('\nLiu & West particle filter: I haven''t been able to solve the model in %u tries.\n',DynareOptions.particle.liu_west_max_resampling_tries) fprintf('\nLiu & West particle filter: I haven''t been able to solve the model in %u tries.\n',options_.particle.liu_west_max_resampling_tries)
fprintf('Liu & West particle filter: The last error message was: %s\n',get_error_message(info)) fprintf('Liu & West particle filter: The last error message was: %s\n',get_error_message(info))
fprintf('Liu & West particle filter: You can try to increase liu_west_max_resampling_tries, but most\n') fprintf('Liu & West particle filter: You can try to increase liu_west_max_resampling_tries, but most\n')
fprintf('Liu & West particle filter: likely there is an issue with the model.\n') fprintf('Liu & West particle filter: likely there is an issue with the model.\n')
@ -301,14 +301,14 @@ for t=1:sample_size
end end
% normalization % normalization
weights = wtilde/sum(wtilde); weights = wtilde/sum(wtilde);
if variance_update && (neff(weights)<DynareOptions.particle.resampling.threshold*sample_size) if variance_update && (neff(weights)<options_.particle.resampling.threshold*sample_size)
variance_update = false; variance_update = false;
end end
% final resampling (not advised) % final resampling (not advised)
if second_resample if second_resample
[~, idmode] = max(weights); [~, idmode] = max(weights);
mode_xparam(:,t) = xparam(:,idmode); mode_xparam(:,t) = xparam(:,idmode);
indx = resample(0, weights,DynareOptions.particle); indx = resample(0, weights,options_.particle);
StateVectors = StateVectors(:,indx) ; StateVectors = StateVectors(:,indx) ;
if pruning if pruning
StateVectors_ = StateVectors_(:,indx); StateVectors_ = StateVectors_(:,indx);
@ -372,24 +372,24 @@ pmedian = median_xparam(:,sample_size) ;
covariance = mat_var_cov; covariance = mat_var_cov;
%% Plot parameters trajectory %% Plot parameters trajectory
TeX = DynareOptions.TeX; TeX = options_.TeX;
nr = ceil(sqrt(number_of_parameters)) ; nr = ceil(sqrt(number_of_parameters)) ;
nc = floor(sqrt(number_of_parameters)); nc = floor(sqrt(number_of_parameters));
nbplt = 1 ; nbplt = 1 ;
if TeX if TeX
fidTeX = fopen([Model.fname '_param_traj.tex'],'w'); fidTeX = fopen([M_.fname '_param_traj.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n'); fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']); fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n'); fprintf(fidTeX,' \n');
end end
for plt = 1:nbplt for plt = 1:nbplt
hh_fig = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories'); hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Trajectories');
for k=1:length(pmean) for k=1:length(pmean)
subplot(nr,nc,k) subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions); [name,texname] = get_the_name(k,TeX,M_,estim_params_,options_);
% Draw the surface for an interval containing 95% of the particles. % Draw the surface for an interval containing 95% of the particles.
area(1:sample_size, ub95_xparam(k,:), 'FaceColor', [.9 .9 .9], 'BaseValue', min(lb95_xparam(k,:))); area(1:sample_size, ub95_xparam(k,:), 'FaceColor', [.9 .9 .9], 'BaseValue', min(lb95_xparam(k,:)));
hold on hold on
@ -405,12 +405,12 @@ for plt = 1:nbplt
axis tight axis tight
drawnow drawnow
end end
dyn_saveas(hh_fig, [Model.fname '_param_traj' int2str(plt)], DynareOptions.nodisplay, DynareOptions.graph_format); dyn_saveas(hh_fig, [M_.fname '_param_traj' int2str(plt)], options_.nodisplay, options_.graph_format);
if TeX if TeX
% TeX eps loader file % TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n'); fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n'); fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParamTraj%s}\n',Model.fname,int2str(plt)); fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParamTraj%s}\n',M_.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Parameters trajectories.}'); fprintf(fidTeX,'\\caption{Parameters trajectories.}');
fprintf(fidTeX,'\\label{Fig:ParametersPlots:%s}\n',int2str(plt)); fprintf(fidTeX,'\\label{Fig:ParametersPlots:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,'\\end{figure}\n');
@ -423,10 +423,10 @@ number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
bandwidth = 0; % Rule of thumb optimal bandwidth parameter. bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation. kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation.
for plt = 1:nbplt for plt = 1:nbplt
hh_fig = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities'); hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Densities');
for k=1:length(pmean) for k=1:length(pmean)
subplot(nr,nc,k) subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions); [name,texname] = get_the_name(k,TeX,M_,estim_params_,options_);
optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function); optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ... [density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ...
number_of_particles, optimal_bandwidth, kernel_function); number_of_particles, optimal_bandwidth, kernel_function);
@ -441,8 +441,8 @@ for plt = 1:nbplt
axis tight axis tight
drawnow drawnow
end end
dyn_saveas(hh_fig,[ Model.fname '_param_density' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format); dyn_saveas(hh_fig,[ M_.fname '_param_density' int2str(plt) ],options_.nodisplay,options_.graph_format);
if TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format))) if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
% TeX eps loader file % TeX eps loader file
fprintf(fidTeX, '\\begin{figure}[H]\n'); fprintf(fidTeX, '\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n'); fprintf(fidTeX,'\\centering \n');

8
matlab/nonlinear-filters/sequential_importance_particle_filter.m
View File

@ -1,4 +1,4 @@
function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model) function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, options_, M_)
% Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling). % Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
@ -37,7 +37,7 @@ steadystate = ReducedForm.steadystate;
constant = ReducedForm.constant; constant = ReducedForm.constant;
state_variables_steady_state = ReducedForm.state_variables_steady_state; state_variables_steady_state = ReducedForm.state_variables_steady_state;
order = DynareOptions.order; order = options_.order;
% Set persistent variables (if needed). % Set persistent variables (if needed).
if isempty(init_flag) if isempty(init_flag)
@ -101,7 +101,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = chol(Q)'; Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn(). % Set seed for randn().
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default'); options_=set_dynare_seed_local_options(options_,'default');
% Initialization of the weights across particles. % Initialization of the weights across particles.
weights = ones(1,number_of_particles)/number_of_particles ; weights = ones(1,number_of_particles)/number_of_particles ;
@ -139,7 +139,7 @@ for t=1:sample_size
end end
else else
if ReducedForm.use_k_order_solver if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
else else
if order == 2 if order == 2
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);

108
matlab/nonlinear-filters/solve_model_for_online_filter.m
View File

@ -1,23 +1,23 @@
function [info, Model, DynareOptions, DynareResults, ReducedForm] = ... function [info, M_, options_, oo_, ReducedForm] = ...
solve_model_for_online_filter(setinitialcondition, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults) solve_model_for_online_filter(setinitialcondition, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_)
% Solves the dsge model for an particular parameters set. % Solves the dsge model for an particular parameters set.
% %
% INPUTS % INPUTS
% - setinitialcondition [logical] return initial condition if true. % - setinitialcondition [logical] return initial condition if true.
% - xparam1 [double] n×1 vector, parameter values. % - xparam1 [double] n×1 vector, parameter values.
% - DynareDataset [struct] Dataset for estimation (dataset_). % - dataset_ [struct] Dataset for estimation.
% - DynareOptions [struct] Dynare options (options_). % - options_ [struct] Dynare options.
% - Model [struct] Model description (M_). % - M_ [struct] Model description.
% - EstimatedParameters [struct] Estimated parameters (estim_params_). % - estim_params_ [struct] Estimated parameters.
% - BayesInfo [struct] Prior definition (bayestopt_). % - bayestopt_ [struct] Prior definition.
% - DynareResults [struct] Dynare results (oo_). % - oo_ [struct] Dynare results.
% %
% OUTPUTS % OUTPUTS
% - info [integer] scalar, nonzero if any problem occur when computing the reduced form. % - info [integer] scalar, nonzero if any problem occur when computing the reduced form.
% - Model [struct] Model description (M_). % - M_ [struct] M_ description.
% - DynareOptions [struct] Dynare options (options_). % - options_ [struct] Dynare options.
% - DynareResults [struct] Dynare results (oo_). % - oo_ [struct] Dynare results.
% - ReducedForm [struct] Reduced form model. % - ReducedForm [struct] Reduced form model.
% Copyright © 2013-2023 Dynare Team % Copyright © 2013-2023 Dynare Team
@ -58,27 +58,27 @@ if any(xparam1>bounds.ub)
end end
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H). % Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
Q = Model.Sigma_e; Q = M_.Sigma_e;
H = Model.H; H = M_.H;
for i=1:EstimatedParameters.nvx for i=1:estim_params_.nvx
k =EstimatedParameters.var_exo(i,1); k =estim_params_.var_exo(i,1);
Q(k,k) = xparam1(i)*xparam1(i); Q(k,k) = xparam1(i)*xparam1(i);
end end
offset = EstimatedParameters.nvx; offset = estim_params_.nvx;
if EstimatedParameters.nvn if estim_params_.nvn
for i=1:EstimatedParameters.nvn for i=1:estim_params_.nvn
H(i,i) = xparam1(i+offset)*xparam1(i+offset); H(i,i) = xparam1(i+offset)*xparam1(i+offset);
end end
offset = offset+EstimatedParameters.nvn; offset = offset+estim_params_.nvn;
else else
H = zeros(size(DynareDataset.data, 2)); H = zeros(size(dataset_.data, 2));
end end
% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite. % Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
if EstimatedParameters.ncx if estim_params_.ncx
for i=1:EstimatedParameters.ncx for i=1:estim_params_.ncx
k1 =EstimatedParameters.corrx(i,1); k1 =estim_params_.corrx(i,1);
k2 =EstimatedParameters.corrx(i,2); k2 =estim_params_.corrx(i,2);
Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2)); Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
Q(k2,k1) = Q(k1,k2); Q(k2,k1) = Q(k1,k2);
end end
@ -89,13 +89,13 @@ if EstimatedParameters.ncx
info = 43; info = 43;
return return
end end
offset = offset+EstimatedParameters.ncx; offset = offset+estim_params_.ncx;
end end
% Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite. % Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite.
if EstimatedParameters.ncn if estim_params_.ncn
corrn_observable_correspondence = EstimatedParameters.corrn_observable_correspondence; corrn_observable_correspondence = estim_params_.corrn_observable_correspondence;
for i=1:EstimatedParameters.ncn for i=1:estim_params_.ncn
k1 = corrn_observable_correspondence(i,1); k1 = corrn_observable_correspondence(i,1);
k2 = corrn_observable_correspondence(i,2); k2 = corrn_observable_correspondence(i,2);
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2)); H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
@ -108,25 +108,25 @@ if EstimatedParameters.ncn
info = 44; info = 44;
return return
end end
offset = offset+EstimatedParameters.ncn; offset = offset+estim_params_.ncn;
end end
% Update estimated structural parameters in Mode.params. % Update estimated structural parameters in Mode.params.
if EstimatedParameters.np > 0 if estim_params_.np > 0
Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end); M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end end
% Update Model.Sigma_e and Model.H. % Update M_.Sigma_e and M_.H.
Model.Sigma_e = Q; M_.Sigma_e = Q;
Model.H = H; M_.H = H;
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
% 2. call model setup & reduction program % 2. call model setup & reduction program
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
warning('off', 'MATLAB:nearlySingularMatrix') warning('off', 'MATLAB:nearlySingularMatrix')
[~, ~, ~, info, DynareResults.dr, Model.params] = ... [~, ~, ~, info, oo_.dr, M_.params] = ...
dynare_resolve(Model, DynareOptions, DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state, 'restrict'); dynare_resolve(M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, 'restrict');
warning('on', 'MATLAB:nearlySingularMatrix') warning('on', 'MATLAB:nearlySingularMatrix')
if info(1)~=0 if info(1)~=0
@ -137,12 +137,12 @@ if info(1)~=0
end end
% Get decision rules and transition equations. % Get decision rules and transition equations.
dr = DynareResults.dr; dr = oo_.dr;
% Set persistent variables (first call). % Set persistent variables (first call).
if isempty(init_flag) if isempty(init_flag)
mf0 = BayesInfo.mf0; mf0 = bayestopt_.mf0;
mf1 = BayesInfo.mf1; mf1 = bayestopt_.mf1;
restrict_variables_idx = dr.restrict_var_list; restrict_variables_idx = dr.restrict_var_list;
state_variables_idx = restrict_variables_idx(mf0); state_variables_idx = restrict_variables_idx(mf0);
number_of_state_variables = length(mf0); number_of_state_variables = length(mf0);
@ -155,17 +155,17 @@ if nargout>4
ReducedForm.ghx = dr.ghx(restrict_variables_idx,:); ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
ReducedForm.ghu = dr.ghu(restrict_variables_idx,:); ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx)); ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
if DynareOptions.order==2 if options_.order==2
ReducedForm.use_k_order_solver = false; ReducedForm.use_k_order_solver = false;
ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:); ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:); ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:); ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx); ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:); ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:);
elseif DynareOptions.order>=3 elseif options_.order>=3
ReducedForm.use_k_order_solver = true; ReducedForm.use_k_order_solver = true;
ReducedForm.dr = dr; ReducedForm.dr = dr;
ReducedForm.udr = folded_to_unfolded_dr(dr, Model, DynareOptions); ReducedForm.udr = folded_to_unfolded_dr(dr, M_, options_);
else else
n_states=size(dr.ghx,2); n_states=size(dr.ghx,2);
n_shocks=size(dr.ghu,2); n_shocks=size(dr.ghu,2);
@ -184,28 +184,28 @@ end
% Set initial condition % Set initial condition
if setinitialcondition if setinitialcondition
switch DynareOptions.particle.initialization switch options_.particle.initialization
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model. case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0); StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
[A,B] = kalman_transition_matrix(dr,dr.restrict_var_list,dr.restrict_columns); [A,B] = kalman_transition_matrix(dr,dr.restrict_var_list,dr.restrict_columns);
StateVectorVariance = lyapunov_symm(A, B*ReducedForm.Q*B', DynareOptions.lyapunov_fixed_point_tol, ... StateVectorVariance = lyapunov_symm(A, B*ReducedForm.Q*B', options_.lyapunov_fixed_point_tol, ...
DynareOptions.qz_criterium, DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); options_.qz_criterium, options_.lyapunov_complex_threshold, [], options_.debug);
StateVectorVariance = StateVectorVariance(mf0,mf0); StateVectorVariance = StateVectorVariance(mf0,mf0);
case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model). case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0); StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
old_DynareOptionsperiods = DynareOptions.periods; old_DynareOptionsperiods = options_.periods;
DynareOptions.periods = 5000; options_.periods = 5000;
old_DynareOptionspruning = DynareOptions.pruning; old_DynareOptionspruning = options_.pruning;
DynareOptions.pruning = DynareOptions.particle.pruning; options_.pruning = options_.particle.pruning;
y_ = simult(dr.ys, dr, Model, DynareOptions); y_ = simult(dr.ys, dr, M_, options_);
y_ = y_(dr.order_var(state_variables_idx),2001:DynareOptions.periods); y_ = y_(dr.order_var(state_variables_idx),2001:options_.periods);
StateVectorVariance = cov(y_'); StateVectorVariance = cov(y_');
DynareOptions.periods = old_DynareOptionsperiods; options_.periods = old_DynareOptionsperiods;
DynareOptions.pruning = old_DynareOptionspruning; options_.pruning = old_DynareOptionspruning;
clear('old_DynareOptionsperiods','y_'); clear('old_DynareOptionsperiods','y_');
case 3% Initial state vector covariance is a diagonal matrix. case 3% Initial state vector covariance is a diagonal matrix.
StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0); StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables); StateVectorVariance = options_.particle.initial_state_prior_std*eye(number_of_state_variables);
otherwise otherwise
error('Unknown initialization option!') error('Unknown initialization option!')
end end

24
matlab/print_moments_implied_prior.m
View File

@ -1,9 +1,9 @@
function print_moments_implied_prior(ModelInfo, mm, vm, mv, vv) function print_moments_implied_prior(M_, mm, vm, mv, vv)
%function print_moments_implied_prior(ModelInfo, mm, vm, mv, vv) %function print_moments_implied_prior(M_, mm, vm, mv, vv)
% This routine prints in the command window some descriptive statistics % This routine prints in the command window some descriptive statistics
% about the endogenous variables implied prior moments. % about the endogenous variables implied prior moments.
% Inputs: % Inputs:
% - ModelInfo [structure] Dynare's model structure % - M_ [structure] Dynare's model structure
% - mm [endo_nbr*1] mean first moments of the endogenous % - mm [endo_nbr*1] mean first moments of the endogenous
% variables % variables
% - vm [endo_nbr*1] variance of the first moments of the % - vm [endo_nbr*1] variance of the first moments of the
@ -14,7 +14,7 @@ function print_moments_implied_prior(ModelInfo, mm, vm, mv, vv)
% endogenous variables % endogenous variables
% Copyright © 2016-2018 Dynare Team % Copyright © 2016-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -39,14 +39,14 @@ disp(printline(64, '-'))
T1 = 'VARIABLE '; T1 = 'VARIABLE ';
T2 = sprintf('Prior mean \t Prior st. dev.'); T2 = sprintf('Prior mean \t Prior st. dev.');
for i=1:ModelInfo.orig_endo_nbr for i=1:M_.orig_endo_nbr
Name = ModelInfo.endo_names{i}; Name = M_.endo_names{i};
T1 = strvcat(T1, Name); T1 = strvcat(T1, Name);
str = sprintf(' %6.4f \t %6.4f', mm(i), sqrt(vm(i))); str = sprintf(' %6.4f \t %6.4f', mm(i), sqrt(vm(i)));
T2 = strvcat(T2, str); T2 = strvcat(T2, str);
end end
T0 = repmat(' ', ModelInfo.orig_endo_nbr+1, 1); T0 = repmat(' ', M_.orig_endo_nbr+1, 1);
TT = [T1, T0, T2]; TT = [T1, T0, T2];
l0 = printline(size(TT, 2)+1, '-'); l0 = printline(size(TT, 2)+1, '-');
@ -64,10 +64,10 @@ T1b = 'VARIABLE-2';
T2a = 'Prior mean'; T2a = 'Prior mean';
T2b = 'Prior st.dev.'; T2b = 'Prior st.dev.';
for i=1:ModelInfo.orig_endo_nbr for i=1:M_.orig_endo_nbr
for j=i:ModelInfo.orig_endo_nbr for j=i:M_.orig_endo_nbr
Name1 = ModelInfo.endo_names{i}; Name1 = M_.endo_names{i};
Name2 = ModelInfo.endo_names{j}; Name2 = M_.endo_names{j};
T1a = strvcat(T1a, Name1); T1a = strvcat(T1a, Name1);
T1b = strvcat(T1b, Name2); T1b = strvcat(T1b, Name2);
sta = sprintf('%12.8f', mv(i,j)); sta = sprintf('%12.8f', mv(i,j));
@ -77,7 +77,7 @@ for i=1:ModelInfo.orig_endo_nbr
end end
end end
T0 = repmat(' ', ModelInfo.orig_endo_nbr*(ModelInfo.orig_endo_nbr+1)/2+1, 1); T0 = repmat(' ', M_.orig_endo_nbr*(M_.orig_endo_nbr+1)/2+1, 1);
TT = [T1a, T0, T1b, T0, T2a, T0, T2b]; TT = [T1a, T0, T1b, T0, T2a, T0, T2b];
l0 = printline(size(TT, 2)+1, '-'); l0 = printline(size(TT, 2)+1, '-');

42
matlab/prior_draw.m
View File

@ -1,5 +1,5 @@
function pdraw = prior_draw(BayesInfo, prior_trunc, uniform) function pdraw = prior_draw(bayestopt_, prior_trunc, uniform)
% pdraw = prior_draw(bayestopt_, prior_trunc, uniform)
% This function generate one draw from the joint prior distribution and % This function generate one draw from the joint prior distribution and
% allows sampling uniformly from the prior support (uniform==1 when called with init==1) % allows sampling uniformly from the prior support (uniform==1 when called with init==1)
% %
@ -48,18 +48,18 @@ persistent uniform_index gaussian_index gamma_index beta_index inverse_gamma_1_i
persistent uniform_draws gaussian_draws gamma_draws beta_draws inverse_gamma_1_draws inverse_gamma_2_draws weibull_draws persistent uniform_draws gaussian_draws gamma_draws beta_draws inverse_gamma_1_draws inverse_gamma_2_draws weibull_draws
if nargin>0 if nargin>0
p6 = BayesInfo.p6; p6 = bayestopt_.p6;
p7 = BayesInfo.p7; p7 = bayestopt_.p7;
p3 = BayesInfo.p3; p3 = bayestopt_.p3;
p4 = BayesInfo.p4; p4 = bayestopt_.p4;
bounds = prior_bounds(BayesInfo, prior_trunc); bounds = prior_bounds(bayestopt_, prior_trunc);
lb = bounds.lb; lb = bounds.lb;
ub = bounds.ub; ub = bounds.ub;
number_of_estimated_parameters = length(p6); number_of_estimated_parameters = length(p6);
if nargin>2 && uniform if nargin>2 && uniform
prior_shape = repmat(5,number_of_estimated_parameters,1); prior_shape = repmat(5,number_of_estimated_parameters,1);
else else
prior_shape = BayesInfo.pshape; prior_shape = bayestopt_.pshape;
end end
beta_index = find(prior_shape==1); beta_index = find(prior_shape==1);
if isempty(beta_index) if isempty(beta_index)
@ -238,23 +238,23 @@ for i=1:14
end end
end end
BayesInfo.pshape = p0; bayestopt_.pshape = p0;
BayesInfo.p1 = p1; bayestopt_.p1 = p1;
BayesInfo.p2 = p2; bayestopt_.p2 = p2;
BayesInfo.p3 = p3; bayestopt_.p3 = p3;
BayesInfo.p4 = p4; bayestopt_.p4 = p4;
BayesInfo.p5 = p5; bayestopt_.p5 = p5;
BayesInfo.p6 = p6; bayestopt_.p6 = p6;
BayesInfo.p7 = p7; bayestopt_.p7 = p7;
ndraws = 1e5; ndraws = 1e5;
m0 = BayesInfo.p1; %zeros(14,1); m0 = bayestopt_.p1; %zeros(14,1);
v0 = diag(BayesInfo.p2.^2); %zeros(14); v0 = diag(bayestopt_.p2.^2); %zeros(14);
% Call the tested routine % Call the tested routine
try try
% Initialization of prior_draws. % Initialization of prior_draws.
prior_draw(BayesInfo, prior_trunc, false); prior_draw(bayestopt_, prior_trunc, false);
% Do simulations in a loop and estimate recursively the mean and teh variance. % Do simulations in a loop and estimate recursively the mean and teh variance.
for i = 1:ndraws for i = 1:ndraws
draw = transpose(prior_draw()); draw = transpose(prior_draw());
@ -269,8 +269,8 @@ catch
end end
if t(1) if t(1)
t(2) = all(abs(m0-BayesInfo.p1)<3e-3); t(2) = all(abs(m0-bayestopt_.p1)<3e-3);
t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3)); t(3) = all(all(abs(v0-diag(bayestopt_.p2.^2))<5e-3));
end end
T = all(t); T = all(t);
%@eof:1 %@eof:1

16
matlab/realtime_shock_decomposition.m
View File

@ -135,8 +135,8 @@ for j=presample+1:nobs
% evalin('base',['options_.nobs=' int2str(j) ';']) % evalin('base',['options_.nobs=' int2str(j) ';'])
options_.nobs=j; options_.nobs=j;
if isequal(fast_realtime,0) if isequal(fast_realtime,0)
[oo,M_,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_); [oo_local,M_,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
gend = size(oo.SmoothedShocks.(M_.exo_names{1}),1); gend = size(oo_local.SmoothedShocks.(M_.exo_names{1}),1);
else else
if j<min(fast_realtime) && gend0<j if j<min(fast_realtime) && gend0<j
options_.nobs=min(fast_realtime); options_.nobs=min(fast_realtime);
@ -146,10 +146,10 @@ for j=presample+1:nobs
end end
if ismember(j,fast_realtime) && gend0<j if ismember(j,fast_realtime) && gend0<j
[oo,M_,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_); [oo_local,M_,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
gend = size(oo.SmoothedShocks.(M_.exo_names{1}),1); gend = size(oo_local.SmoothedShocks.(M_.exo_names{1}),1);
gend0 = gend; gend0 = gend;
oo0=oo; oo0=oo_local;
Smoothed_Variables_deviation_from_mean0=Smoothed_Variables_deviation_from_mean; Smoothed_Variables_deviation_from_mean0=Smoothed_Variables_deviation_from_mean;
else else
if j>gend0 if j>gend0
@ -164,13 +164,13 @@ for j=presample+1:nobs
end end
gend = j; gend = j;
oo=oo0; oo_local=oo0;
Smoothed_Variables_deviation_from_mean = Smoothed_Variables_deviation_from_mean0(:,1:gend); Smoothed_Variables_deviation_from_mean = Smoothed_Variables_deviation_from_mean0(:,1:gend);
end end
end end
% reduced form % reduced form
dr = oo.dr; dr = oo_local.dr;
% data reordering % data reordering
order_var = dr.order_var; order_var = dr.order_var;
@ -194,7 +194,7 @@ for j=presample+1:nobs
% initialization % initialization
epsilon=NaN(nshocks,gend); epsilon=NaN(nshocks,gend);
for i = 1:nshocks for i = 1:nshocks
epsilon(i,:) = oo.SmoothedShocks.(M_.exo_names{i})(1:gend); epsilon(i,:) = oo_local.SmoothedShocks.(M_.exo_names{i})(1:gend);
end end
epsilon=[epsilon zeros(nshocks,forecast_)]; epsilon=[epsilon zeros(nshocks,forecast_)];

15
matlab/remove_aux_variables_from_expression.m
View File

@ -1,6 +1,7 @@
function expression = remove_aux_variables_from_expression(expression, DynareModel) function expression = remove_aux_variables_from_expression(expression, M_)
% expression = remove_aux_variables_from_expression(expression, M_)
% Copyright © 2022 Dynare Team % Copyright © 2022-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -19,17 +20,17 @@ function expression = remove_aux_variables_from_expression(expression, DynareMod
% Get list of endogenous variables in expression % Get list of endogenous variables in expression
list_of_words = regexp(expression, '\<\w*\>', 'match'); list_of_words = regexp(expression, '\<\w*\>', 'match');
list_of_words = setdiff(list_of_words, DynareModel.param_names); list_of_words = setdiff(list_of_words, M_.param_names);
list_of_words = setdiff(list_of_words, DynareModel.exo_names); list_of_words = setdiff(list_of_words, M_.exo_names);
isnotanumber = isnan(str2double(list_of_words)); isnotanumber = isnan(str2double(list_of_words));
list_of_words = list_of_words(isnotanumber); list_of_words = list_of_words(isnotanumber);
list_of_words = setdiff(list_of_words, {'diff','log','exp'}); list_of_words = setdiff(list_of_words, {'diff','log','exp'});
for i=1:length(list_of_words) for i=1:length(list_of_words)
id = find(strcmp(list_of_words{i}, DynareModel.endo_names)); id = find(strcmp(list_of_words{i}, M_.endo_names));
if isempty(id) || id<=DynareModel.orig_endo_nbr if isempty(id) || id<=M_.orig_endo_nbr
continue continue
end end
auxinfo = DynareModel.aux_vars(get_aux_variable_id(id)); auxinfo = M_.aux_vars(get_aux_variable_id(id));
expression = regexprep(expression, sprintf('\\<%s\\>', list_of_words{i}), auxinfo.orig_expr); expression = regexprep(expression, sprintf('\\<%s\\>', list_of_words{i}), auxinfo.orig_expr);
end end

16
matlab/set_exogenous_variables_for_simulation.m
View File

@ -1,17 +1,17 @@
function DynareModel = set_exogenous_variables_for_simulation(DynareModel) function M_ = set_exogenous_variables_for_simulation(M_)
% M_ = set_exogenous_variables_for_simulation(M_)
% Appends the list of observed exogenous variables in Dynare's model structure (if any). % Appends the list of observed exogenous variables in Dynare's model structure (if any).
% %
% INPUTS % INPUTS
% - DynareModel [struct] Dynare's model global structure, M_. % - M_ [struct] Dynare's model global structure.
% %
% OUTPUTS % OUTPUTS
% - DynareModel [struct] Dynare's model global structure, M_. % - M_ [struct] Dynare's model global structure.
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2019 Dynare Team % Copyright © 2019-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -28,8 +28,8 @@ function DynareModel = set_exogenous_variables_for_simulation(DynareModel)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if isfield(DynareModel, 'exo_partitions') if isfield(M_, 'exo_partitions')
if isfield(DynareModel.exo_partitions, 'used') if isfield(M_.exo_partitions, 'used')
DynareModel.simulation_exo_names = DynareModel.exo_names(~strcmpi('estimationonly', DynareModel.exo_partitions.used)); M_.simulation_exo_names = M_.exo_names(~strcmpi('estimationonly', M_.exo_partitions.used));
end end
end end

16
matlab/set_observed_exogenous_variables.m
View File

@ -1,17 +1,17 @@
function DynareModel = set_observed_exogenous_variables(DynareModel) function M_ = set_observed_exogenous_variables(M_)
% M_ = set_observed_exogenous_variables(M_)
% Appends the list of observed exogenous variables in Dynare's model structure (if any). % Appends the list of observed exogenous variables in Dynare's model structure (if any).
% %
% INPUTS % INPUTS
% - DynareModel [struct] Dynare's model global structure, M_. % - M_ [struct] Dynare's model global structure.
% %
% OUTPUTS % OUTPUTS
% - DynareModel [struct] Dynare's model global structure, M_. % - M_ [struct] Dynare's model global structure.
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2019 Dynare Team % Copyright © 2019-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -28,8 +28,8 @@ function DynareModel = set_observed_exogenous_variables(DynareModel)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if isfield(DynareModel, 'exo_partitions') if isfield(M_, 'exo_partitions')
if isfield(DynareModel.exo_partitions, 'status') if isfield(M_.exo_partitions, 'status')
DynareModel.observed_exo_names = DynareModel.exo_names(strcmpi('observed', DynareModel.exo_partitions.status)); M_.observed_exo_names = M_.exo_names(strcmpi('observed', M_.exo_partitions.status));
end end
end end

6
matlab/simul_static_model.m
View File

@ -10,8 +10,8 @@ function simulation = simul_static_model(samplesize, innovations)
% - simulation [dseries] Simulated endogenous and exogenous variables. % - simulation [dseries] Simulated endogenous and exogenous variables.
% %
% REMARKS % REMARKS
% [1] The innovations used for the simulation are saved in DynareOutput.exo_simul, and the resulting paths for the endogenous % [1] The innovations used for the simulation are saved in oo_.exo_simul, and the resulting paths for the endogenous
% variables are saved in DynareOutput.endo_simul. % variables are saved in oo_.endo_simul.
% [2] The last input argument is not mandatory. If absent we use random draws and rescale them with the informations provided % [2] The last input argument is not mandatory. If absent we use random draws and rescale them with the informations provided
% through the shocks block. % through the shocks block.
@ -62,7 +62,7 @@ if nargin<2 || isempty(innovations)
otherwise otherwise
error('%s distribution for the structural innovations is not (yet) implemented!', options_.bnlms.innovation_distribution) error('%s distribution for the structural innovations is not (yet) implemented!', options_.bnlms.innovation_distribution)
end end
% Put the simulated innovations in DynareOutput.exo_simul. % Put the simulated innovations in oo_.exo_simul.
oo_.exo_simul = zeros(samplesize, number_of_shocks); oo_.exo_simul = zeros(samplesize, number_of_shocks);
oo_.exo_simul(:,positive_var_indx) = oo_.bnlms.shocks; oo_.exo_simul(:,positive_var_indx) = oo_.bnlms.shocks;
innovations = []; innovations = [];

10
matlab/update_all_parameters_in_workspace.m
View File

@ -1,8 +1,8 @@
function update_all_parameters_in_workspace(DynareModel) function update_all_parameters_in_workspace(M_)
% update_all_parameters_in_workspace(M_)
% Updates all parameter values in Matlab/Octave base workspace. % Updates all parameter values in Matlab/Octave base workspace.
% Copyright © 2018 Dynare Team % Copyright © 2018-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -19,6 +19,6 @@ function update_all_parameters_in_workspace(DynareModel)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
for i=1:length(DynareModel.params) for i=1:length(M_.params)
assignin('base', DynareModel.param_names{i}, DynareModel.params(i)); assignin('base', M_.param_names{i}, M_.params(i));
end end

7
matlab/utilities/estimation/get_estimated_parameters_indices.m
View File

@ -1,6 +1,7 @@
function ipnames = get_estimated_parameters_indices(params, pnames, eqname, DynareModel) function ipnames = get_estimated_parameters_indices(params, pnames, eqname, M_)
% ipnames = get_estimated_parameters_indices(params, pnames, eqname, M_)
% Copyright © 2021 Dynare Team % Copyright © 2021-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -36,5 +37,5 @@ end
ipnames = zeros(size(list_of_parameters)); ipnames = zeros(size(list_of_parameters));
for i=1:length(ipnames) for i=1:length(ipnames)
ipnames(i) = find(strcmp(list_of_parameters{i}, DynareModel.param_names)); ipnames(i) = find(strcmp(list_of_parameters{i}, M_.param_names));
end end

2
matlab/utilities/estimation/rewrite_equation_with_tables.m
View File

@ -36,7 +36,7 @@ objTypes = objTypes(I);
for i=1:length(objNames) for i=1:length(objNames)
switch objTypes(i) switch objTypes(i)
case 1 case 1
expression = strrep(expression, objNames{i}, sprintf('DynareModel.params(%u)', objIndex(i))); expression = strrep(expression, objNames{i}, sprintf('M_.params(%u)', objIndex(i)));
case 2 case 2
k = find(strcmp(objNames{i}, data.name)); k = find(strcmp(objNames{i}, data.name));
if isempty(k) if isempty(k)

14
matlab/utilities/estimation/write_residuals_routine.m
View File

@ -1,8 +1,8 @@
function write_residuals_routine(lhs, rhs, eqname, ipnames, DynareModel, pacmodl) function write_residuals_routine(lhs, rhs, eqname, ipnames, M_, pacmodl)
% write_residuals_routine(lhs, rhs, eqname, ipnames, M_, pacmodl)
% Creates a routine for evaluating the residuals of the nonlinear equation. % Creates a routine for evaluating the residuals of the nonlinear equation.
% Copyright © 2021 Dynare Team % Copyright © 2021-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -20,18 +20,18 @@ function write_residuals_routine(lhs, rhs, eqname, ipnames, DynareModel, pacmodl
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
fun = sprintf('r_%s', eqname); fun = sprintf('r_%s', eqname);
fid = fopen(['+' DynareModel.fname filesep() fun '.m'], 'w'); fid = fopen(['+' M_.fname filesep() fun '.m'], 'w');
fprintf(fid, 'function r = %s(params, data, DynareModel, DynareOutput)\n', fun); fprintf(fid, 'function r = %s(params, data, M_, oo_)\n', fun);
fprintf(fid, '\n'); fprintf(fid, '\n');
fprintf(fid, '%% Evaluates the residuals for equation %s.\n', eqname); fprintf(fid, '%% Evaluates the residuals for equation %s.\n', eqname);
fprintf(fid, '%% File created by Dynare (%s).\n', datetime); fprintf(fid, '%% File created by Dynare (%s).\n', datetime);
fprintf(fid, '\n'); fprintf(fid, '\n');
for i=1:length(ipnames) for i=1:length(ipnames)
fprintf(fid, 'DynareModel.params(%u) = params(%u);\n', ipnames(i), i); fprintf(fid, 'M_.params(%u) = params(%u);\n', ipnames(i), i);
end end
fprintf(fid, '\n'); fprintf(fid, '\n');
if nargin>5 if nargin>5
fprintf(fid, 'DynareModel = pac.update.parameters(''%s'', DynareModel, DynareOutput, false);\n', pacmodl); fprintf(fid, 'M_ = pac.update.parameters(''%s'', M_, oo_, false);\n', pacmodl);
fprintf(fid, '\n'); fprintf(fid, '\n');
end end
fprintf(fid, 'r = %s-(%s);\n', lhs, rhs); fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);

14
matlab/utilities/estimation/write_ssr_routine.m
View File

@ -1,8 +1,8 @@
function write_ssr_routine(lhs, rhs, eqname, ipnames, DynareModel, pacmodl) function write_ssr_routine(lhs, rhs, eqname, ipnames, M_, pacmodl)
% write_ssr_routine(lhs, rhs, eqname, ipnames, M_, pacmodl)
% Creates a routine for evaluating the sum of squared residuals of the nonlinear equation. % Creates a routine for evaluating the sum of squared residuals of the nonlinear equation.
% Copyright © 2021 Dynare Team % Copyright © 2021-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -20,8 +20,8 @@ function write_ssr_routine(lhs, rhs, eqname, ipnames, DynareModel, pacmodl)
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
fun = sprintf('ssr_%s', eqname); fun = sprintf('ssr_%s', eqname);
fid = fopen(['+' DynareModel.fname filesep() fun '.m'], 'w'); fid = fopen(['+' M_.fname filesep() fun '.m'], 'w');
fprintf(fid, 'function [s, fake1, fake2, fake3, fake4] = %s(params, data, DynareModel, DynareOutput)\n', fun); fprintf(fid, 'function [s, fake1, fake2, fake3, fake4] = %s(params, data, M_, oo_)\n', fun);
fprintf(fid, '\n'); fprintf(fid, '\n');
fprintf(fid, '%% Evaluates the sum of square residuals for equation %s.\n', eqname); fprintf(fid, '%% Evaluates the sum of square residuals for equation %s.\n', eqname);
fprintf(fid, '%% File created by Dynare (%s).\n', datetime); fprintf(fid, '%% File created by Dynare (%s).\n', datetime);
@ -32,11 +32,11 @@ fprintf(fid, 'fake3 = [];\n');
fprintf(fid, 'fake4 = [];\n'); fprintf(fid, 'fake4 = [];\n');
fprintf(fid, '\n'); fprintf(fid, '\n');
for i=1:length(ipnames) for i=1:length(ipnames)
fprintf(fid, 'DynareModel.params(%u) = params(%u);\n', ipnames(i), i); fprintf(fid, 'M_.params(%u) = params(%u);\n', ipnames(i), i);
end end
fprintf(fid, '\n'); fprintf(fid, '\n');
if nargin>5 if nargin>5
fprintf(fid, 'DynareModel = pac.update.parameters(''%s'', DynareModel, DynareOutput, false);\n', pacmodl); fprintf(fid, 'M_ = pac.update.parameters(''%s'', M_, oo_, false);\n', pacmodl);
fprintf(fid, '\n'); fprintf(fid, '\n');
end end
fprintf(fid, 'r = %s-(%s);\n', lhs, rhs); fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);

6
matlab/write_latex_prior_table.m
View File

@ -45,7 +45,7 @@ if (size(estim_params_.var_endo,1) || size(estim_params_.corrn,1))
end end
% Fill or update bayestopt_ structure % Fill or update bayestopt_ structure
[xparam1, EstimatedParameters, BayesOptions, lb, ub, Model] = set_prior(estim_params_, M_, options_); [xparam1, estim_params_, BayesOptions, lb, ub, M_] = set_prior(estim_params_, M_, options_);
% Get untruncated bounds % Get untruncated bounds
bounds = prior_bounds(BayesOptions, options_.prior_trunc); bounds = prior_bounds(BayesOptions, options_.prior_trunc);
@ -57,7 +57,7 @@ PriorNames = { 'Beta' , 'Gamma' , 'Gaussian' , 'Inv. Gamma' , 'Uniform' , 'Inv.
if ~exist([M_.dname '/latex'],'dir') if ~exist([M_.dname '/latex'],'dir')
mkdir(M_.dname,'latex'); mkdir(M_.dname,'latex');
end end
fidTeX = fopen([M_.dname, '/latex/' Model.fname '_priors_table.tex'],'w+'); fidTeX = fopen([M_.dname, '/latex/' M_.fname '_priors_table.tex'],'w+');
fprintf(fidTeX,'%% TeX-table generated by Dynare write_latex_prior_table.m.\n'); fprintf(fidTeX,'%% TeX-table generated by Dynare write_latex_prior_table.m.\n');
fprintf(fidTeX,'%% Prior Information\n'); fprintf(fidTeX,'%% Prior Information\n');
fprintf(fidTeX,['%% ' datestr(now,0)]); fprintf(fidTeX,['%% ' datestr(now,0)]);
@ -112,7 +112,7 @@ fprintf(fidTeX,'\\endlastfoot\n');
% Column 8: the upper bound of the interval containing 90% of the prior mass. % Column 8: the upper bound of the interval containing 90% of the prior mass.
PriorIntervals = prior_bounds(BayesOptions,(1-options_.prior_interval)/2) ; PriorIntervals = prior_bounds(BayesOptions,(1-options_.prior_interval)/2) ;
for i=1:size(BayesOptions.name,1) for i=1:size(BayesOptions.name,1)
[tmp,TexName] = get_the_name(i, 1, Model, EstimatedParameters, options_); [tmp,TexName] = get_the_name(i, 1, M_, estim_params_, options_);
PriorShape = PriorNames{ BayesOptions.pshape(i) }; PriorShape = PriorNames{ BayesOptions.pshape(i) };
PriorMean = BayesOptions.p1(i); PriorMean = BayesOptions.p1(i);
PriorMode = BayesOptions.p5(i); PriorMode = BayesOptions.p5(i);