Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 4.3
  • 4.4
  • 4.5
  • 4.6
  • 5.x
  • 6.x
  • asm
  • aux_func
  • clang+openmp
  • dates-and-dseries-improvements
  • declare_vars_in_model_block
  • dmm
  • dragonfly
  • dynare_minreal
  • eigen
  • error_msg_undeclared_model_vars
  • estim_params
  • exo_steady_state
  • gpm-optimal-policy
  • julia
  • madysson
  • master
  • mex-GetPowerDeriv
  • penalty
  • separateM_
  • slice
  • sphinx-doc-experimental
  • static_aux_vars
  • time-varying-information-set
  • various_fixes
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
  • 4.5.7
  • 4.6-beta1
  • 4.6.0
  • 4.6.0-rc1
  • 4.6.0-rc2
  • 4.6.1
  • 4.6.2
  • 4.6.3
  • 4.6.4
  • 4.7-beta1
  • 4.7-beta2
  • 4.7-beta3
  • 5.0
  • 5.0-rc1
  • 5.1
  • 5.2
  • 5.3
  • 5.4
  • 5.5
  • 6-beta1
  • 6-beta2
  • 6.0
  • 6.1
  • 6.2
  • 6.3
  • 6.4
91 results

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • ebenetce/dynare
  • chskcau/dynare-doc-fixes
28 results
Select Git revision
  • bayesian_irf_matching
  • cet-octave-one-more-time
  • docker-meson
  • identificationME
  • irf-matching-example
  • irf_match_perfect_foresight
  • main
  • mom_fname
  • new-desktop-matlab
  • optimizers
  • pskf
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
46 results
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 1980 additions and 784 deletions
matlab/+mom/print_info_on_estimation_settings.m
View file @ 4a32acc2
function print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters)
% function print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters)
function print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, do_bayesian_estimation)
% print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, do_bayesian_estimation)
% -------------------------------------------------------------------------
% Print information on the method of moments estimation settings to the console
% =========================================================================
% -------------------------------------------------------------------------
% INPUTS
% options_mom_ [struct] Options for the method of moments estimation
% number_of_estimated_parameters [integer] Number of estimated parameters
% options_mom_ [struct] options for the method of moments estimation
% number_of_estimated_parameters [integer] number of estimated parameters
% do_bayesian_estimation [boolean] true if the estimation is Bayesian
% -------------------------------------------------------------------------
% OUTPUT
% No output, just displays the chosen settings
......@@ -15,7 +16,8 @@ function print_info_on_estimation_settings(options_mom_, number_of_estimated_par
% -------------------------------------------------------------------------
% This function calls
% o skipline
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
......@@ -32,7 +34,8 @@ function print_info_on_estimation_settings(options_mom_, number_of_estimated_par
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
fprintf('\n---------------------------------------------------\n')
if strcmp(options_mom_.mom.mom_method,'SMM')
fprintf('Simulated method of moments with');
......@@ -49,7 +52,16 @@ if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_meth
fprintf('\n - penalized estimation using deviation from prior mean and weighted with prior precision');
end
end
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
if do_bayesian_estimation
fprintf('Bayesian Impulse Response Function Matching with');
else
fprintf('Frequentist Impulse Response Function Matching with');
end
if ~isempty(options_mom_.mom.irf_matching_file.name)
fprintf('\n - irf_matching_file: %s',[options_mom_.mom.irf_matching_file.path filesep options_mom_.mom.irf_matching_file.name '.m']);
end
end
for i = 1:length(options_mom_.optimizer_vec)
if i == 1
str = '- optimizer (mode_compute';
......@@ -118,6 +130,8 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
fprintf('\n - standard errors: numerical derivatives');
end
fprintf('\n - number of matched moments: %d', options_mom_.mom.mom_nbr);
elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
fprintf('\n - number of matched IRFs: %d', options_mom_.mom.mom_nbr);
end
fprintf('\n - number of parameters: %d', number_of_estimated_parameters);
fprintf('\n\n');
\ No newline at end of file
matlab/+mom/run.m
View file @ 4a32acc2
......@@ -6,31 +6,47 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
% o Preparing local options_mom_ structure
% o Checking the options and the compatibility of the settings
% o Initializations of variables, orderings and state space representation
% o Checks and transformations for matched moments structure
% o Checks and transformations for matched_moments structure
% o Checks and transformations for matched_irfs and matched_irfs_weights structure
% o Checks and transformations for estimated parameters, priors, and bounds
% o Checks and transformations for data
% o Checks for objective function at initial parameters
% o GMM/SMM: iterated method of moments estimation
% o GMM/SMM: J-Test and fit of moments%
% o Mode computation: optimization
% - GMM/SMM: iterated optimization
% - IRF_MATCHING: optimization
% o Bayesian MCMC estimation
% o Display of results
% - GMM/SMM: J-Test and fit of moments
% - IRF_MATCHING: fit of IRFs
% o Clean up
% -------------------------------------------------------------------------
% Note that we call a "mode" the minimum of the objective function, i.e.
% the parameter vector that minimizes the distance between the moments/IRFs
% computed from the model and the moments/IRFs computed from the data.
% -------------------------------------------------------------------------
% This function is inspired by replication codes accompanied to the following papers:
% GMM/SMM:
% o Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49.
% o Born, Pfeifer (2014): "Risk Matters: Comment", American Economic Review, 104(12):4231-4239.
% o Mutschler (2018): "Higher-order statistics for DSGE models", Econometrics and Statistics, 6:44-56.
% =========================================================================
% o Ruge-Murcia (2007): "Methods to Estimate Dynamic Stochastic General Equilibrium Models", Journal of Economic Dynamics and Control, 31(8):2599-2636.
% IRF MATCHING:
% o Christiano, Trabandt, Walentin (2010): "DSGE Models for Monetary Policy Analysis." In Handbook of Monetary Economics, 3:285–367.
% o Christiano, Eichenbaum, Trabandt (2016): "Unemployment and Business Cycles." Econometrica, 84: 1523-1569.
% o Ruge-Murcia (2020): "Estimating Nonlinear Dynamic Equilibrium Models by Matching Impulse Responses", Economics Letters, 197.
% -------------------------------------------------------------------------
% INPUTS
% o bayestopt_: [structure] information about priors
% o options_: [structure] information about global options
% o oo_: [structure] storage for results
% o oo_: [structure] results
% o estim_params_: [structure] information about estimated parameters
% o M_: [structure] information about model with
% o matched_moments: [cell] information about selected moments to match in GMM/SMM estimation
% vars: matched_moments{:,1});
% lead/lags: matched_moments{:,2};
% powers: matched_moments{:,3};
% o matched_irfs: [cell] information about selected IRFs to match in IRF_MATCHING estimation
% o matched_irfs_weights: [cell] information about entries in weight matrix for an IRF_MATCHING estimation
% o options_mom_: [structure] information about settings specified by the user
% -------------------------------------------------------------------------
% OUTPUTS
......@@ -42,25 +58,26 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
% o driver.m
% -------------------------------------------------------------------------
% This function calls
% o cellofchararraymaxlength
% o check_for_calibrated_covariances
% o check_mode_file
% o check_posterior_sampler_options
% o check_prior_bounds
% o check_prior_stderr_corr
% o check_steady_state_changes_parameters
% o check_varobs_are_endo_and_declared_once
% o check_hessian_at_the_mode
% o display_estimation_results_table
% o do_parameter_initialization
% o dyn_latex_table
% o dynare_minimize_objective
% o dyntable
% o get_all_parameters
% o get_dynare_random_generator_state
% o get_matrix_entries_for_psd_check
% o M_.fname '_prior_restrictions'
% o makedataset
% o mom.check_plot
% o mode_check
% o mom.check_irf_matching_file
% o mom.default_option_mom_values
% o mom.get_data_moments
% o mom.matched_irfs_blocks
% o mom.matched_moments_block
% o mom.objective_function
% o mom.optimal_weighting_matrix
......@@ -77,8 +94,8 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
% o set_state_space
% o skipline
% o test_for_deep_parameters_calibration
% o warning_config
% =========================================================================
% o transform_prior_to_laplace_prior
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
......@@ -95,34 +112,8 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
% Maintaining Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (johannes.pfeifer@unibw.de)
% =========================================================================
% -------------------------------------------------------------------------
% TO DO LISTS
% -------------------------------------------------------------------------
% GENERAL
% - document all options in manual
% - document analytic_jacobian better
% - make endogenous_prior_restrictions work
% - dirname option to save output to different directory not yet implemented
% - create test for prior restrictions file
% - add mode_file option
% - implement penalty objective
% - test optimizers
% GMM/SMM
% - speed up pruned_state_space_system (by using doubling with old initial values, hardcoding zeros, other "tricks" used in e.g. nlma)
% - add option to use autocorrelations (we have useautocorr in identification toolbox already)
% - SMM with extended path
% - deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
% - display scaled moments
% - enable first moments despite prefilter
% - do "true" Bayesian GMM and SMM not only penalized
fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom_method)
fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom_method);
% -------------------------------------------------------------------------
......@@ -130,29 +121,33 @@ fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom
% -------------------------------------------------------------------------
if isempty(estim_params_) % structure storing the info about estimated parameters in the estimated_params block
if ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
error('method_of_moments: You need to provide an ''estimated_params'' block!')
error('method_of_moments: You need to provide an ''estimated_params'' block!');
else
error('method_of_moments: The ''estimated_params'' block must not be empty!')
error('method_of_moments: The ''estimated_params'' block must not be empty!');
end
end
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
if ~isfield(M_,'matched_moments') || isempty(M_.matched_moments) % structure storing the moments used for GMM and SMM estimation
error('method_of_moments: You need to provide a ''matched_moments'' block for ''mom_method=%s''!',options_mom_.mom.mom_method)
error('method_of_moments: You need to provide a ''matched_moments'' block for ''mom_method=%s''!',options_mom_.mom.mom_method);
end
elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
if ~isfield(M_,'matched_irfs') || isempty(M_.matched_irfs) % structure storing the irfs used for matching
error('method_of_moments: You need to provide a ''matched_irfs'' block for ''mom_method=%s''!',options_mom_.mom.mom_method);
end
end
if (~isempty(estim_params_.var_endo) || ~isempty(estim_params_.corrn)) && strcmp(options_mom_.mom.mom_method, 'GMM')
error('method_of_moments: GMM estimation does not support measurement error(s) yet. Please specify them as a structural shock!')
error('method_of_moments: GMM estimation does not support measurement error(s) yet. Please specify them as a structural shock!');
end
doBayesianEstimation = [estim_params_.var_exo(:,5); estim_params_.var_endo(:,5); estim_params_.corrx(:,6); estim_params_.corrn(:,6); estim_params_.param_vals(:,5)];
if all(doBayesianEstimation~=0)
doBayesianEstimation = true;
elseif all(doBayesianEstimation==0)
doBayesianEstimation = false;
do_bayesian_estimation = [estim_params_.var_exo(:,5); estim_params_.var_endo(:,5); estim_params_.corrx(:,6); estim_params_.corrn(:,6); estim_params_.param_vals(:,5)];
if all(do_bayesian_estimation~=0)
do_bayesian_estimation = true;
elseif all(do_bayesian_estimation==0)
do_bayesian_estimation = false;
else
error('method_of_moments: Estimation must be either fully Frequentist or fully Bayesian. Maybe you forgot to specify a prior distribution!')
error('method_of_moments: Estimation must be either fully Frequentist or fully Bayesian. Maybe you forgot to specify a prior distribution!');
end
if ~isfield(options_,'varobs')
error('method_of_moments: VAROBS statement is missing!')
error('method_of_moments: VAROBS statement is missing!');
end
check_varobs_are_endo_and_declared_once(options_.varobs,M_.endo_names);
......@@ -166,28 +161,12 @@ check_varobs_are_endo_and_declared_once(options_.varobs,M_.endo_names);
% The idea is to be independent of options_ and have full control of the
% estimation instead of possibly having to deal with options chosen somewhere
% else in the mod file.
options_mom_ = mom.default_option_mom_values(options_mom_, options_, M_.dname, doBayesianEstimation);
options_mom_ = mom.default_option_mom_values(options_mom_, options_, M_.dname, M_.fname, do_bayesian_estimation);
% -------------------------------------------------------------------------
% workarounds
% -------------------------------------------------------------------------
% The TeX option crashes MATLAB R2014a run with "-nodisplay" option
% (as is done from the testsuite).
% Since we can’t directly test whether "-nodisplay" has been passed,
% we test for the "TOP_TEST_DIR" environment variable, which is set
% by the testsuite.
% Note that it was not tested whether the crash happens with more
% recent MATLAB versions, so when OLD_MATLAB_VERSION is increased,
% one should make a test before removing this workaround.
if options_.TeX && ~isoctave && matlab_ver_less_than('8.4') && ~isempty(getenv('TOP_TEST_DIR'))
warning('Disabling TeX option due to a bug in MATLAB R2014a with -nodisplay')
options_.TeX = false;
end
if isfield(options_mom_, 'TeX') && options_mom_.TeX && ~isoctave && matlab_ver_less_than('8.4') && ~isempty(getenv('TOP_TEST_DIR'))
warning('Disabling TeX option due to a bug in MATLAB R2014a with -nodisplay')
options_mom_.TeX = false;
end
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
% temporary workaround for https://git.dynare.org/Dynare/dseries/-/issues/51
if options_mom_.xls_sheet~=1
......@@ -199,60 +178,29 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
end
% -------------------------------------------------------------------------
% checks on settings
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
if numel(options_mom_.nobs) > 1
error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''nobs''!');
end
if numel(options_mom_.first_obs) > 1
error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''first_obs''!');
end
end
if options_mom_.order < 1
error('method_of_moments: The order of the Taylor approximation cannot be 0!')
end
if options_mom_.order > 2
fprintf('Dynare will use ''k_order_solver'' as the order>2\n');
options_mom_.k_order_solver = true;
end
if strcmp(options_mom_.mom.mom_method,'SMM')
if options_mom_.mom.simulation_multiple < 1
fprintf('The simulation horizon is shorter than the data. Dynare resets the simulation_multiple to 7.\n')
options_mom_.mom.simulation_multiple = 7;
end
end
if strcmp(options_mom_.mom.mom_method,'GMM')
% require pruning with GMM at higher order
if options_mom_.order > 1 && ~options_mom_.pruning
fprintf('GMM at higher order only works with pruning, so we set pruning option to 1.\n');
options_mom_.pruning = true;
end
if options_mom_.order > 3
error('method_of_moments: Perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM!');
end
end
if options_mom_.mom.analytic_jacobian && ~strcmp(options_mom_.mom.mom_method,'GMM')
options_mom_.mom.analytic_jacobian = false;
fprintf('\n''analytic_jacobian'' option will be dismissed as it only works with GMM.\n');
end
% -------------------------------------------------------------------------
% initializations
% -------------------------------------------------------------------------
%save warning state for restoring later on
orig_warning_state = warning;
% create output directories to store results
M_.dname = options_mom_.dirname;
CheckPath(M_.dname,'.');
CheckPath('method_of_moments',M_.dname);
CheckPath('graphs',options_mom_.dirname);
% initialize options that might change
options_mom_.mom.compute_derivs = false; % flag to compute derivs in objective function (might change for GMM with either analytic_standard_errors or analytic_jacobian (dependent on optimizer))
options_mom_.mom.vector_output = false; % specifies whether the objective function returns a vector
CheckPath('graphs',M_.dname);
if do_bayesian_estimation
oo_.mom.posterior.optimization.mode = [];
oo_.mom.posterior.optimization.Variance = [];
oo_.mom.posterior.optimization.log_density=[];
end
do_bayesian_estimation_mcmc = do_bayesian_estimation && ( (options_mom_.mh_replic>0) || options_mom_.load_mh_file );
invhess = [];
% decision rule
oo_.dr = set_state_space(oo_.dr,M_,options_mom_); % get state-space representation
oo_.mom.obs_var = []; % create index of observed variables in DR order
oo_.dr = set_state_space(oo_.dr,M_); % get state-space representation
options_mom_.mom.obs_var = []; % create index of observed variables in DR order
for i = 1:options_mom_.obs_nbr
oo_.mom.obs_var = [oo_.mom.obs_var; find(strcmp(options_mom_.varobs{i}, M_.endo_names(oo_.dr.order_var)))];
options_mom_.mom.obs_var = [options_mom_.mom.obs_var; find(strcmp(options_mom_.varobs{i}, M_.endo_names(oo_.dr.order_var)))];
end
......@@ -271,37 +219,76 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
% Get maximum lag number for autocovariances/autocorrelations
options_mom_.ar = max(cellfun(@max,M_.matched_moments(:,2))) - min(cellfun(@min,M_.matched_moments(:,2)));
% Check that only observed variables are involved in moments
not_observed_variables=setdiff(oo_.dr.inv_order_var([M_.matched_moments{:,1}]),oo_.mom.obs_var);
not_observed_variables=setdiff(oo_.dr.inv_order_var([M_.matched_moments{:,1}]),options_mom_.mom.obs_var);
if ~isempty(not_observed_variables)
skipline;
error('method_of_moments: You specified moments involving %s, but it is not a varobs!',M_.endo_names{oo_.dr.order_var(not_observed_variables)})
error('method_of_moments: You specified moments involving %s, but it is not a varobs!',M_.endo_names{oo_.dr.order_var(not_observed_variables)});
end
end
% -------------------------------------------------------------------------
% matched_irfs: checks and transformations
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
[oo_.mom.data_moments, oo_.mom.weighting_info.W, options_mom_.mom.irfIndex, options_mom_.irf] = mom.matched_irfs_blocks(M_.matched_irfs, M_.matched_irfs_weights, options_mom_.varobs_id, options_mom_.obs_nbr, M_.exo_nbr, M_.endo_names, M_.exo_names);
% compute inverse of weighting matrix
try
oo_.mom.weighting_info.Winv = inv(oo_.mom.weighting_info.W);
catch
error('method_of_moments: Something wrong while computing inv(W), check your weighting matrix!');
end
if any(isnan(oo_.mom.weighting_info.Winv(:))) || any(isinf(oo_.mom.weighting_info.Winv(:)))
error('method_of_moments: There are NaN or Inf values in inv(W), check your weighting matrix!');
end
% compute log determinant of inverse of weighting matrix in a robust way to avoid Inf or NaN
try
oo_.mom.weighting_info.Winv_logdet = 2*sum(log(diag(chol(oo_.mom.weighting_info.Winv))));
catch
error('method_of_moments: Something wrong while computing log(det(inv(W))), check your weighting matrix!');
end
if any(isnan(oo_.mom.weighting_info.Winv_logdet(:))) || any(isinf(oo_.mom.weighting_info.Winv_logdet(:)))
error('method_of_moments: There are NaN or Inf values in log(det(inv(W))), check your weighting matrix!');
end
options_mom_.mom.mom_nbr = length(options_mom_.mom.irfIndex);
end
% -------------------------------------------------------------------------
% irf_matching_file: checks and transformations
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
[options_mom_.mom.irf_matching_file.name, options_mom_.mom.irf_matching_file.path] = mom.check_irf_matching_file(options_mom_.mom.irf_matching_file.name);
% check for irf_matching_file
if ~( isempty(options_mom_.mom.irf_matching_file.path) || strcmp(options_mom_.mom.irf_matching_file.path,'.') )
fprintf('\nAdding %s to MATLAB''s path.\n',options_mom_.mom.irf_matching_file.path);
addpath(options_mom_.mom.irf_matching_file.path);
end
end
% -------------------------------------------------------------------------
% estimated parameters: checks and transformations on values, priors, bounds
% estimated parameters: checks and transformations on values, priors, bounds, posterior options
% -------------------------------------------------------------------------
% set priors and bounds over the estimated parameters
[xparam0, estim_params_, bayestopt_, lb, ub, M_] = set_prior(estim_params_, M_, options_mom_);
number_of_estimated_parameters = length(xparam0);
hessian_xparam0 = []; % initialize hessian
hessian_xparam0 = []; % initialize Hessian
% check if enough moments for estimation
if strcmp(options_mom_.mom.mom_method, 'GMM') || strcmp(options_mom_.mom.mom_method, 'SMM')
if options_mom_.mom.mom_nbr < length(xparam0)
skipline;
error('method_of_moments: There must be at least as many moments as parameters for a %s estimation!',options_mom_.mom.mom_method);
end
skipline(2);
end
% check if a _prior_restrictions.m file exists
if exist([M_.fname '_prior_restrictions.m'],'file')
options_mom_.prior_restrictions.status = 1;
options_mom_.prior_restrictions.routine = str2func([M_.fname '_prior_restrictions']);
end
% check that the provided mode_file is compatible with the current estimation settings
if ~isempty(options_mom_.mode_file) && ( ~do_bayesian_estimation || (do_bayesian_estimation && ~options_mom_.mh_posterior_mode_estimation) )
[xparam0, hessian_xparam0] = check_mode_file(xparam0, hessian_xparam0, options_mom_, bayestopt_);
end
% check on specified priors and penalized estimation (which uses Laplace approximated priors)
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
bayestopt_orig = bayestopt_;
......@@ -314,9 +301,8 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
bayestopt_ = mom.transform_prior_to_laplace_prior(bayestopt_);
end
end
% check for calibrated covariances before updating parameters
estim_params_ = check_for_calibrated_covariances(xparam0,estim_params_,M_);
estim_params_ = check_for_calibrated_covariances(estim_params_,M_);
% checks on parameter calibration and initialization
xparam_calib = get_all_parameters(estim_params_,M_); % get calibrated parameters
......@@ -326,89 +312,104 @@ else
estim_params_.full_calibration_detected = false;
end
if options_mom_.use_calibration_initialization % set calibration as starting values
if ~isempty(bayestopt_) && ~doBayesianEstimation && any(all(isnan([xparam_calib xparam0]),2))
if ~isempty(bayestopt_) && ~do_bayesian_estimation && any(all(isnan([xparam_calib xparam0]),2))
error('method_of_moments: When using the use_calibration option with %s without prior, the parameters must be explicitly initialized!',options_mom_.mom.mom_method);
else
[xparam0,estim_params_] = do_parameter_initialization(estim_params_,xparam_calib,xparam0); % get explicitly initialized parameters that have precedence over calibrated values
end
end
% check initialization
if ~isempty(bayestopt_) && ~doBayesianEstimation && any(isnan(xparam0))
if ~isempty(bayestopt_) && ~do_bayesian_estimation && any(isnan(xparam0))
error('method_of_moments: Frequentist %s requires all estimated parameters to be initialized, either in an estimated_params or estimated_params_init-block!',options_mom_.mom.mom_method);
end
% set and check parameter bounds
if ~isempty(bayestopt_) && doBayesianEstimation
if ~isempty(bayestopt_) && do_bayesian_estimation
% plot prior densities
if ~options_mom_.nograph && options_mom_.plot_priors
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
plot_priors(bayestopt_orig,M_,estim_params_,options_mom_,'Original priors'); % only for visual inspection (not saved to disk, because overwritten in next call to plot_priors)
plot_priors(bayestopt_,M_,estim_params_,options_mom_,'Laplace approximated priors');
clear('bayestopt_orig'); % make sure stale structure cannot be used
else
plot_priors(bayestopt_,M_,estim_params_,options_mom_,'Priors');
end
end
% set prior bounds
Bounds = prior_bounds(bayestopt_, options_mom_.prior_trunc);
Bounds.lb = max(Bounds.lb,lb);
Bounds.ub = min(Bounds.ub,ub);
BoundsInfo = prior_bounds(bayestopt_, options_mom_.prior_trunc);
BoundsInfo.lb = max(BoundsInfo.lb,lb);
BoundsInfo.ub = min(BoundsInfo.ub,ub);
else
% no priors are declared so Dynare will estimate the parameters with
% Frequentist methods using inequality constraints for the parameters
Bounds.lb = lb;
Bounds.ub = ub;
if options_mom_.mom.penalized_estimator
fprintf('Penalized estimation turned off as you did not declare priors\n')
% no priors are declared so Dynare will estimate the parameters with Frequentist methods using inequality constraints for the parameters
BoundsInfo.lb = lb;
BoundsInfo.ub = ub;
if (strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) && options_mom_.mom.penalized_estimator
fprintf('Penalized estimation turned off as you did not declare priors\n');
options_mom_.mom.penalized_estimator = 0;
else
if isfield(options_mom_,'mh_replic') && options_mom_.mh_replic > 0
fprintf('Setting ''mh_replic=0'' as you did not declare priors.\n');
options_mom_.mh_replic = 0;
end
end
end
% set correct bounds for standard deviations and correlations
Bounds = mom.set_correct_bounds_for_stderr_corr(estim_params_,Bounds);
BoundsInfo = mom.set_correct_bounds_for_stderr_corr(estim_params_,BoundsInfo);
% test if initial values of the estimated parameters are all between the prior lower and upper bounds
if options_mom_.use_calibration_initialization
try
check_prior_bounds(xparam0,Bounds,M_,estim_params_,options_mom_,bayestopt_);
check_prior_bounds(xparam0,BoundsInfo,M_,estim_params_,options_mom_,bayestopt_);
catch last_error
fprintf('Cannot use parameter values from calibration as they violate the prior bounds.')
fprintf('Cannot use parameter values from calibration as they violate the prior bounds.');
rethrow(last_error);
end
else
check_prior_bounds(xparam0,Bounds,M_,estim_params_,options_mom_,bayestopt_);
check_prior_bounds(xparam0,BoundsInfo,M_,estim_params_,options_mom_,bayestopt_);
end
% check for positive definiteness
estim_params_ = get_matrix_entries_for_psd_check(M_,estim_params_);
% set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file)
M_.sigma_e_is_diagonal = true;
if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances')
M_.sigma_e_is_diagonal = false;
end
% storing prior parameters in results
% storing prior parameters in results structure
if do_bayesian_estimation || ( (strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) && options_mom_.mom.penalized_estimator)
oo_.mom.prior.mean = bayestopt_.p1;
oo_.mom.prior.mode = bayestopt_.p5;
oo_.mom.prior.variance = diag(bayestopt_.p2.^2);
oo_.mom.prior.hyperparameters.first = bayestopt_.p6;
oo_.mom.prior.hyperparameters.second = bayestopt_.p7;
end
% set all parameters
M_ = set_all_parameters(xparam0,estim_params_,M_);
% provide warning if there is NaN in parameters
test_for_deep_parameters_calibration(M_);
if doBayesianEstimation
% set jscale
if do_bayesian_estimation_mcmc
if ~strcmp(options_mom_.posterior_sampler_options.posterior_sampling_method,'slice')
if isempty(options_mom_.mh_jscale)
options_mom_.mh_jscale = 2.38/sqrt(number_of_estimated_parameters); % use optimal value for univariate normal distribution (check_posterior_sampler_options and mode_compute=6 may overwrite this setting)
end
bayestopt_.jscale(find(isnan(bayestopt_.jscale))) = options_mom_.mh_jscale;
end
end
% initialization of posterior sampler options
if do_bayesian_estimation_mcmc
[current_options, options_mom_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_mom_, BoundsInfo, bayestopt_);
options_mom_.posterior_sampler_options.current_options = current_options;
if strcmp(current_options.posterior_sampling_method,'slice') && current_options.use_mh_covariance_matrix && ~current_options.rotated
error('method_of_moments: Using the slice sampler with the ''use_mh_covariance_matrix'' option requires also setting the ''rotated'' option!');
end
end
% warning if prior allows that stderr parameters are negative or corr parameters are outside the unit circle
if do_bayesian_estimation
check_prior_stderr_corr(estim_params_,bayestopt_);
% check value of prior density
[~,~,~,info] = priordens(xparam0,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
if any(info)
fprintf('The prior density evaluated at the initial values is Inf for the following parameters: %s\n',bayestopt_.name{info,1})
error('The initial value of the prior is -Inf!')
fprintf('The prior density evaluated at the initial values is Inf for the following parameters: %s\n',bayestopt_.name{info,1});
error('The initial value of the prior is -Inf!');
end
end
......@@ -416,34 +417,34 @@ end
% -------------------------------------------------------------------------
% datafile: checks and transformations
% -------------------------------------------------------------------------
% Build dataset
% build dataset
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
% Check if datafile has same name as mod file
[~,name,~] = fileparts(options_mom_.datafile);
% check if datafile has same name as mod file
[~,name] = fileparts(options_mom_.datafile);
if strcmp(name,M_.fname)
error('method_of_moments: ''datafile'' and mod file are not allowed to have the same name; change the name of the ''datafile''!')
error('method_of_moments: ''datafile'' and mod file are not allowed to have the same name; change the name of the ''datafile''!');
end
dataset_ = makedataset(options_mom_);
% set options for old interface from the ones for new interface
if ~isempty(dataset_)
options_mom_.nobs = dataset_.nobs;
end
% Check length of data for estimation of second moments
% check length of data for estimation of second moments
if options_mom_.ar > options_mom_.nobs+1
error('method_of_moments: Dataset is too short to compute higher than first moments!');
end
% Provide info on data moments handling
% provide info on data moments handling
fprintf('Computing data moments. Note that NaN values in the moments (due to leads and lags or missing data) are replaced by the mean of the corresponding moment.\n');
% Get data moments for the method of moments
[oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_);
% get data moments for the method of moments
[oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, options_mom_.mom.obs_var, oo_.dr.inv_order_var, M_.matched_moments, options_mom_);
if ~isreal(dataset_.data)
error('method_of_moments: The data moments contain complex values!')
error('method_of_moments: The data moments contain complex values!');
end
end
% -------------------------------------------------------------------------
% SMM: Get shock series fand set variance correction factor
% SMM: Get shock series and set variance correction factor
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'SMM')
options_mom_.mom.long = round(options_mom_.mom.simulation_multiple*options_mom_.nobs);
......@@ -477,7 +478,6 @@ end
% -------------------------------------------------------------------------
% checks for steady state at initial parameters
% -------------------------------------------------------------------------
% check if steady state solves static model and if steady-state changes estimated parameters
if options_mom_.steadystate.nocheck
steadystate_check_flag_vec = [0 1];
......@@ -486,7 +486,7 @@ else
end
[oo_.steady_state, info, steady_state_changes_parameters] = check_steady_state_changes_parameters(M_, estim_params_, oo_, options_mom_, steadystate_check_flag_vec);
if info(1)
fprintf('\nThe steady state at the initial parameters cannot be computed.\n')
fprintf('\nThe steady state at the initial parameters cannot be computed.\n');
print_info(info, 0, options_mom_);
end
if steady_state_changes_parameters && strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
......@@ -494,7 +494,6 @@ if steady_state_changes_parameters && strcmp(options_mom_.mom.mom_method,'GMM')
fprintf('because the steady-state changes estimated parameters. Option ''analytic_derivation_mode'' reset to -2.');
options_mom_.analytic_derivation_mode = -2;
end
% display warning if some parameters are still NaN
test_for_deep_parameters_calibration(M_);
......@@ -504,17 +503,25 @@ test_for_deep_parameters_calibration(M_);
% -------------------------------------------------------------------------
objective_function = str2func('mom.objective_function');
try
% Check for NaN or complex values of moment-distance-funtion evaluated at initial parameters
% check for NaN or complex values of moment-distance-function evaluated at initial parameters
if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
oo_.mom.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix
oo_.mom.weighting_info.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix
end
tic_id = tic;
[fval, info, ~, ~, ~, oo_, M_] = feval(objective_function, xparam0, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_);
[fval, info] = feval(objective_function, xparam0, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
elapsed_time = toc(tic_id);
if isnan(fval)
error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!')
if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!');
else
error('method_of_moments: The initial value of the objective function is NaN!');
end
elseif imag(fval)
error('method_of_moments: The initial value of the objective function with identity weighting matrix is complex!')
if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
error('method_of_moments: The initial value of the objective function with identity weighting matrix is complex!');
else
error('method_of_moments: The initial value of the objective function is complex!');
end
end
if info(1) > 0
disp('method_of_moments: Error in computing the objective function for initial parameter values')
......@@ -529,10 +536,10 @@ try
catch last_error % if check fails, provide info on using calibration if present
if estim_params_.full_calibration_detected %calibrated model present and no explicit starting values
skipline(1);
fprintf('There was an error in computing the moments for initial parameter values.\n')
fprintf('If this is not a problem with the setting of options (check the error message below),\n')
fprintf('you should try using the calibrated version of the model as starting values. To do\n')
fprintf('this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n')
fprintf('There was an error in computing the moments for initial parameter values.\n');
fprintf('If this is not a problem with the setting of options (check the error message below),\n');
fprintf('you should try using the calibrated version of the model as starting values. To do\n');
fprintf('this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n');
fprintf('command (and after the estimated_params-block so that it does not get overwritten):\n');
skipline(2);
end
......@@ -543,60 +550,282 @@ end
% -------------------------------------------------------------------------
% print some info to console
% -------------------------------------------------------------------------
mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters);
mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, do_bayesian_estimation);
% -------------------------------------------------------------------------
% compute mode for GMM/SMM
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
[xparam1, oo_.mom.weighting_info, oo_.mom.verbose] = mom.mode_compute_gmm_smm(xparam0, objective_function, oo_.mom.m_data, oo_.mom.data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
end
% -------------------------------------------------------------------------
% compute mode for IRF matching
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
if ~do_bayesian_estimation || (do_bayesian_estimation && ~options_mom_.mh_posterior_mode_estimation)
[xparam1, hessian_xparam1, fval, oo_.mom.verbose] = mom.mode_compute_irf_matching(xparam0, hessian_xparam0, objective_function, do_bayesian_estimation, oo_.mom.weighting_info, oo_.mom.data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
else
xparam1 = xparam0;
hessian_xparam1 = hessian_xparam0;
end
end
% -------------------------------------------------------------------------
% GMM/SMM: iterated estimation
% compute standard errors and initialize covariance of the proposal distribution
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
% compute mode
[xparam1, oo_, Woptflag] = mom.mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds);
% compute standard errors at mode
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
options_mom_.mom.compute_derivs = true; % for GMM we compute derivatives analytically in the objective function with this flag
end
[~, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_); % compute model moments and oo_.mom.model_moments_params_derivs
[~, ~, ~, ~, ~, oo_.mom.Q, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs] = feval(objective_function, xparam1, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization
[stdh,hessian_xparam1] = mom.standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag);
[stdh, invhess] = mom.standard_errors(xparam1, objective_function, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs, oo_.mom.m_data, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if options_mom_.cova_compute
hessian_xparam1 = inv(invhess);
end
else
if ~do_bayesian_estimation || ~options_mom_.mh_posterior_mode_estimation
if do_bayesian_estimation
oo_.mom.posterior.optimization.mode = xparam1;
if exist('fval','var')
oo_.mom.posterior.optimization.log_density = -fval;
end
end
if options_mom_.cova_compute
hsd = sqrt(diag(hessian_xparam1)); % represent the curvature (or second derivatives) of the likelihood with respect to each parameter being estimated.
invhess = inv(hessian_xparam1./(hsd*hsd'))./(hsd*hsd'); % before taking the inverse scale the Hessian matrix by dividing each of its elements by the outer product of hsd such that the diagonal of the resulting matrix is approximately 1. This kind of scaling can help in regularizing the matrix and potentially improves its condition number, which in turn can make the matrix inversion more stable.
stdh = sqrt(diag(invhess));
if do_bayesian_estimation
oo_.mom.posterior.optimization.Variance = invhess;
end
end
else
variances = bayestopt_.p2.*bayestopt_.p2;
id_Inf = isinf(variances);
variances(id_Inf) = 1;
invhess = options_mom_.mh_posterior_mode_estimation*diag(variances);
xparam1 = bayestopt_.p5;
id_NaN = isnan(xparam1);
xparam1(id_NaN) = bayestopt_.p1(id_NaN);
outside_bound_pars=find(xparam1 < BoundsInfo.lb | xparam1 > BoundsInfo.ub);
xparam1(outside_bound_pars) = bayestopt_.p1(outside_bound_pars);
end
if ~options_mom_.cova_compute
stdh = NaN(length(xparam1),1);
end
end
% -------------------------------------------------------------------------
% display estimation results at mode
% -------------------------------------------------------------------------
if do_bayesian_estimation && ~options_mom_.mom.penalized_estimator && ~options_mom_.mh_posterior_mode_estimation
% display table with Bayesian mode estimation results and store parameter estimates and standard errors in oo_
oo_.mom = display_estimation_results_table(xparam1, stdh, M_, options_mom_, estim_params_, bayestopt_, oo_.mom, prior_dist_names, 'Posterior', 'posterior');
% Laplace approximation to the marginal log density
if options_mom_.cova_compute
estim_params_nbr = size(xparam1,1);
if ispd(invhess)
log_det_invhess = log(det(invhess./(stdh*stdh')))+2*sum(log(stdh));
likelihood = feval(objective_function, xparam1, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
oo_.mom.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
else
oo_.mom.MarginalDensity.LaplaceApproximation = NaN;
end
fprintf('\nLog data density [Laplace approximation] is %f.\n',oo_.mom.MarginalDensity.LaplaceApproximation);
end
elseif ~do_bayesian_estimation || (do_bayesian_estimation && options_mom_.mom.penalized_estimator)
% display table with Frequentist estimation results and store parameter estimates and standard errors in oo_
oo_.mom = display_estimation_results_table(xparam1, stdh, M_, options_mom_, estim_params_, bayestopt_, oo_.mom, prior_dist_names, options_mom_.mom.mom_method, lower(options_mom_.mom.mom_method));
end
% -------------------------------------------------------------------------
% checks for mode and hessian at the mode
% -------------------------------------------------------------------------
if (~do_bayesian_estimation && options_mom_.cova_compute) || (do_bayesian_estimation && ~options_mom_.mh_posterior_mode_estimation && options_mom_.cova_compute)
check_hessian_at_the_mode(hessian_xparam1, xparam1, M_, estim_params_, options_, BoundsInfo);
end
if options_mom_.mode_check.status
mode_check(objective_function, xparam1, hessian_xparam1, options_mom_, M_, estim_params_, bayestopt_, Bounds, true,...
Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_);
if ~do_bayesian_estimation || (do_bayesian_estimation && ~options_mom_.mh_posterior_mode_estimation)
mode_check(objective_function, xparam1, diag(stdh), options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, true,... % use diag(stdh) instead of hessian_xparam1 as mode_check uses diagonal elements
oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
end
end
% -------------------------------------------------------------------------
% Bayesian MCMC estimation
% -------------------------------------------------------------------------
if do_bayesian_estimation_mcmc
invhess = set_mcmc_jumping_covariance(invhess, length(xparam1), options_mom_.MCMC_jumping_covariance, bayestopt_, 'method_of_moments');
% reset bounds as lb and ub must only be operational during mode-finding
BoundsInfo = set_mcmc_prior_bounds(xparam1, bayestopt_, options_mom_, 'method_of_moments');
% tunes the jumping distribution's scale parameter
if isfield(options_mom_,'mh_tune_jscale') && options_mom_.mh_tune_jscale.status
if strcmp(options_mom_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings')
options_mom_.mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_mom_, M_, objective_function, xparam1, BoundsInfo,...
oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
bayestopt_.jscale(:) = options_mom_.mh_jscale;
fprintf('mh_tune_jscale: mh_jscale has been set equal to %s.\n', num2str(options_mom_.mh_jscale));
else
warning('mh_tune_jscale is only available with ''random_walk_metropolis_hastings''!')
end
end
% run MCMC sampling
posterior_sampler_options = options_mom_.posterior_sampler_options.current_options;
posterior_sampler_options.invhess = invhess;
[posterior_sampler_options, options_mom_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, M_.fname, M_.dname, options_mom_, BoundsInfo, bayestopt_,'method_of_moments');
options_mom_.posterior_sampler_options.current_options = posterior_sampler_options; % store current options
if options_mom_.mh_replic>0
posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,BoundsInfo,oo_.mom.data_moments,oo_.mom.weighting_info,options_mom_,M_,estim_params_,bayestopt_,oo_,'method_of_moments::mcmc');
end
CutSample(M_, options_mom_, 'method_of_moments::mcmc'); % discard first mh_drop percent of the draws
if options_mom_.mh_posterior_mode_estimation
% skip optimizer-based mode-finding and instead compute the mode based on a run of a MCMC
[~,~,posterior_mode,~] = compute_mh_covariance_matrix(bayestopt_,M_.fname,M_.dname,'method_of_moments');
oo_.mom = fill_mh_mode(posterior_mode',NaN(length(posterior_mode),1),M_,options_mom_,estim_params_,bayestopt_,oo_.mom,'posterior');
warning(orig_warning_state);
return
else
% get stored results if required
if options_mom_.load_mh_file && options_mom_.load_results_after_load_mh
oo_load_mh = load([M_.dname filesep 'method_of_moments' filesep M_.fname '_mom_results'],'oo_');
end
% convergence diagnostics
if ~options_mom_.nodiagnostic
if (options_mom_.mh_replic>0 || (options_mom_.load_mh_file && ~options_mom_.load_results_after_load_mh))
oo_.mom = mcmc_diagnostics(options_mom_, estim_params_, M_, oo_.mom);
elseif options_mom_.load_mh_file && options_mom_.load_results_after_load_mh
if isfield(oo_load_mh.oo_.mom,'convergence')
oo_.mom.convergence = oo_load_mh.oo_.mom.convergence;
end
end
end
% statistics and plots for posterior draws
if options_mom_.mh_replic || (options_mom_.load_mh_file && ~options_mom_.load_results_after_load_mh)
[~,oo_.mom] = marginal_density(M_, options_mom_, estim_params_, oo_.mom, bayestopt_, 'method_of_moments');
oo_.mom = GetPosteriorParametersStatistics(estim_params_, M_, options_mom_, bayestopt_, oo_.mom, prior_dist_names);
if ~options_mom_.nograph
oo_.mom = PlotPosteriorDistributions(estim_params_, M_, options_mom_, bayestopt_, oo_.mom);
end
[oo_.mom.posterior.metropolis.mean,oo_.mom.posterior.metropolis.Variance] = GetPosteriorMeanVariance(options_mom_,M_);
elseif options_mom_.load_mh_file && options_mom_.load_results_after_load_mh
% load fields from previous MCMC run stored in results-file
field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density
'posterior_mean','posterior_hpdinf','posterior_hpdsup','posterior_median','posterior_variance','posterior_std','posterior_deciles','posterior_density',...% fields set by GetPosteriorParametersStatistics
'prior_density',...% fields set by PlotPosteriorDistributions
};
for field_iter=1:size(field_names,2)
if isfield(oo_load_mh.oo_.mom,field_names{1,field_iter})
oo_.mom.(field_names{1,field_iter}) = oo_load_mh.oo_.mom.(field_names{1,field_iter});
end
end
if isfield(oo_load_mh.oo_.mom,'MarginalDensity') && isfield(oo_load_mh.oo_.mom.MarginalDensity,'ModifiedHarmonicMean') % field set by marginal_density
oo_.mom.MarginalDensity.ModifiedHarmonicMean = oo_load_mh.oo_.mom.MarginalDensity.ModifiedHarmonicMean;
end
if isfield(oo_load_mh.oo_.mom,'posterior') && isfield(oo_load_mh.oo_.mom.posterior,'metropolis') % field set by GetPosteriorMeanVariance
oo_.mom.posterior.metropolis = oo_load_mh.oo_.mom.posterior.metropolis;
end
end
[options_mom_.sub_draws, error_flag]=set_number_of_subdraws(M_,options_mom_); %check whether number of sub_draws is feasible
if ~(~isempty(options_mom_.sub_draws) && options_mom_.sub_draws==0)
% THIS IS PROBABLY NOT USEFUL HERE AND CAN BE REMOVED (PREPROCESSOR: REMOVE bayesian_irf, moments_varendo)
%if options_mom_.bayesian_irf
% if error_flag
% error('method_of_moments: Cannot compute the posterior IRFs!');
% end
% PosteriorIRF('posterior','method_of_moments::mcmc');
%end
% if options_mom_.moments_varendo
% if error_flag
% error('method_of_moments: Cannot compute the posterior moments for the endogenous variables!');
% end
% if options_mom_.load_mh_file && options_mom_.mh_replic==0 %user wants to recompute results
% [MetropolisFolder, info] = CheckPath('metropolis',options_mom_.dirname);
% if ~info
% generic_post_data_file_name={'Posterior2ndOrderMoments','decomposition','PosteriorVarianceDecomposition','correlation','PosteriorCorrelations','conditional decomposition','PosteriorConditionalVarianceDecomposition'};
% for ii=1:length(generic_post_data_file_name)
% delete_stale_file([MetropolisFolder filesep M_.fname '_' generic_post_data_file_name{1,ii} '*']);
% end
% % restore compatibility for loading pre-4.6.2 runs where estim_params_ was not saved; see 6e06acc7 and !1944
% NumberOfDrawsFiles = length(dir([M_.dname '/metropolis/' M_.fname '_posterior_draws*' ]));
% if NumberOfDrawsFiles>0
% temp=load([M_.dname '/metropolis/' M_.fname '_posterior_draws1']);
% if ~isfield(temp,'estim_params_')
% for file_iter=1:NumberOfDrawsFiles
% save([M_.dname '/metropolis/' M_.fname '_posterior_draws' num2str(file_iter)],'estim_params_','-append')
% end
% end
% end
% end
% end
% oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
% end
else
fprintf('''sub_draws'' was set to 0. Skipping posterior computations.');
end
xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_.mom,options_);
end
% MAYBE USEFUL????
% % Posterior correlations
% extreme_corr_bound = 0.7;
% if ~isnan(extreme_corr_bound)
% tril_para_correlation_matrix=tril(para_correlation_matrix,-1);
% [row_index,col_index]=find(abs(tril_para_correlation_matrix)>extreme_corr_bound);
% extreme_corr_params=cell(length(row_index),3);
% for i=1:length(row_index)
% extreme_corr_params{i,1}=char(parameter_names(row_index(i),:));
% extreme_corr_params{i,2}=char(parameter_names(col_index(i),:));
% extreme_corr_params{i,3}=tril_para_correlation_matrix(row_index(i),col_index(i));
% end
% end
% disp(' ');
% disp(['Correlations of Parameters (at Posterior Mode) > ',num2str(extreme_corr_bound)]);
% disp(extreme_corr_params)
end
% -------------------------------------------------------------------------
% display final estimation results
% -------------------------------------------------------------------------
M_ = set_all_parameters(xparam1,estim_params_,M_); % update parameters
[~, ~, ~, ~, ~, oo_.mom.Q, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs, oo_.mom.irf_model_varobs] = objective_function(xparam1, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); % store final results in oo_.mom
if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
% Store results in output structure
oo_.mom = display_estimation_results_table(xparam1,stdh,M_,options_mom_,estim_params_,bayestopt_,oo_.mom,prior_dist_names,options_mom_.mom.mom_method,lower(options_mom_.mom.mom_method));
% J test
oo_ = mom.Jtest(xparam1, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, dataset_.nobs);
% display comparison of model moments and data moments
mom.display_comparison_moments(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments);
oo_.mom.J_test = mom.Jtest(xparam1, objective_function, oo_.mom.Q, oo_.mom.model_moments, oo_.mom.m_data, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
if ~options_mom_.nograph
mom.graph_comparison_irfs(M_.matched_irfs,oo_.mom.irf_model_varobs,options_mom_.varobs_id,options_mom_.irf,options_mom_.relative_irf,M_.endo_names,M_.endo_names_tex,M_.exo_names,M_.exo_names_tex,M_.dname,M_.fname,options_mom_.graph_format,options_mom_.TeX,options_mom_.nodisplay,options_mom_.figures.textwidth)
end
end
% display comparison of model moments/IRFs and data moments/IRFs
mom.display_comparison_moments_irfs(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments);
% save results to _mom_results.mat
save([M_.dname filesep 'method_of_moments' filesep M_.fname '_mom_results.mat'], 'oo_', 'options_mom_', 'M_', 'estim_params_', 'bayestopt_');
fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mom_.mom.mom_method);
% -------------------------------------------------------------------------
% clean up
% -------------------------------------------------------------------------
fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mom_.mom.mom_method)
%reset warning state
warning_config;
warning(orig_warning_state); %reset warning state
if isoctave && isfield(options_mom_, 'prior_restrictions') && ...
isfield(options_mom_.prior_restrictions, 'routine')
% Octave crashes if it tries to save function handles (to the _results.mat file)
% See https://savannah.gnu.org/bugs/?43215
options_mom_.prior_restrictions.routine = [];
end
if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
if isfield(oo_.mom,'irf_model_varobs') && isempty(oo_.mom.irf_model_varobs)
oo_.mom = rmfield(oo_.mom,'irf_model_varobs'); % remove empty field
end
end
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && ~isempty(options_mom_.mom.irf_matching_file.path) && ~strcmp(options_mom_.mom.irf_matching_file.path,'.')
rmpath(options_mom_.irf_matching_file.path); % remove path to irf_matching_file
end
\ No newline at end of file
matlab/+mom/set_correct_bounds_for_stderr_corr.m
View file @ 4a32acc2
function Bounds = set_correct_bounds_for_stderr_corr(estim_params_,Bounds)
% function Bounds = set_correct_bounds_for_stderr_corr(estim_params_,Bounds)
function BoundsInfo = set_correct_bounds_for_stderr_corr(estim_params_,BoundsInfo)
% BoundsInfo = set_correct_bounds_for_stderr_corr(estim_params_,BoundsInfo)
% -------------------------------------------------------------------------
% Set correct bounds for standard deviation and corrrelation parameters
% =========================================================================
% -------------------------------------------------------------------------
% INPUTS
% o estim_params_ [struct] information on estimated parameters
% o Bounds [struct] information on bounds
% o BoundsInfo [struct] information on bounds
% -------------------------------------------------------------------------
% OUTPUT
% o Bounds [struct] updated bounds
% o BoundsInfo [struct] updated bounds
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
......@@ -29,15 +30,15 @@ function Bounds = set_correct_bounds_for_stderr_corr(estim_params_,Bounds)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
number_of_estimated_parameters = estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np;
% set correct bounds for standard deviations and corrrelations
param_of_interest = (1:number_of_estimated_parameters)'<=estim_params_.nvx+estim_params_.nvn;
LB_below_0 = (Bounds.lb<0 & param_of_interest);
Bounds.lb(LB_below_0) = 0;
LB_below_0 = (BoundsInfo.lb<0 & param_of_interest);
BoundsInfo.lb(LB_below_0) = 0;
param_of_interest = (1:number_of_estimated_parameters)'> estim_params_.nvx+estim_params_.nvn & (1:number_of_estimated_parameters)'<estim_params_.nvx+estim_params_.nvn +estim_params_.ncx + estim_params_.ncn;
LB_below_minus_1 = (Bounds.lb<-1 & param_of_interest);
UB_above_1 = (Bounds.ub>1 & param_of_interest);
Bounds.lb(LB_below_minus_1) = -1;
Bounds.ub(UB_above_1) = 1;
\ No newline at end of file
LB_below_minus_1 = (BoundsInfo.lb<-1 & param_of_interest);
UB_above_1 = (BoundsInfo.ub>1 & param_of_interest);
BoundsInfo.lb(LB_below_minus_1) = -1;
BoundsInfo.ub(UB_above_1) = 1;
\ No newline at end of file
matlab/+mom/standard_errors.m
View file @ 4a32acc2
function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Wopt_flag)
% [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Wopt_flag)
function [stderr_values, asympt_cov_mat] = standard_errors(xparam, objective_function, model_moments, model_moments_params_derivs, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [stderr_values, asympt_cov_mat] = standard_errors(xparam, objective_function, model_moments, model_moments_params_derivs, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% This function computes standard errors to the method of moments estimates
% Adapted from replication codes of
% o Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49.
% =========================================================================
% Adapted from replication codes of Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018):
% "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications",
% Review of Economic Studies, 85(1):1-49.
% -------------------------------------------------------------------------
% INPUTS
% o xparam: value of estimated parameters as returned by set_prior()
% o objective_function string of objective function
% o Bounds: structure containing parameter bounds
% o oo_: structure for results
% o estim_params_: structure describing the estimated_parameters
% o M_ structure describing the model
% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
% o Wopt_flag: indicator whether the optimal weighting is actually used
% - xparam: [vector] value of estimated parameters as returned by set_prior()
% - objective_function [func] function handle with string of objective function
% - model_moments: [vector] model moments
% - model_moments_params_derivs: [matrix] analytical Jacobian of the model moments wrt estimated parameters (currently for GMM only)
% - m_data [matrix] selected empirical moments at each point in time
% - data_moments: [vector] data with moments/IRFs to match
% - weighting_info: [structure] storing information on weighting matrices
% - options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% - M_ [structure] model information
% - estim_params_: [structure] information from estimated_params block
% - bayestopt_: [structure] information on the prior distributions
% - BoundsInfo: [structure] parameter bounds
% - dr: [structure] reduced form model
% - endo_steady_state: [vector] steady state value for endogenous variables (initval)
% - exo_steady_state: [vector] steady state value for exogenous variables (initval)
% - exo_det_steady_state: [vector] steady state value for exogenous deterministic variables (initval)
% -------------------------------------------------------------------------
% OUTPUTS
% o SE_values [nparam x 1] vector of standard errors
% o Asympt_Var [nparam x nparam] asymptotic covariance matrix
% o stderr_values [nparam x 1] vector of standard errors
% o asympt_cov_mat [nparam x nparam] asymptotic covariance matrix
% -------------------------------------------------------------------------
% This function is called by
% o mom.run.m
......@@ -27,8 +36,9 @@ function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, B
% o get_error_message
% o mom.objective_function
% o mom.optimal_weighting_matrix
% =========================================================================
% Copyright © 2020-2021 Dynare Team
% -------------------------------------------------------------------------
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -44,84 +54,74 @@ function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, B
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
% Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (johannes.pfeifer@unibw.de)
% =========================================================================
% Some dimensions
num_mom = size(oo_.mom.model_moments,1);
num_mom = size(model_moments,1);
dim_params = size(xparam,1);
D = zeros(num_mom,dim_params);
eps_value = options_mom_.mom.se_tolx;
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
fprintf('\nComputing standard errors using analytical derivatives of moments\n');
D = oo_.mom.model_moments_params_derivs; %already computed in objective function via get_perturbation_params.m
D = model_moments_params_derivs; % already computed in objective function via get_perturbation_params.m
idx_nan = find(any(isnan(D)));
if any(idx_nan)
for i = idx_nan
fprintf('No standard errors available for parameter %s\n',get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_))
fprintf('No standard errors available for parameter %s\n',get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_.varobs))
end
warning('There are NaN in the analytical Jacobian of Moments. Check your bounds and/or priors, or use a different optimizer.')
Asympt_Var = NaN(length(xparam),length(xparam));
SE_values = NaN(length(xparam),1);
asympt_cov_mat = NaN(length(xparam),length(xparam));
stderr_values = NaN(length(xparam),1);
return
end
else
fprintf('\nComputing standard errors using numerical derivatives of moments\n');
for i=1:dim_params
%Positive step
% positive step
xparam_eps_p = xparam;
xparam_eps_p(i,1) = xparam_eps_p(i) + eps_value;
[~, info_p, ~, ~,~, oo__p] = feval(objective_function, xparam_eps_p, Bounds, oo_, estim_params_, M_, options_mom_);
% Negative step
[~, info_p, ~, ~, ~, ~, model_moments_p] = feval(objective_function, xparam_eps_p, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
% negative step
xparam_eps_m = xparam;
xparam_eps_m(i,1) = xparam_eps_m(i) - eps_value;
[~, info_m, ~, ~,~, oo__m] = feval(objective_function, xparam_eps_m, Bounds, oo_, estim_params_, M_, options_mom_);
% The Jacobian:
[~, info_m, ~, ~, ~, ~, model_moments_m] = feval(objective_function, xparam_eps_m, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
% the Jacobian
if nnz(info_p)==0 && nnz(info_m)==0
D(:,i) = (oo__p.mom.model_moments - oo__m.mom.model_moments)/(2*eps_value);
D(:,i) = (model_moments_p - model_moments_m)/(2*eps_value);
else
problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_);
problematic_parameter = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_.varobs);
if info_p(1)==42
warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s due to hitting the upper bound - no standard errors available.\n',problpar)
warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s due to hitting the upper bound - no standard errors available.\n',problematic_parameter)
else
message_p = get_error_message(info_p, options_mom_);
end
if info_m(1)==41
warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s due to hitting the lower bound - no standard errors available.\n',problpar)
warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s due to hitting the lower bound - no standard errors available.\n',problematic_parameter)
else
message_m = get_error_message(info_m, options_mom_);
end
if info_m(1)~=41 && info_p(1)~=42
warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s - no standard errors available\n %s %s\nCheck your priors or use a different optimizer.\n',problpar, message_p, message_m)
warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s - no standard errors available\n %s %s\nCheck your priors or use a different optimizer.\n',problematic_parameter, message_p, message_m)
end
Asympt_Var = NaN(length(xparam),length(xparam));
SE_values = NaN(length(xparam),1);
asympt_cov_mat = NaN(length(xparam),length(xparam));
stderr_values = NaN(length(xparam),1);
return
end
end
end
T = options_mom_.nobs; %Number of observations
T = options_mom_.nobs;
if isfield(options_mom_,'variance_correction_factor')
T = T*options_mom_.variance_correction_factor;
end
WW = oo_.mom.Sw'*oo_.mom.Sw;
if Wopt_flag
% We have the optimal weighting matrix
Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params));
WW = weighting_info.Sw'*weighting_info.Sw;
if weighting_info.Woptflag
% we already have the optimal weighting matrix
asympt_cov_mat = 1/T*((D'*WW*D)\eye(dim_params));
else
% We do not have the optimal weighting matrix yet
WWopt = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
% we do not have the optimal weighting matrix yet
WWopt = mom.optimal_weighting_matrix(m_data, model_moments, options_mom_.mom.bartlett_kernel_lag);
S = WWopt\eye(size(WWopt,1));
AA = (D'*WW*D)\eye(dim_params);
Asympt_Var = 1/T*AA*D'*WW*S*WW*D*AA;
asympt_cov_mat = 1/T*AA*D'*WW*S*WW*D*AA;
end
SE_values = sqrt(diag(Asympt_Var));
stderr_values = sqrt(diag(asympt_cov_mat));
\ No newline at end of file
matlab/+mom/transform_prior_to_laplace_prior.m
View file @ 4a32acc2
function bayestopt_ = transform_prior_to_laplace_prior(bayestopt_)
% function bayestopt_ = transform_prior_to_laplace_prior(bayestopt_)
% bayestopt_ = transform_prior_to_laplace_prior(bayestopt_)
% -------------------------------------------------------------------------
% Transforms the prior specification to a Laplace type of approximation:
% only the prior mean and standard deviation are relevant, all other shape
% information, except for the parameter bounds, is ignored.
% =========================================================================
% -------------------------------------------------------------------------
% INPUTS
% bayestopt_ [structure] prior information
% -------------------------------------------------------------------------
......@@ -13,7 +13,8 @@ function bayestopt_ = transform_prior_to_laplace_prior(bayestopt_)
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
......@@ -30,7 +31,7 @@ function bayestopt_ = transform_prior_to_laplace_prior(bayestopt_)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
if any(setdiff([0;bayestopt_.pshape],[0,3]))
fprintf('\nNon-normal priors specified. Penalized estimation uses a Laplace type of approximation:');
fprintf('\nOnly the prior mean and standard deviation are relevant, all other shape information, except for the parameter bounds, is ignored.\n\n');
......
matlab/+occbin/DSGE_smoother.m
View file @ 4a32acc2
function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
%function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
% [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
% Runs a DSGE smoother with occasionally binding constraints
%
% INPUTS
......@@ -8,13 +8,13 @@ function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,P
% - Y [double] (n*T) matrix of data.
% - data_index [cell] 1*smpl cell of column vectors of indices.
% - missing_value [boolean] 1 if missing values, 0 otherwise
% - M_ [structure] Matlab's structure describing the model (M_).
% - oo_ [structure] Matlab's structure containing the results (oo_).
% - options_ [structure] Matlab's structure describing the current options (options_).
% - M_ [structure] MATLAB's structure describing the model (M_).
% - oo_ [structure] MATLAB's structure containing the results (oo_).
% - options_ [structure] MATLAB's structure describing the current options (options_).
% - bayestopt_ [structure] describing the priors
% - estim_params_ [structure] characterizing parameters to be estimated
% - dataset_ [structure] the dataset after required transformation
% - dataset_info [structure] Various informations about the dataset (descriptive statistics and missing observations)
% - dataset_info [structure] Various information about the dataset (descriptive statistics and missing observations)
%
% OUTPUTS
% - alphahat [double] (m*T) matrix, smoothed endogenous variables (a_{t|T}) (decision-rule order)
......@@ -35,7 +35,7 @@ function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,P
% - Trend [double] (n*T) pure trend component; stored in options_.varobs order
% - state_uncertainty [double] (K,K,T) array, storing the uncertainty
% about the smoothed state (decision-rule order)
% - M_ [structure] decribing the model
% - M_ [structure] describing the model
% - oo_ [structure] storing the results
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
......@@ -65,8 +65,10 @@ regime_history=[];
if options_.occbin.smoother.linear_smoother && nargin==12
%% linear smoother
options_.occbin.smoother.status=false;
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
tmp_smoother=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = ...
DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
tmp_smoother=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,...
aK,P,PK,decomp,Trend,state_uncertainty);
for jf=1:length(smoother_field_list)
oo_.occbin.linear_smoother.(smoother_field_list{jf}) = tmp_smoother.(smoother_field_list{jf});
end
......@@ -80,8 +82,10 @@ if options_.occbin.smoother.linear_smoother && nargin==12
oo_.occbin.linear_smoother.T0=T0;
oo_.occbin.linear_smoother.R0=R0;
oo_.occbin.linear_smoother.decomp=decomp;
oo_.occbin.linear_smoother.alphahat0=alphahat0;
oo_.occbin.linear_smoother.state_uncertainty0=state_uncertainty0;
fprintf('\nOccbin: linear smoother done.\n')
disp_verbose('OccBin: linear smoother done.',options_.verbosity)
options_.occbin.smoother.status=true;
end
% if init_mode
......@@ -112,15 +116,46 @@ opts_simul.max_check_ahead_periods = options_.occbin.smoother.max_check_ahead_pe
opts_simul.periodic_solution = options_.occbin.smoother.periodic_solution;
opts_simul.full_output = options_.occbin.smoother.full_output;
opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only;
% init_mode = options_.occbin.smoother.init_mode; % 0 = standard; 1 = unconditional frcsts zero shocks+smoothed states in each period
% init_mode = 0;
occbin_options = struct();
occbin_options.first_period_occbin_update = options_.occbin.smoother.first_period_occbin_update;
occbin_options.opts_regime = opts_simul; % this builds the opts_simul options field needed by occbin.solver
occbin_options.opts_regime.binding_indicator = options_.occbin.likelihood.init_binding_indicator;
occbin_options.opts_regime.regime_history=options_.occbin.likelihood.init_regime_history;
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
occbin_options.opts_simul = opts_simul; % this builds the opts_simul options field needed by occbin.solver
occbin_options.opts_regime.binding_indicator = options_.occbin.smoother.init_binding_indicator;
occbin_options.opts_regime.regime_history=options_.occbin.smoother.init_regime_history;
options_.noprint = true;
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0,~,error_indicator] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
if error_indicator(1) || isempty(alphahat0)
if ~options_.occbin.smoother.linear_smoother || nargin~=12 %make sure linear smoother results are set before using them
options_.occbin.smoother.status=false;
[~,etahat,~,~,~,~,~,~,~,~,~,~,~,~,~,~,alphahat0] = ...
DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
options_.occbin.smoother.status=true;
else
etahat= oo_.occbin.linear_smoother.etahat;
alphahat0= oo_.occbin.linear_smoother.alphahat0;
end
base_regime = struct();
if M_.occbin.constraint_nbr==1
base_regime.regime = 0;
base_regime.regimestart = 1;
else
base_regime.regime1 = 0;
base_regime.regimestart1 = 1;
base_regime.regime2 = 0;
base_regime.regimestart2 = 1;
end
oo_.occbin.smoother.regime_history = [];
for jper=1:size(alphahat,2)+1
if jper == 1
oo_.occbin.smoother.regime_history = base_regime;
else
oo_.occbin.smoother.regime_history(jper) = base_regime;
end
end
end
oo_.occbin.smoother.realtime_regime_history = oo_.occbin.smoother.regime_history;
regime_history = oo_.occbin.smoother.regime_history;
......@@ -136,17 +171,16 @@ occbin_options.first_period_occbin_update = inf;
opts_regime.binding_indicator=[];
regime_history0 = regime_history;
fprintf('Occbin smoother iteration 1.\n')
disp_verbose('OccBin smoother iteration 1.',options_.verbosity)
opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
opts_simul.exo_pos = 1:M_.exo_nbr;
opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
opts_simul.periods = size(opts_simul.SHOCKS,1);
options_.occbin.simul=opts_simul;
options_.noprint = true;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
print_info(out.error_flag, options_.noprint, options_)
disp_verbose('OccBin smoother:: simulation within smoother did not converge.',options_.verbosity)
oo_.occbin.smoother.error_flag=321;
return;
end
......@@ -156,15 +190,10 @@ if options_.smoother_redux
[T0,R0] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
oo_.occbin.linear_smoother.T0=T0;
oo_.occbin.linear_smoother.R0=R0;
% oo_.occbin.linear_smoother.T0=ss.T(oo_.dr.order_var,oo_.dr.order_var,1);
% oo_.occbin.linear_smoother.R0=ss.R(oo_.dr.order_var,:,1);
end
TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:);
RR = ss.R(oo_.dr.order_var,:,:);
CC = ss.C(oo_.dr.order_var,:);
% TT = cat(3,TT(:,:,1),TT);
% RR = cat(3,RR(:,:,1),RR);
% CC = cat(2,CC(:,1),CC);
opts_regime.regime_history = regime_history;
opts_regime.binding_indicator = [];
......@@ -186,6 +215,7 @@ sto_etahat={etahat};
sto_CC = CC;
sto_RR = RR;
sto_TT = TT;
sto_eee=NaN(size(TT,1),size(TT,3));
for k=1:size(TT,3)
sto_eee(:,k) = eig(TT(:,:,k));
end
......@@ -193,13 +223,14 @@ end
while is_changed && maxiter>iter && ~is_periodic
iter=iter+1;
fprintf('Occbin smoother iteration %u.\n', iter)
disp_verbose(sprintf('OccBin smoother iteration %u.', iter),options_.verbosity)
occbin_options.opts_regime.regime_history=regime_history;
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT;
[alphahat,etahat,epsilonhat,~,SteadyState,trend_coeff,~,T0,R0,P,~,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0]...
= DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);
sto_etahat(iter)={etahat};
regime_history0(iter,:) = regime_history;
if occbin_smoother_debug
save('info1','regime_history0');
save('Occbin_smoother_debug_regime_history','regime_history0');
end
sto_CC = CC;
......@@ -211,8 +242,7 @@ while is_changed && maxiter>iter && ~is_periodic
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
print_info(out.error_flag, false, options_)
disp_verbose('OccBin smoother:: simulation within smoother did not converge.',options_.verbosity)
oo_.occbin.smoother.error_flag=321;
return;
end
......@@ -220,13 +250,13 @@ while is_changed && maxiter>iter && ~is_periodic
TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:);
RR = ss.R(oo_.dr.order_var,:,:);
CC = ss.C(oo_.dr.order_var,:);
% TT = cat(3,TT(:,:,1),TT);
% RR = cat(3,RR(:,:,1),RR);
% CC = cat(2,CC(:,1),CC);
opts_regime.regime_history = regime_history;
[TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, options_ , oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
is_changed = ~isequal(regime_history0(iter,:),regime_history);
isdiff_regime=NaN(size(regime_history0,2),M_.occbin.constraint_nbr);
isdiff_start=NaN(size(isdiff_regime));
isdiff_=NaN(size(isdiff_regime));
if M_.occbin.constraint_nbr==2
for k=1:size(regime_history0,2)
isdiff_regime(k,1) = ~isequal(regime_history0(end,k).regime1,regime_history(k).regime1);
......@@ -236,16 +266,16 @@ while is_changed && maxiter>iter && ~is_periodic
isdiff_start(k,2) = ~isequal(regime_history0(end,k).regimestart2,regime_history(k).regimestart2);
isdiff_(k,2) = isdiff_regime(k,2) || isdiff_start(k,2);
end
is_changed_regime = ~isempty(find(isdiff_regime(:,1))) || ~isempty(find(isdiff_regime(:,2)));
is_changed_start = ~isempty(find(isdiff_start(:,1))) || ~isempty(find(isdiff_start(:,2)));
is_changed_regime = any(isdiff_regime(:,1)) || any(isdiff_regime(:,2));
is_changed_start = any(isdiff_start(:,1)) || any(isdiff_start(:,2));
else
for k=1:size(regime_history0,2)
isdiff_regime(k,1) = ~isequal(regime_history0(end,k).regime,regime_history(k).regime);
isdiff_start(k,1) = ~isequal(regime_history0(end,k).regimestart,regime_history(k).regimestart);
isdiff_(k,1) = isdiff_regime(k,1) || isdiff_start(k,1);
end
is_changed_regime = ~isempty(find(isdiff_regime(:,1)));
is_changed_start = ~isempty(find(isdiff_start(:,1)));
is_changed_regime = any(isdiff_regime(:,1));
is_changed_start = any(isdiff_start(:,1));
end
if occbin_smoother_fast
is_changed = is_changed_regime;
......@@ -261,9 +291,11 @@ while is_changed && maxiter>iter && ~is_periodic
end
if is_changed
eee=NaN(size(TT,1),size(TT,3));
for k=1:size(TT,3)
eee(:,k) = eig(TT(:,:,k));
end
if options_.debug
err_eig(iter-1) = max(max(abs(sort(eee)-sort(sto_eee))));
err_alphahat(iter-1) = max(max(max(abs(alphahat-sto_alphahat))));
err_etahat(iter-1) = max(max(max(abs(etahat-sto_etahat{iter-1}))));
......@@ -271,6 +303,7 @@ while is_changed && maxiter>iter && ~is_periodic
err_RR(iter-1) = max(max(max(abs(RR-sto_RR))));
err_TT(iter-1) = max(max(max(abs(TT-sto_TT))));
end
end
if occbin_smoother_debug || is_periodic
regime_ = cell(0);
......@@ -350,26 +383,26 @@ regime_history0(max(iter+1,1),:) = regime_history;
oo_.occbin.smoother.regime_history=regime_history0(end,:);
oo_.occbin.smoother.regime_history_iter=regime_history0;
if occbin_smoother_debug
save('info1','regime_history0')
save('Occbin_smoother_debug_regime_history','regime_history0')
end
if (maxiter==iter && is_changed) || is_periodic
disp('occbin.DSGE_smoother: smoother did not converge.')
fprintf('occbin.DSGE_smoother: The algorithm did not reach a fixed point for the smoothed regimes.\n')
disp_verbose('occbin.DSGE_smoother: smoother did not converge.',options_.verbosity)
disp_verbose('occbin.DSGE_smoother: The algorithm did not reach a fixed point for the smoothed regimes.',options_.verbosity)
if is_periodic
oo_.occbin.smoother.error_flag=0;
fprintf('occbin.DSGE_smoother: For the periods indicated above, regimes loops between the "regime_" and the "regime_new_" pattern displayed above.\n')
fprintf('occbin.DSGE_smoother: We provide smoothed shocks consistent with "regime_" in oo_.\n')
disp_verbose('occbin.DSGE_smoother: For the periods indicated above, regimes loops between the "regime_" and the "regime_new_" pattern displayed above.',options_.verbosity)
disp_verbose('occbin.DSGE_smoother: We provide smoothed shocks consistent with "regime_" in oo_.',options_.verbosity)
else
fprintf('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.\n')
disp_verbose('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.',options_.verbosity)
oo_.occbin.smoother=[];
oo_.occbin.smoother.error_flag=322;
end
else
disp('occbin.DSGE_smoother: smoother converged.')
disp_verbose('occbin.DSGE_smoother: smoother converged.',options_.verbosity)
oo_.occbin.smoother.error_flag=0;
if occbin_smoother_fast && is_changed_start
disp('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge')
disp_verbose('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge',options_.verbosity)
end
end
if (~is_changed || occbin_smoother_debug) && nargin==12
......@@ -393,15 +426,27 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
oo_.occbin.smoother.T0=TT;
oo_.occbin.smoother.R0=RR;
oo_.occbin.smoother.C0=CC;
oo_.occbin.smoother.simul.piecewise = out.piecewise(1:end-1,:);
if ~options_.occbin.simul.piecewise_only
oo_.occbin.smoother.simul.linear = out.linear(1:end-1,:);
end
if options_.occbin.smoother.plot
GraphDirectoryName = CheckPath('graphs',M_.fname);
latexFolder = CheckPath('latex',M_.dname);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([latexFolder filesep M_.fname '_OccBin_smoother_plots.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by occbin.DSGE_smoother.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
j1=0;
ifig=0;
for j=1:M_.exo_nbr
if M_.Sigma_e(j,j)
if max(abs(oo_.occbin.smoother.etahat(j,:)))>1.e-8
j1=j1+1;
if mod(j1,9)==1
hh_fig = dyn_figure(options_.nodisplay,'name','Occbin smoothed shocks');
hh_fig = dyn_figure(options_.nodisplay,'name','OccBin smoothed shocks');
ifig=ifig+1;
isub=0;
end
......@@ -414,13 +459,29 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
plot(oo_.occbin.smoother.etahat(j,:)','r--','linewidth',2)
hold on, plot([0 options_.nobs],[0 0],'k--')
set(gca,'xlim',[0 options_.nobs])
title(deblank(M_.exo_names(j,:)),'interpreter','none')
if options_.TeX
title(['$' M_.exo_names_tex{j,:} '$'],'interpreter','latex')
else
title(M_.exo_names{j,:},'interpreter','none')
end
if mod(j1,9)==0
if options_.occbin.smoother.linear_smoother
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
end
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
dyn_saveas(gcf,[GraphDirectoryName filesep M_.fname,'_smoothedshocks_occbin',int2str(ifig)],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
% TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_smoothedshocks_occbin%s}\n',options_.figures.textwidth*min(j1/3,1),[GraphDirectoryName '/' M_.fname],int2str(ifig)); % don't use filesep as it will create issues with LaTeX on Windows
fprintf(fidTeX,'\\caption{OccBin smoothed shocks.}');
fprintf(fidTeX,'\\label{Fig:smoothedshocks_occbin:%s}\n',int2str(ifig));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
end
end
end
......@@ -428,9 +489,19 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
dyn_saveas(hh_fig,[GraphDirectoryName filesep M_.fname,'_smoothedshocks_occbin',int2str(ifig)],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
% TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_smoothedshocks_occbin%s}\n',options_.figures.textwidth*min(j1/3,1),[GraphDirectoryName '/' M_.fname],int2str(ifig)); % don't use filesep as it will create issues with LaTeX on Windows
fprintf(fidTeX,'\\caption{OccBin smoothed shocks.}');
fprintf(fidTeX,'\\label{Fig:smoothedshocks_occbin:%s}\n',int2str(ifig));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fclose(fidTeX);
end
end
end
matlab/+occbin/IVF_core.m
View file @ 4a32acc2
function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% Computes thre
function [filtered_errs, resids, Emat, stateval, error_code, regime_history] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% Calls the solver to get the shocks explaining the data for the inversion filter
%
% Outputs:
% - filtered_errs [T by N_obs] filtered shocks
......@@ -10,9 +10,12 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o
% - error_code [4 by 1] error code
%
% Inputs
% - M_ [structure] Matlab's structure describing the model (M_).
% - oo_ [structure] Matlab's structure containing the results (oo_).
% - options_ [structure] Matlab's structure describing the current options (options_).
% - M_ [structure] MATLAB's structure describing the model (M_).
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - options_ [structure] MATLAB's structure describing the current options (options_).
% - err_index [double] index of shocks with strictly positive variance in M_.exo_names
% - filtered_errs_init [T by N_obs] initial values for the shocks
% - my_obs_list [cell] names of observables
......@@ -25,7 +28,7 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o
% Adapted for Dynare by Dynare Team.
%
% This code is in the public domain and may be used freely.
% However the authors would appreciate acknowledgement of the source by
% However the authors would appreciate acknowledgment of the source by
% citation of any of the following papers:
%
% Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello, and Molin Zhong (2019): "Likelihood evaluation of models
......@@ -39,7 +42,7 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o
[sample_length, n_obs]= size(obs);
nerrs = size(err_index,1);
if nargin<8
if nargin<11
init_val = zeros(M_.endo_nbr,1);
end
......@@ -83,12 +86,15 @@ for this_period=1:sample_length
opts_simul.exo_pos=err_index(inan); %err_index is predefined mapping from observables to shocks
opts_simul.endo_init = init_val_old;
opts_simul.SHOCKS = filtered_errs_init(this_period,inan);
[err_vals_out, exitflag] = dynare_solve(@occbin.match_function, filtered_errs_init(this_period,inan)', options_.occbin.solver.maxit, options_.occbin.solver.solve_tolf, options_.occbin.solver.solve_tolx, options_, obs_list, current_obs, opts_simul, M_, oo_, options_);
[err_vals_out, exitflag] = dynare_solve(@occbin.match_function, filtered_errs_init(this_period,inan)', options_.occbin.solver.maxit, options_.occbin.solver.solve_tolf, options_.occbin.solver.solve_tolx, options_, obs_list, current_obs, opts_simul, M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_);
if exitflag
filtered_errs=NaN;
error_code(1) = 304;
error_code(4) = 1000;
if this_period == 1
regime_history(this_period) = [];
end
if options_.occbin.likelihood.waitbar; dyn_waitbar_close(hh_fig); end
return
end
......@@ -96,9 +102,10 @@ for this_period=1:sample_length
opts_simul.SHOCKS = err_vals_out;
[ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_] = occbin.match_function(...
err_vals_out,obs_list,current_obs,opts_simul, M_,oo_,options_);
[ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_, out] = occbin.match_function(...
err_vals_out,obs_list,current_obs,opts_simul, M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_);
init_val = stateval(this_period,:); %update
regime_history(this_period) = out.regime_history(1);
if max(abs(err_vals_out))>1e8
error_code(1) = 306;
error_code(4) = max(abs(err_vals_out))/1000;
......
matlab/+occbin/IVF_posterior.m
View file @ 4a32acc2
function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults, atT, innov] = IVF_posterior(xparam1,...
dataset_,obs_info,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults)
% function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults, atT, innov] = IVF_posterior(xparam1,...
% dataset_,obs_info,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults)
function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov, regime_history] = IVF_posterior(xparam1,...
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
% dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% Computes Likelihood with inversion filter
%
% INPUTS
% - xparam1 [double] current values for the estimated parameters.
% - dataset_ [structure] dataset after transformations
% - DynareOptions [structure] Matlab's structure describing the current options (options_).
% - Model [structure] Matlab's structure describing the model (M_).
% - EstimatedParameters [structure] characterizing parameters to be estimated
% - BayesInfo [structure] describing the priors
% - dataset_info [structure] storing information about the
% sample; not used but required for interface
% - options_ [structure] MATLAB's structure describing the current options
% - M_ [structure] MATLAB's structure describing the model
% - estim_params_ [structure] characterizing parameters to be estimated
% - bayestopt_ [structure] describing the priors
% - BoundsInfo [structure] containing prior bounds
% - DynareResults [structure] Matlab's structure containing the results (oo_).
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - fval [double] scalar, value of the likelihood or posterior kernel.
% - info [integer] 4×1 vector, informations resolution of the model and evaluation of the likelihood.
% - info [integer] 4×1 vector, information resolution of the model and evaluation of the likelihood.
% - exit_flag [integer] scalar, equal to 1 (no issues when evaluating the likelihood) or 0 (not able to evaluate the likelihood).
% - DLIK [double] Empty array.
% - Hess [double] Empty array.
% - SteadyState [double] Empty array.
% - trend [double] Empty array.
% - Model [struct] Updated Model structure described in INPUTS section.
% - DynareOptions [struct] Updated DynareOptions structure described in INPUTS section.
% - BayesInfo [struct] See INPUTS section.
% - DynareResults [struct] Updated DynareResults structure described in INPUTS section.
% - M_ [struct] Updated M_ structure described in INPUTS section.
% - options_ [struct] Updated options_ structure described in INPUTS section.
% - bayestopt_ [struct] See INPUTS section.
% - dr [structure] Reduced form model.
% - atT [double] (m*T) matrix, smoothed endogenous variables (a_{t|T}) (decision-rule order)
% - innov [double] (r*T) matrix, smoothed structural shocks (r>n is the umber of shocks).
% Copyright © 2021 Dynare Team
% Copyright © 2021-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -51,8 +56,10 @@ DLIK=[];
Hess=[];
trend_coeff = [];
obs = dataset_.data;
obs_list = DynareOptions.varobs(:);
obs_list = options_.varobs(:);
exit_flag = 1;
% initialize output argument in case of early return
regime_history = [];
if size(xparam1,1)<size(xparam1,2)
......@@ -65,18 +72,18 @@ end
%------------------------------------------------------------------------------
if ~isempty(xparam1)
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
[fval,info,exit_flag,Q,H]=check_bounds_and_definiteness_estimation(xparam1, Model, EstimatedParameters, BoundsInfo);
M_ = set_all_parameters(xparam1,estim_params_,M_);
[fval,info,exit_flag]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, BoundsInfo);
if info(1)
return
end
end
err_index=DynareOptions.occbin.likelihood.IVF_shock_observable_mapping; % err_index= find(diag(Model.Sigma_e)~=0);
COVMAT1 = Model.Sigma_e(err_index,err_index);
err_index=options_.occbin.likelihood.IVF_shock_observable_mapping; % err_index= find(diag(M_.Sigma_e)~=0);
COVMAT1 = M_.Sigma_e(err_index,err_index);
% Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
[T,R,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');
[~,~,SteadyState,info,dr, M_.params] = dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
......@@ -85,7 +92,11 @@ if info(1)
info(1) == 81 || info(1) == 84 || info(1) == 85 || info(1) == 86
%meaningful second entry of output that can be used
fval = Inf;
if ~isfinite(info(2))
info(4) = 0.1;
else
info(4) = info(2);
end
exit_flag = 0;
return
else
......@@ -102,17 +113,17 @@ end
sample_length = size(obs,1);
filtered_errs_init = zeros(sample_length,sum(err_index));
[filtered_errs, resids, Emat, stateval, info] = occbin.IVF_core(Model,DynareResults,DynareOptions,err_index,filtered_errs_init,obs_list,obs);
[filtered_errs, resids, Emat, stateval, info, regime_history] = occbin.IVF_core(M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,err_index,filtered_errs_init,obs_list,obs);
if info(1)
fval = Inf;
exit_flag = 0;
atT=NaN(size(stateval(:,DynareResults.dr.order_var)'));
innov=NaN(Model.exo_nbr,sample_length);
atT=NaN(size(stateval(:,dr.order_var)'));
innov=NaN(M_.exo_nbr,sample_length);
return
else
atT = stateval(:,DynareResults.dr.order_var)';
innov = zeros(Model.exo_nbr,sample_length);
innov(diag(Model.Sigma_e)~=0,:)=filtered_errs';
atT = stateval(:,dr.order_var)';
innov = zeros(M_.exo_nbr,sample_length);
innov(diag(M_.Sigma_e)~=0,:)=filtered_errs';
end
nobs=size(filtered_errs,1);
......@@ -121,7 +132,7 @@ nobs=size(filtered_errs,1);
% Calculate the selection matrix
%-------------------------------------
iobs=DynareOptions.varobs_id';
iobs=options_.varobs_id';
likei=NaN(nobs,1);
if ~any(any(isnan(obs)))
......@@ -145,42 +156,42 @@ else
end
end
like = 0.5*sum(likei(DynareOptions.presample+1:end));
like = 0.5*sum(likei(options_.presample+1:end));
if isinf(like)
fval = Inf;
info(1) = 301;
info(4) = 1000;
exit_flag = 0;
fval = Inf; info(1) = 301; info(4) = 1000; exit_flag = 0;
return
elseif isnan(like)
fval = Inf;
info(1) = 302;
info(4) = 1000;
exit_flag = 0;
fval = Inf; info(1) = 302; info(4) = 1000; exit_flag = 0;
return
elseif imag(like)~=0
fval = Inf; info(1) = 300; info(4) = 1000; exit_flag = 0;
return
end
maxresid = max(abs(resids(:)));
if maxresid>1e-3
disp_verbose('Penalize failure of residuals to be zero',options_.verbosity)
fval = Inf;
info(1) = 303;
info(4) = sum(resids(:).^2);
exit_flag = 0;
fval = Inf; info(1) = 303; info(4) = sum(resids(:).^2); exit_flag = 0;
return
end
if ~isempty(xparam1)
prior = -priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
prior = -priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
else
prior = 0;
end
if prior == Inf
if isinf(prior)
% If parameters outside prior bound, minus prior density is very large
fval = Inf;
info(4) = 1000;
exit_flag = 0;
fval = Inf; info(1) = 40; info(4) = 1000; exit_flag = 0;
return
end
if isnan(prior)
fval = Inf; info(1) = 47; info(4) = 1000; exit_flag = 0;
return
end
if imag(prior)~=0
fval = Inf; info(1) = 48; info(4) = 1000; exit_flag = 0;
return
end
......
matlab/+occbin/backward_map_regime.m 0 → 100644
View file @ 4a32acc2
function [binding_indicator, A, regime_string] = backward_map_regime(regime, regime_start)
% [binding_indicator, A, regime_string] = backward_map_regime(regime, regime_start)
% Map regime information into regime indicator
%
% Inputs:
% - regime [integer] [1 by n_transitions] vector of regime number indices
% - regime_start [integer] [1 by n_transitions] vectors with period numbers in which regime starts
%
% Outputs:
% - binding_indicator [integer] [nperiods by 1] vector of regime indices
% - A [bin] binary representation of binding indicator
% - error_flag [boolean] 1 if regime never leaves 1 or is still there at the end of nperiods
% 0 otherwise
% Copyright © 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
if nargin ==1
% polymorphism
if ~isstruct(regime)
disp('error::backward_map_regime')
disp('input arguments may be 1 structure with regime info')
disp('or two arrays: regime and regimestart ')
error('wrong input')
end
fnam = fieldnames(regime);
if length(fnam) == 2
[binding_indicator, A, regime_string] = occbin.backward_map_regime(regime.regime, regime.regimestart);
else
for k=1:2
nperiods(k) = regime.(['regimestart' int2str(k)])(end);
number_of_binary_tokens(k) = ceil((nperiods(k)-1)/50);
end
binding_indicator = false(max(nperiods),2);
A = int64(zeros(max(number_of_binary_tokens),2));
for k=1:2
[binding_indicator(1:nperiods(k),k), A(1:number_of_binary_tokens(k),k), tmp{k}] = ...
occbin.backward_map_regime(regime.(['regime' int2str(k)]), regime.(['regimestart' int2str(k)]));
end
regime_string = char(tmp{1},tmp{2});
end
return
else
if isstruct(regime)
disp('error::backward_map_regime')
disp('input arguments may be ONE structure with regime info')
disp('or TWO arrays: regime and regimestart ')
error('wrong input')
end
end
regime_string = char(mat2str(double(regime)),mat2str(regime_start));
nperiods_0 = regime_start(end);
number_of_binary_tokens = max(1,ceil((nperiods_0-1)/50));
A = int64(zeros(number_of_binary_tokens,1));
binding_indicator = false(nperiods_0,1);
if length(regime)>1
for ir=1:length(regime)-1
binding_indicator(regime_start(ir):regime_start(ir+1)-1,1) = regime(ir);
for k=regime_start(ir):regime_start(ir+1)-1
this_token = ceil(k/50);
A(this_token) = int64(bitset(A(this_token),k-50*(this_token-1),regime(ir)));
end
end
end
binding_indicator = logical(binding_indicator);
% to convert regime in a readable string array
% a = dec2bin(A);
% bindicator = [a(end:-1:1) '0'];
\ No newline at end of file
matlab/+occbin/check_regimes.m
View file @ 4a32acc2
......@@ -6,8 +6,8 @@ function [TT, RR, CC, regime_history] = check_regimes(TT, RR, CC, opts_regime, M
% - RR [N by N_exo] shock impact matrix of state space
% - CC [N by 1] constant of state space
% - opts_regime_ [structure] structure describing the regime
% - M_ [structure] Matlab's structure describing the model
% - options_ [structure] Matlab's structure describing the current options
% - M_ [structure] MATLAB's structure describing the model
% - options_ [structure] MATLAB's structure describing the current options
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
......
matlab/+occbin/cost_function.m 0 → 100644
View file @ 4a32acc2
function [cost, out] = cost_function(err_0, current_obs, weights, opts_simul,...
M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_)
% [cost, out] = cost_function(err_0, current_obs, opts_simul,...
% M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_)
% Outputs:
% - cost [double] penalty
% - out [structure] OccBin's results structure
%
% Inputs
% - err_0 [double] value of shocks
% - current_obs [double] [1 by n_obs] current value of observables
% - weights [double] [1 by n_obs] variance of observables,
% - opts_simul [structure] Structure with simulation options
% used in cost function
% - M_ [structure] MATLAB's structure describing the model (M_).
% - dr_ [structure] model information structure
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - options_ [structure] MATLAB's structure describing the current options (options_).
% Copyright © 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
opts_simul.SHOCKS = err_0';
options_.occbin.simul=opts_simul;
options_.occbin.simul.full_output=1;
options_.noprint = 1;
[~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
cost = 0;
if ~out.error_flag
cost = mean((out.piecewise(1,opts_simul.varobs_id)'-current_obs').^2./weights);
else
cost = cost+1.e10;
end
\ No newline at end of file
matlab/+occbin/dynare_resolve.m
View file @ 4a32acc2
function [A,B,ys,info,dr,params,TT, RR, CC, A0, B0] ...
= dynare_resolve(M_,options_,oo_,regime_history, reduced_state_space, A, B)
% function [A,B,ys,info,M_,options_,oo_,TT, RR, CC, A0, B0] ...
% = dynare_resolve(M_,options_,oo_,regime_history, reduced_state_space, A, B)
= dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,regime_history, reduced_state_space, A, B)
% [A,B,ys,info,M_,options_,oo_,TT, RR, CC, A0, B0] ...
% = dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,regime_history, reduced_state_space, A, B)
% Computes the linear approximation and the matrices A and B of the
% transition equation. Mirrors dynare_resolve
%
% Inputs:
% - M_ [structure] Matlab's structure describing the model
% - options_ [structure] Matlab's structure containing the options
% - oo_ [structure] Matlab's structure containing the results
% - M_ [structure] MATLAB's structure describing the model
% - options_ [structure] MATLAB's structure containing the options
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - reduced_state_space [string]
% - A [double] State transition matrix
% - B [double] shock impact matrix
......@@ -26,7 +29,7 @@ function [A,B,ys,info,dr,params,TT, RR, CC, A0, B0] ...
% - A0 [double] State transition matrix (unrestricted state space)
% - B0 [double] shock impact matrix (unrestricted state space)
% Copyright © 2001-2021 Dynare Team
% Copyright © 2001-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -43,19 +46,18 @@ function [A,B,ys,info,dr,params,TT, RR, CC, A0, B0] ...
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if nargin<6
[A,B,ys,info,dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if nargin<9
[A,B,ys,info,dr,M_.params] = dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
else
ys = oo_.dr.ys;
dr = oo_.dr;
ys = dr.ys;
info = 0;
end
params=M_.params;
if ~info(1) && nargin>4 && ~isempty(regime_history)
if ~info(1) && nargin>7 && ~isempty(regime_history)
opts_regime.regime_history=regime_history;
opts_regime.binding_indicator=[];
[TT, RR, CC] = ...
occbin.check_regimes(A, B, [], opts_regime, M_, options_, dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
occbin.check_regimes(A, B, [], opts_regime, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
else
TT=A;
RR=B;
......@@ -64,7 +66,7 @@ end
A0=A;
B0=B;
if ~info(1) && nargin>4 && ischar(reduced_state_space) && ~isempty(reduced_state_space)
if ~info(1) && nargin>7 && ischar(reduced_state_space) && ~isempty(reduced_state_space)
iv = dr.restrict_var_list;
A=A(iv,iv);
B=B(iv,:);
......
matlab/+occbin/findmin.m 0 → 100644
View file @ 4a32acc2
function [y, out, cost] = findmin(d_index, a0, P1, Qt, Y, ZZ, opts_simul,M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_)
% [y, out, cost] = findmin(d_index, a0, P1, Qt, Y, ZZ, opts_simul,M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_)
% Outputs:
% - cost [double] penalty
% - out [structure] OccBin's results structure
%
% Inputs
% - opts_simul [structure] Structure with simulation options
% used in cost function
% - M_ [structure] MATLAB's structure describing the model
% - dr_ [structure] model information structure
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - options_ [structure] MATLAB's structure describing the current options
% Copyright © 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
current_obs = Y(d_index,2)'+dr.ys(options_.varobs_id(d_index))';
err_index = find(diag(Qt(:,:,2))~=0);
F = ZZ(d_index,:)*P1(:,:,2)*ZZ(d_index,:)' ;
weights=diag(F);
filtered_errs_init = zeros(1,length(err_index));
opts_simul.varobs_id=options_.varobs_id(d_index)';
opts_simul.exo_pos=err_index; %err_index is predefined mapping from observables to shocks
opts_simul.SHOCKS = filtered_errs_init;
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(dr.restrict_var_list,1)=a0(:,1); %updated state
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = a0(dr.inv_order_var,1);
end
[y] = fminsearch(@cost_function,filtered_errs_init');
[cost, out] = occbin.cost_function(y, current_obs, weights, opts_simul,...
M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_);
function cost = cost_function(x)
cost = occbin.cost_function(x, current_obs, weights, opts_simul,...
M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_);
end
end
\ No newline at end of file
matlab/+occbin/forecast.m 0 → 100644
View file @ 4a32acc2
function [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,forecast_horizon)
% [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,forecast_horizon)
% Occbin forecasts
%
% INPUTS
% - options_ [structure] MATLAB's structure describing the current options
% - M_ [structure] MATLAB's structure describing the model
% - dr_in [structure] model information structure
% - endo_steady_state [double] steady state value for endogenous variables
% - exo_steady_state [double] steady state value for exogenous variables
% - exo_det_steady_state [double] steady state value for exogenous deterministic variables
% - forecast_horizon [integer] forecast horizon
%
% OUTPUTS
% - forecast [structure] forecast results
% - error_flag [integer] error code
%
% SPECIAL REQUIREMENTS
% none.
% Copyright © 2022-2024 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/>.
opts = options_.occbin.forecast;
options_.occbin.simul.maxit = opts.maxit;
options_.occbin.simul.check_ahead_periods = opts.check_ahead_periods;
options_.occbin.simul.periods = forecast_horizon;
shocks_input = opts.SHOCKS0;
if ~isempty(shocks_input)
n_shocks=size(shocks_input,1);
if iscell(shocks_input)
inds=NaN(n_shocks,1);
periods=length(shocks_input{1}{2});
shock_mat=NaN(n_shocks,periods);
for j=1:n_shocks
exo_pos=strmatch(shocks_input{j}{1},M_.exo_names,'exact');
if isempty(exo_pos)
error('occbin.forecast: unknown exogenous shock %s',shocks_input{j}{1})
else
inds(j)=exo_pos;
end
if length(shocks_input{j}{2})~=periods
error('occbin.forecast: unknown exogenous shock %s',shocks_input{j}{1})
else
shock_mat(j,:)=shocks_input{j}{2};
end
end
elseif isreal(shocks_input)
shock_mat=shocks_input;
inds = (1:M_.exo_nbr)';
end
end
options_.occbin.simul.endo_init = M_.endo_histval(:,1)-endo_steady_state; %initial condition
options_.occbin.simul.init_regime = opts.frcst_regimes;
options_.occbin.simul.init_binding_indicator = [];
shocks_base = zeros(forecast_horizon,M_.exo_nbr);
if ~isempty(shocks_input)
for j=1:n_shocks
shocks_base(:,inds(j))=shock_mat(j,:);
end
end
if opts.replic
options_.noprint=true;
h = dyn_waitbar(0,'Please wait occbin forecast replic ...');
ishock = find(sqrt(diag((M_.Sigma_e))));
options_.occbin.simul.exo_pos=ishock;
effective_exo_nbr= length(ishock);
effective_Sigma_e = M_.Sigma_e(ishock,ishock); % does not take heteroskedastic shocks into account
[U,S] = svd(effective_Sigma_e);
% draw random shocks
if opts.qmc
opts.replic =2^(round(log2(opts.replic+1)))-1;
SHOCKS_add = qmc_sequence(forecast_horizon*effective_exo_nbr, int64(1), 1, opts.replic);
else
SHOCKS_add = randn(forecast_horizon*effective_exo_nbr,opts.replic);
end
SHOCKS_add=reshape(SHOCKS_add,effective_exo_nbr,forecast_horizon,opts.replic);
z.linear=NaN(forecast_horizon,M_.endo_nbr,opts.replic);
z.piecewise=NaN(forecast_horizon,M_.endo_nbr,opts.replic);
error_flag=true(opts.replic,1);
simul_SHOCKS=NaN(forecast_horizon,M_.exo_nbr,opts.replic);
for iter=1:opts.replic
options_.occbin.simul.SHOCKS = shocks_base(:,ishock)+transpose(U*sqrt(S)*SHOCKS_add(:,:,iter));
options_.occbin.simul.waitbar=0;
[~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
error_flag(iter)=out.error_flag;
if ~error_flag(iter)
z.linear(:,:,iter)=out.linear;
z.piecewise(:,:,iter)=out.piecewise;
frcst_regime_history(iter,:)=out.regime_history;
error_flag(iter)=out.error_flag;
simul_SHOCKS(:,:,iter) = shocks_base;
else
if options_.debug
save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag','out','shocks_base')
end
end
dyn_waitbar(iter/opts.replic,h,['OccBin MC forecast replic ',int2str(iter),'/',int2str(opts.replic)])
end
dyn_waitbar_close(h);
if options_.debug
save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag')
end
inx=find(error_flag==0);
if length(inx)<opts.replic
fprintf('\noccbin.forecast: forecast simulation was only successful in %u of %u simulations.\n',length(inx),opts.replic);
end
z.linear=z.linear(:,:,inx);
z.piecewise=z.piecewise(:,:,inx);
z.min.piecewise = min(z.piecewise,[],3);
z.max.piecewise = max(z.piecewise,[],3);
z.min.linear = min(z.linear,[],3);
z.max.linear = max(z.linear,[],3);
field_names={'linear','piecewise'};
post_mean=NaN(forecast_horizon,1);
post_median=NaN(forecast_horizon,1);
post_var=NaN(forecast_horizon,1);
hpd_interval=NaN(forecast_horizon,2);
post_deciles=NaN(forecast_horizon,9);
for field_iter=1:2
for i=1:M_.endo_nbr
for j=1:forecast_horizon
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(z.(field_names{field_iter})(j,i,:)),options_.forecasts.conf_sig);
end
forecast.(field_names{field_iter}).Mean.(M_.endo_names{i})=post_mean;
forecast.(field_names{field_iter}).Median.(M_.endo_names{i})=post_median;
forecast.(field_names{field_iter}).Var.(M_.endo_names{i})=post_var;
forecast.(field_names{field_iter}).HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
forecast.(field_names{field_iter}).HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
forecast.(field_names{field_iter}).Deciles.(M_.endo_names{i})=post_deciles;
forecast.(field_names{field_iter}).Min.(M_.endo_names{i})=z.min.(field_names{field_iter})(:,i);
forecast.(field_names{field_iter}).Max.(M_.endo_names{i})=z.max.(field_names{field_iter})(:,i);
end
end
else
options_.occbin.simul.irfshock = M_.exo_names;
options_.occbin.simul.SHOCKS = shocks_base;
[~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
error_flag=out.error_flag;
if out.error_flag
fprintf('occbin.forecast: forecast simulation failed.')
return;
end
frcst_regime_history=out.regime_history;
error_flag=out.error_flag;
for i=1:M_.endo_nbr
forecast.linear.Mean.(M_.endo_names{i})= out.linear(:,i);
forecast.piecewise.Mean.(M_.endo_names{i})= out.piecewise(:,i);
end
end
forecast.regimes=frcst_regime_history;
\ No newline at end of file
matlab/+occbin/get_deriv.m
View file @ 4a32acc2
......@@ -13,7 +13,7 @@ function [h_minus_1, h, h_plus_1, h_exo, resid] = get_deriv(M_, ys_)
% - h_exo [N by N_exo] derivative matrix with respect to exogenous variables
% - resid [N by 1] vector of residuals
% Copyright © 2021 Dynare Team
% Copyright © 2021-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -30,40 +30,16 @@ function [h_minus_1, h, h_plus_1, h_exo, resid] = get_deriv(M_, ys_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
x = zeros(M_.maximum_lag + M_.maximum_lead + 1,M_.exo_nbr);
dyn_endo_ss = repmat(ys_, 3, 1);
x = zeros(M_.exo_nbr);
iyv = M_.lead_lag_incidence';
iyr0 = find(iyv(:)) ;
z=repmat(ys_,1,M_.maximum_lag + M_.maximum_lead + 1);
[resid, T_order, T] = feval([M_.fname '.sparse.dynamic_resid'], dyn_endo_ss, x, M_.params, ys_);
[resid,g1]=feval([M_.fname,'.dynamic'],z(iyr0),x, M_.params, ys_, M_.maximum_exo_lag+1);
g1 = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, x, M_.params, ys_, ...
M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
M_.dynamic_g1_sparse_colptr, T_order, T);
% Initialize matrices
h_minus_1=zeros(M_.endo_nbr);
h = h_minus_1;
h_plus_1 = h_minus_1;
% build h_minus_1
if M_.maximum_lag
lag_columns=find(iyv(:,1));
n_lag_columns=length(lag_columns);
h_minus_1(:,lag_columns) = g1(:,1:n_lag_columns);
else
n_lag_columns=0;
end
% build h
contemporaneous_columns=find(iyv(:,M_.maximum_lag+1));
n_contemporaneous_columns = length(contemporaneous_columns);
h(:,contemporaneous_columns) = g1(:,1+n_lag_columns:n_lag_columns+n_contemporaneous_columns);
%build h_plus_1
if M_.maximum_lead
lead_columns=find(iyv(:,end));
n_lead_columns = length(lead_columns);
h_plus_1(:,lead_columns) = g1(:,n_lag_columns+n_contemporaneous_columns+1:n_lag_columns+n_contemporaneous_columns+n_lead_columns);
else
n_lead_columns=0;
end
h_exo =g1(:,n_lag_columns+n_contemporaneous_columns+n_lead_columns+1:end);
h_minus_1 = full(g1(:,1:M_.endo_nbr));
h = full(g1(:,M_.endo_nbr + (1:M_.endo_nbr)));
h_plus_1 = full(g1(:,2*M_.endo_nbr + (1:M_.endo_nbr)));
h_exo = full(g1(:,3*M_.endo_nbr+1:end));
matlab/+occbin/graph.m
View file @ 4a32acc2
......@@ -2,10 +2,10 @@ function graph(M_, options_, options_occbin_, oo_, var_list)
% function graph(M_, options_, options_occbin_, oo_, var_list)
%
% Inputs:
% - M_ [structure] Matlab's structure describing the model
% - options_ [structure] Matlab's structure containing the options
% - options_occbin_ [structure] Matlab's structure containing Occbin options
% - oo_ [structure] Matlab's structure containing the results
% - M_ [structure] MATLAB's structure describing the model
% - options_ [structure] MATLAB's structure containing the options
% - options_occbin_ [structure] MATLAB's structure containing OccBin options
% - oo_ [structure] MATLAB's structure containing the results
% - var_list [char] list of the variables to plot
% Copyright © 2021-2023 Dynare Team
......@@ -69,7 +69,7 @@ if number_of_plots_to_draw_exo>0
if ~isempty(temp_index)
exo_index(ii)=temp_index;
else
error('%s was not part of the shocks for Occbin.', var_list{i_var_exo(ii)});
error('%s was not part of the shocks for OccBin.', var_list{i_var_exo(ii)});
end
end
data_to_plot(:,end+1:end+number_of_plots_to_draw_exo,1)=[oo_.occbin.simul.shocks_sequence(:,exo_index); zeros(nperiods-size(oo_.occbin.simul.shocks_sequence,1),number_of_plots_to_draw_exo)];
......@@ -79,7 +79,7 @@ end
[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw_endo+number_of_plots_to_draw_exo);
for fig = 1:nbplt
hh_fig = dyn_figure(options_.nodisplay,'Name',['Occbin simulated paths, figure ' int2str(fig)]);
hh_fig = dyn_figure(options_.nodisplay,'Name',['OccBin simulated paths, figure ' int2str(fig)]);
for plt = 1:nstar
if fig==nbplt && ~lr==0
subplot(lr,lc,plt);
......@@ -109,7 +109,7 @@ for fig = 1:nbplt
if ndim==2
legend([h1,h2],legend_list,'box','off')
else
legend([h1],legend_list,'box','off')
legend(h1,legend_list,'box','off')
end
end
if options_.TeX
......
matlab/+occbin/irf.m 0 → 100644
View file @ 4a32acc2
function irfs = irf(M_,oo_,options_)
% irfs = irf(M_,oo_,options_)
% Calls a minimizer
%
% INPUTS
% - M_ [structure] MATLAB's structure describing the model
% - oo_ [structure] MATLAB's structure containing the results
% - options_ [structure] MATLAB's structure describing the current options
%
% OUTPUTS
% - irfs [structure] IRF results
%
% SPECIAL REQUIREMENTS
% none.
%
%
% Copyright © 2022-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/>.
shocknames = options_.occbin.irf.exo_names;
shocksigns = options_.occbin.irf.shocksigns; %'pos','neg'
shocksize = options_.occbin.irf.shocksize;
t_0 = options_.occbin.irf.t0;
%% set simulation options based on IRF options
options_.occbin.simul.init_regime = options_.occbin.irf.init_regime;
options_.occbin.simul.check_ahead_periods = options_.occbin.irf.check_ahead_periods;
options_.occbin.simul.maxit = options_.occbin.irf.maxit;
options_.occbin.simul.periods = options_.irf;
%% Run initial conditions + other shocks
if t_0 == 0
shocks_base = zeros(options_.occbin.simul.periods+1,M_.exo_nbr);
options_.occbin.simul.endo_init = [];
else
if ~isfield(oo_.occbin,'smoother')
error('occbin.irfs: smoother must be run before requesting GIRFs based on smoothed results')
end
% GIRF conditional on smoothed states in t_0 and shocks in t_0+1
shocks_base= [oo_.occbin.smoother.etahat(:,t_0+1)'; zeros(options_.occbin.simul.periods,M_.exo_nbr)];
options_.occbin.simul.SHOCKS=shocks_base;
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
end
options_.occbin.simul.SHOCKS=shocks_base;
[~, out_base] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out_base.error_flag
error('occbin.irfs: could not compute the solution')
end
irfs.linear = struct();
irfs.piecewise = struct();
% Get indices of shocks of interest
exo_index =zeros(size(shocknames,1),1);
for i=1:length(shocknames)
exo_index(i) = strmatch(shocknames{i},M_.exo_names,'exact');
end
% cs=get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
% irf_shocks_indx = getIrfShocksIndx(M_, options_);
% Set shock size
if isempty(shocksize)
shocksize = sqrt(diag(M_.Sigma_e(exo_index,exo_index)));
if any(shocksize < 1.e-9)
shocksize(shocksize < 1.e-9) = 0.01;
end
end
if numel(shocksize)==1
shocksize=repmat(shocksize,[length(shocknames),1]);
end
% Run IRFs
for sign_iter=1:length(shocksigns)
for IRF_counter = 1:length(exo_index)
jexo = exo_index(IRF_counter);
if ~options_.noprint && options_.debug
fprintf('occbin.irf: Producing GIRFs for shock %s. Simulation %d out of %d. \n',M_.exo_names{jexo},IRF_counter,size(exo_index,1));
end
shocks1=shocks_base;
if ismember('pos',shocksigns{sign_iter})
shocks1(1,jexo)=shocks_base(1,jexo)+shocksize(IRF_counter);
elseif ismember('neg',shocksigns{sign_iter})
shocks1(1,jexo)=shocks_base(1,jexo)-shocksize(IRF_counter);
end
options_.occbin.simul.SHOCKS=shocks1;
if t_0 == 0
options_.occbin.simul.endo_init = [];
else
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
end
[~, out_sim] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out_sim.error_flag
warning('occbin.irfs: simulation failed')
skip
end
% Substract inital conditions + other shocks
zdiff.linear.(shocksigns{sign_iter}) = out_sim.linear-out_base.linear;
zdiff.piecewise.(shocksigns{sign_iter}) = out_sim.piecewise-out_base.piecewise;
for j_endo=1:M_.endo_nbr
if ismember('pos',shocksigns)
irfs.piecewise.([M_.endo_names{j_endo} '_' M_.exo_names{jexo} '_' shocksigns{sign_iter}]) = zdiff.piecewise.(shocksigns{sign_iter})(:,j_endo);
irfs.linear.([M_.endo_names{j_endo} '_' M_.exo_names{jexo} '_' shocksigns{sign_iter}]) = zdiff.linear.(shocksigns{sign_iter})(:,j_endo);
end
end
end
end
\ No newline at end of file
matlab/+occbin/kalman_update_algo_1.m
View file @ 4a32acc2
function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,oo_,options_,occbin_options)
% function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,oo_,options_,occbin_options)
function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat, alphahat, V] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options)
% function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options)
% INPUTS
% - a [N by 1] t-1's state estimate
% - a1 [N by 2] state predictions made at t-1:t
......@@ -17,10 +17,13 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal
% - RR [N by N_exo by 2] shock impact matrix at t-1:t
% - CC [N by 2] state space constant state transition matrix at t-1:t
% - regimes0 [structure] regime info at t-1:t
% - M_ [structure] Matlab's structure describing the model (M_).
% - oo_ [structure] Matlab's structure containing the results (oo_).
% - options_ [structure] Matlab's structure describing the current options (options_).
% - occbin_options_ [structure] Matlab's structure describing the Occbin options.
% - M_ [structure] MATLAB's structure describing the model
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - options_ [structure] MATLAB's structure describing the current options
% - occbin_options_ [structure] MATLAB's structure describing the OccBin options.
% - kalman_tol [double] tolerance for reciprocal condition number
%
% Outputs
......@@ -34,7 +37,7 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal
% - C [N by 2] state space constant state transition matrix at t-1:t
% - regimes_ [structure] regime info at t-1:t
% - error_flag [integer] error code
% - M_ [structure] Matlab's structure describing the model (M_).
% - M_ [structure] MATLAB's structure describing the model
% - lik [double] likelihood
% - etahat: smoothed shocks
%
......@@ -42,7 +45,7 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal
% Philipp Pfeiffer, Marco Ratto (2021), Efficient and robust inference of models with occasionally binding
% constraints, Working Papers 2021-03, Joint Research Centre, European Commission
% Copyright © 2021-2022 Dynare Team
% Copyright © 2021-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -59,7 +62,7 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
warning off
orig_warning_state = warning;
options_.noprint = true;
R=NaN(size(RR));
......@@ -83,6 +86,8 @@ else
base_regime.regimestart2 = 1;
end
regimes_ = [base_regime base_regime base_regime];
opts_simul = occbin_options.opts_simul;
options_.occbin.simul=opts_simul;
mm=size(a,1);
%% store info in t=1
......@@ -100,45 +105,86 @@ PZI = P1(:,:,t)*ZZ'*iF(di,di,t);
% L(:,:,t) = T-K(:,di,t)*ZZ;
L(:,:,t) = eye(mm)-PZI*ZZ;
if ~options_.occbin.filter.use_relaxation
[a, a1, P, P1, v, alphahat, etahat, lik, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
if ~isempty(fieldnames(regimes0))
if options_.occbin.filter.guess_regime
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, regimes0(1),'reduced_state_space',T0,R0);
if M_.occbin.constraint_nbr==1
bindx = occbin.backward_map_regime(regimes0(1).regime, regimes0(1).regimestart);
bindx = bindx(2:end);
[regimes0(2).regime, regimes0(2).regimestart, error_flag]=occbin.map_regime(bindx);
bindx = bindx(2:end);
[regimes0(3).regime, regimes0(3).regimestart, error_flag]=occbin.map_regime(bindx);
else
bindx1 = occbin.backward_map_regime(regimes0(1).regime1, regimes0(1).regimestart1);
bindx2 = occbin.backward_map_regime(regimes0(1).regime2, regimes0(1).regimestart2);
bindx1 = bindx1(2:end);
bindx2 = bindx2(2:end);
[regimes0(2).regime1, regimes0(2).regimestart1, error_flag]=occbin.map_regime(bindx1);
[regimes0(2).regime2, regimes0(2).regimestart2, error_flag]=occbin.map_regime(bindx2);
bindx1 = bindx1(2:end);
bindx2 = bindx2(2:end);
[regimes0(3).regime1, regimes0(3).regimestart1, error_flag]=occbin.map_regime(bindx1);
[regimes0(3).regime2, regimes0(3).regimestart2, error_flag]=occbin.map_regime(bindx2);
end
% regimes0=[regimes0 base_regime base_regime];
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
end
[a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
else
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, base_regime,'reduced_state_space',T0,R0);
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, base_regime,'reduced_state_space',T0,R0);
if isempty(fieldnames(regimes0))
regimes0 = regimes_;
else
regimes0(1)=base_regime;
end
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
[a, a1, P, P1, v, alphahat, etahat, lik, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
[a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
end
if error_flag
etahat=NaN(size(QQQ,1),1);
warning(orig_warning_state);
return;
end
%% run here the occbin simul
opts_simul = occbin_options.opts_simul;
opts_simul.SHOCKS = zeros(3,M_.exo_nbr);
opts_simul.exo_pos=1:M_.exo_nbr;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
my_order_var = oo_.dr.order_var(oo_.dr.restrict_var_list);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
my_order_var = dr.order_var(dr.restrict_var_list);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
my_order_var = oo_.dr.order_var;
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
my_order_var = dr.order_var;
end
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if options_.occbin.filter.guess_regime
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
else
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
end
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
warning(orig_warning_state);
return;
end
......@@ -156,66 +202,34 @@ end
lik_hist=lik;
niter=1;
is_periodic=0;
if options_.occbin.filter.use_relaxation || isequal(regimes0(1),base_regime)
nguess=1;
else
nguess=0;
end
newguess=0;
if any(myregime) || ~isequal(regimes_(1),regimes0(1))
while ~isequal(regimes_(1),regimes0(1)) && ~is_periodic && ~out.error_flag && niter<=options_.occbin.likelihood.max_number_of_iterations
niter=niter+1;
oldstart=1;
if M_.occbin.constraint_nbr==1 && length(regimes0(1).regimestart)>1
oldstart = regimes0(1).regimestart(end);
end
newstart=1;
if M_.occbin.constraint_nbr==1 && length(regimes_(1).regimestart)>1
newstart = regimes_(1).regimestart(end);
end
if M_.occbin.constraint_nbr==1 && (newstart-oldstart)>2 && options_.occbin.filter.use_relaxation
regimestart = max(oldstart+2,round(0.5*(newstart+oldstart)));
regimestart = min(regimestart,oldstart+4);
if regimestart<=regimes_(1).regimestart(end-1)
if length(regimes_(1).regimestart)<=3
regimestart = max(regimestart, min(regimes_(1).regimestart(end-1)+2,newstart));
else
regimes_(1).regime = regimes_(1).regime(1:end-2);
regimes_(1).regimestart = regimes_(1).regimestart(1:end-2);
regimestart = max(regimestart, regimes_(1).regimestart(end-1)+1);
end
end
regimes_(1).regimestart(end)=regimestart;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],'reduced_state_space', T0, R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
elseif newguess==0
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
end
newguess=0;
regime_hist(niter) = {regimes_(1)};
if M_.occbin.constraint_nbr==1
regime_end(niter) = regimes_(1).regimestart(end);
end
[a, a1, P, P1, v, alphahat, etahat, lik] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
[a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
if error_flag
etahat=NaN(size(QQQ,1),1);
warning(orig_warning_state);
return;
end
etahat_hist(niter) = {etahat};
lik_hist(niter) = lik;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
if not(options_.occbin.filter.use_relaxation)
opts_simul.init_regime=regimes_(1);
end
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
else
......@@ -223,10 +237,16 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul.init_regime=[];
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
warning(orig_warning_state);
return;
end
regimes0=regimes_;
......@@ -235,36 +255,14 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
for kiter=1:niter-1
is_periodic(kiter) = isequal(regime_hist{kiter}, regimes_(1));
end
is_periodic_iter = find(is_periodic);
is_periodic = any(is_periodic);
if is_periodic
if nguess<3 && M_.occbin.constraint_nbr==1
newguess=1;
is_periodic=0;
nguess=nguess+1;
if nguess==1
% change starting regime
regimes_(1).regime=0;
regimes_(1).regimestart=1;
elseif nguess==2
% change starting regime
regimes_(1).regime=[0 1 0];
regimes_(1).regimestart=[1 2 3];
else
regimes_(1).regime=[1 0];
regimes_(1).regimestart=[1 2];
end
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],'reduced_state_space',T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
regime_hist = regime_hist(1);
niter=1;
else
% re-set to previous regime
regimes_ = regimes0;
% force projection conditional on previous regime
opts_simul.init_regime=regimes0(1);
if options_.occbin.filter.periodic_solution
% force projection conditional on most likely regime
[m, im]=min(lik_hist(is_periodic_iter:end));
opts_simul.init_regime=regime_hist{is_periodic_iter+im-1};
if M_.occbin.constraint_nbr==1
myregimestart = [regimes0.regimestart];
else
......@@ -273,79 +271,42 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
opts_simul.maxit=1;
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
warning(orig_warning_state);
return;
end
end
end
end
end
end
error_flag = out.error_flag;
if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1))
error_flag = 1;
if M_.occbin.constraint_nbr==1 % try some other regime
[ll, il]=sort(lik_hist);
[ll, il]=sort(regime_end);
rr=regime_hist(il(2:3));
newstart=1;
if length(rr{1}.regimestart)>1
newstart = rr{1}.regimestart(end)-rr{1}.regimestart(end-1)+1;
end
oldstart=1;
if length(rr{2}.regimestart)>1
oldstart = rr{2}.regimestart(end)-rr{2}.regimestart(end-1)+1;
end
nstart=sort([newstart oldstart]);
regimes_=rr{1}(1);
for k=(nstart(1)+1):(nstart(2)-1)
niter=niter+1;
regimes_(1).regimestart(end)=k;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],'reduced_state_space',T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
[a, a1, P, P1, v, alphahat, etahat, lik] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
etahat_hist(niter) = {etahat};
lik_hist(niter) = lik;
regime_hist(niter) = {regimes_(1)};
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
regimes_ = out.regime_history;
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
[a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
if error_flag
etahat=NaN(size(QQQ,1),1);
warning(orig_warning_state);
return;
end
% opts_simul.init_regime=regimes_(1);
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
else
myregimestart = [regimes_.regimestart1 regimes_.regimestart2];
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
else
error_flag = 330;
etahat=etahat(:,2);
lik=inf;
warning(orig_warning_state);
return;
end
if isequal(out.regime_history(1),regimes_(1))
error_flag=0;
break
end
end
regimes_ = out.regime_history;
end
end
error_flag = out.error_flag;
if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1))
error_flag = 331;
end
if ~error_flag
a = out.piecewise(1:2,my_order_var)' - repmat(out.ys(my_order_var),1,2);
regimes_=regimes_(1:3);
......@@ -356,18 +317,28 @@ C = ss.C(my_order_var,1:2);
QQ = R(:,:,2)*QQQ(:,:,3)*transpose(R(:,:,2));
P(:,:,1) = P(:,:,2);
P(:,:,2) = T(:,:,2)*P(:,:,1)*transpose(T(:,:,2))+QQ;
if nargout>=13
etahat=etahat(:,2);
end
warning_config;
%reset warning state
warning(orig_warning_state);
end
function [a, a1, P, P1, v, alphahat, etahat, lik, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, rescale_prediction_error_covariance, IF_likelihood)
function [a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, rescale_prediction_error_covariance, IF_likelihood, state_uncertainty_flag)
alphahat=NaN(size(a));
etahat=NaN(size(QQQ,1),2);
lik=Inf;
error_flag=0;
warning off
if state_uncertainty_flag
V = zeros(mm,mm,2);
N = zeros(mm,mm,3);
else
V=[];
end
orig_warning_state = warning;
if nargin<18
IF_likelihood=0;
end
......@@ -396,7 +367,8 @@ else
F = ZZ*P(:,:,t)*ZZ' + H(di,di);
sig=sqrt(diag(F));
if any(any(isnan(F)))
error_flag=1;
error_flag=325;
warning(orig_warning_state);
return;
end
if rank(F)<size(F,1)
......@@ -442,9 +414,15 @@ while t > 1
if isempty(di)
% in this case, L is simply T due to Z=0, so that DK (2012), eq. 4.93 obtains
r(:,t) = L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38 with Z=0
if state_uncertainty_flag
N(:,:,t)=L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42 with Z=0
end
else
ZZ = Z(di,:);
r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38
if state_uncertainty_flag
N(:,:,t)=ZZ'*iF(di,di,t)*ZZ+L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42
end
end
Q=QQQ(:,:,t);
QRt = Q*transpose(RR(:,:,t));
......@@ -452,6 +430,10 @@ while t > 1
alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t); %DK (2012), eq. 4.35
etahat(:,t) = QRt*r(:,t); %DK (2012), eq. 4.63
r(:,t) = T'*r(:,t); % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
if state_uncertainty_flag
V(:,:,t) = P1(:,:,t)-P1(:,:,t)*N(:,:,t)*P1(:,:,t); %DK (2012), eq. 4.43
N(:,:,t) = T'*N(:,:,t)*T; %DK (2012), eq. 4.43
end
if IF_likelihood && t==2 && not(isempty(di))
ishocks = any(ZZ*RR(:,:,t));
......@@ -469,5 +451,6 @@ while t > 1
end
end
warning_config;
%reset warning state
warning(orig_warning_state);
end
matlab/+occbin/kalman_update_algo_3.m
View file @ 4a32acc2
function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,oo_,options_,occbin_options,kalman_tol,nk)
% function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,oo_,options_,occbin_options,kalman_tol,nk)
function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, lik, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options,kalman_tol,nk)
% function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options,kalman_tol,nk)
%
% INPUTS
% - a [N by 1] t-1's state estimate
......@@ -20,10 +20,13 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat,
% - RR [N by N_exo by 2] shock impact matrix at t-1:t
% - CC [N by 2] state space constant state transition matrix at t-1:t
% - regimes0 [structure] regime info at t-1:t
% - M_ [structure] Matlab's structure describing the model (M_).
% - options_ [structure] Matlab's structure describing the current options (options_).
% - oo_ [structure] Matlab's structure containing the results (oo_).
% - occbin_options_ [structure] Matlab's structure describing the Occbin options.
% - M_ [structure] MATLAB's structure describing the model (M_).
% - options_ [structure] MATLAB's structure describing the current options (options_).
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - occbin_options_ [structure] MATLAB's structure describing the OccBin options.
% - kalman_tol [double] tolerance for reciprocal condition number
% - nk [double] number of forecasting periods
%
......@@ -40,7 +43,7 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat,
% - CC [N by 2] state space constant state transition matrix at t-1:t
% - regimes_ [structure] regime info at t-1:t
% - error_flag [structure] error flag
% - M_ [structure] Matlab's structure describing the model (M_).
% - M_ [structure] MATLAB's structure describing the model (M_).
% - alphahat: smoothed variables (a_{t|T})
% - etahat: smoothed shocks
% - TT [N by N by 2] state transition matrix at t-1:t
......@@ -52,7 +55,7 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat,
% constraints, Working Papers 2021-03, Joint Research Centre, European Commission
% Copyright © 2021 Dynare Team
% Copyright © 2021-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -82,7 +85,8 @@ if isempty(nk)
end
nk=max(nk,1);
opts_simul = occbin_options.opts_regime;
opts_simul = occbin_options.opts_simul;
options_.occbin.simul=opts_simul;
base_regime = struct();
if M_.occbin.constraint_nbr==1
base_regime.regime = 0;
......@@ -93,6 +97,7 @@ else
base_regime.regime2 = 0;
base_regime.regimestart2 = 1;
end
regimes_ = [base_regime base_regime base_regime];
myrestrict=[];
if options_.smoother_redux
opts_simul.restrict_state_space =1;
......@@ -100,15 +105,46 @@ if options_.smoother_redux
end
mm=size(a,1);
if ~options_.occbin.filter.use_relaxation
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
if ~isempty(fieldnames(regimes0))
if options_.occbin.filter.guess_regime
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, regimes0(1),myrestrict,T0,R0);
if M_.occbin.constraint_nbr==1
bindx = occbin.backward_map_regime(regimes0(1).regime, regimes0(1).regimestart);
bindx = bindx(2:end);
[regimes0(2).regime, regimes0(2).regimestart, error_flag]=occbin.map_regime(bindx);
bindx = bindx(2:end);
[regimes0(3).regime, regimes0(3).regimestart, error_flag]=occbin.map_regime(bindx);
else
bindx1 = occbin.backward_map_regime(regimes0(1).regime1, regimes0(1).regimestart1);
bindx2 = occbin.backward_map_regime(regimes0(1).regime2, regimes0(1).regimestart2);
bindx1 = bindx1(2:end);
bindx2 = bindx2(2:end);
[regimes0(2).regime1, regimes0(2).regimestart1, error_flag]=occbin.map_regime(bindx1);
[regimes0(2).regime2, regimes0(2).regimestart2, error_flag]=occbin.map_regime(bindx2);
bindx1 = bindx1(2:end);
bindx2 = bindx2(2:end);
[regimes0(3).regime1, regimes0(3).regimestart1, error_flag]=occbin.map_regime(bindx1);
[regimes0(3).regime2, regimes0(3).regimestart2, error_flag]=occbin.map_regime(bindx2);
end
% regimes0=[regimes0 base_regime base_regime];
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
end
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
else
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, base_regime,myrestrict,T0,R0);
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, base_regime,myrestrict,T0,R0);
if isempty(fieldnames(regimes0))
regimes0 = regimes_;
else
regimes0(1)=base_regime;
end
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
regimes0(1)=base_regime;
end
......@@ -120,18 +156,32 @@ opts_simul.exo_pos=1:M_.exo_nbr;
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
my_order_var = oo_.dr.order_var(oo_.dr.restrict_var_list);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
my_order_var = dr.order_var(dr.restrict_var_list);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
my_order_var = oo_.dr.order_var;
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
my_order_var = dr.order_var;
end
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if options_.occbin.filter.guess_regime
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
else
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
end
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
return;
end
......@@ -145,76 +195,31 @@ regime_hist = {regimes0(1)};
if M_.occbin.constraint_nbr==1
regime_end = regimes0(1).regimestart(end);
end
lik_hist=lik;
niter=1;
is_periodic=0;
if options_.occbin.filter.use_relaxation || isequal(regimes0(1),base_regime)
nguess=1;
else
nguess=0;
end
newguess=0;
if any(myregime) || ~isequal(regimes_(1),regimes0(1))
while ~isequal(regimes_(1),regimes0(1)) && ~is_periodic && ~out.error_flag && niter<=options_.occbin.likelihood.max_number_of_iterations
niter=niter+1;
newstart=1;
if M_.occbin.constraint_nbr==1 && length(regimes_(1).regimestart)>1
newstart = regimes_(1).regimestart(end);
end
oldstart=1;
if M_.occbin.constraint_nbr==1 && length(regimes0(1).regimestart)>1
oldstart = regimes0(1).regimestart(end);
end
if M_.occbin.constraint_nbr==1 && (newstart-oldstart)>2 && options_.occbin.filter.use_relaxation
regimestart = max(oldstart+2,round(0.5*(newstart+oldstart)));
regimestart = min(regimestart,oldstart+4);
if regimestart<=regimes_(1).regimestart(end-1)
if length(regimes_(1).regimestart)<=3
regimestart = max(regimestart, min(regimes_(1).regimestart(end-1)+2,newstart));
else
regimes_(1).regime = regimes_(1).regime(1:end-2);
regimes_(1).regimestart = regimes_(1).regimestart(1:end-2);
regimestart = max(regimestart, regimes_(1).regimestart(end-1)+1);
end
end
% % if (newstart-oldstart)>3
% % regimestart = regimes_(1).regimestart(end-1)+oldstart+2;
% % % regimestart = regimes_(1).regimestart(end-1)+round(0.5*(newstart+oldstart))-1;
regimes_(1).regimestart(end)=regimestart;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
elseif newguess==0
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
end
newguess=0;
regime_hist(niter) = {regimes_(1)};
if M_.occbin.constraint_nbr==1
regime_end(niter) = regimes_(1).regimestart(end);
end
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
lik_hist(niter) = lik;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
% if opts_simul.restrict_state_space
% tmp=zeros(M_.endo_nbr,1);
% tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
% opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
% else
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
% end
if not(options_.occbin.filter.use_relaxation)
opts_simul.init_regime=regimes_(1);
end
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
else
......@@ -222,9 +227,11 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
return;
end
regimes0=regimes_;
......@@ -233,36 +240,14 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
for kiter=1:niter-1
is_periodic(kiter) = isequal(regime_hist{kiter}, regimes_(1));
end
is_periodic_iter = find(is_periodic);
is_periodic = any(is_periodic);
if is_periodic
if nguess<3 && M_.occbin.constraint_nbr==1
newguess=1;
is_periodic=0;
nguess=nguess+1;
if nguess==1
% change starting regime
regimes_(1).regime=0;
regimes_(1).regimestart=1;
elseif nguess==2
% change starting regime
regimes_(1).regime=[0 1 0];
regimes_(1).regimestart=[1 2 3];
else
regimes_(1).regime=[1 0];
regimes_(1).regimestart=[1 2];
end
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
regime_hist = regime_hist(1);
niter=1;
else
% re-set to previous regime
regimes_ = regimes0;
% force projection conditional on previous regime
opts_simul.init_regime=regimes0(1);
if options_.occbin.filter.periodic_solution
% force projection conditional on most likely regime
[m, im]=min(lik_hist(is_periodic_iter:end));
opts_simul.init_regime=regime_hist{is_periodic_iter+im-1};
if M_.occbin.constraint_nbr==1
myregimestart = [regimes0.regimestart];
else
......@@ -271,74 +256,45 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
opts_simul.maxit=1;
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
end
end
end
end
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
if error_flag==0 && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1)) %fixed point algorithm did not converge
error_flag = 1;
if M_.occbin.constraint_nbr==1
% try some other regime before giving up
[ll, il]=sort(regime_end);
rr=regime_hist(il(2:3));
newstart=1;
if length(rr{1}(1).regimestart)>1
newstart = rr{1}(1).regimestart(end)-rr{1}(1).regimestart(end-1)+1;
end
oldstart=1;
if length(rr{2}(1).regimestart)>1
oldstart = rr{2}(1).regimestart(end)-rr{2}(1).regimestart(end-1)+1;
end
nstart=sort([newstart oldstart]);
regimes_=rr{1}(1);
for k=(nstart(1)+1):(nstart(2)-1)
niter=niter+1;
regimes_(1).regimestart(end)=k;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
etahat=etahat(:,2);
lik=inf;
return;
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
regimes_ = out.regime_history;
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
end
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
else
myregimestart = [regimes_.regimestart1 regimes_.regimestart2];
error_flag = 330;
etahat=etahat(:,2);
lik=inf;
return;
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_);
if isequal(out.regime_history(1),regimes_(1))
error_flag=0;
break
end
end
regimes_ = out.regime_history;
end
end
regimes_=regimes_(1:3);
error_flag = out.error_flag;
if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1)) %fixed point algorithm did not converge
error_flag = 331;
end
if ~error_flag
a = out.piecewise(1:nk+1,my_order_var)' - repmat(out.ys(my_order_var),1,nk+1);
regimes_=regimes_(1:3);
end
T = ss.T(my_order_var,my_order_var,:);
R = ss.R(my_order_var,:,:);
C = ss.C(my_order_var,:);
TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,1);
RR = ss.R(oo_.dr.order_var,:,1);
CC = ss.C(oo_.dr.order_var,1);
TT = ss.T(dr.order_var,dr.order_var,1);
RR = ss.R(dr.order_var,:,1);
CC = ss.C(dr.order_var,1);
QQ = R(:,:,2)*QQQ(:,:,3)*transpose(R(:,:,2));
P(:,:,1) = P(:,:,2);
for j=1:nk
......@@ -349,7 +305,7 @@ end
end
function [a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol)
function [a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol)
% [a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol)
% - a
% - a1
......@@ -374,9 +330,18 @@ a(:,t) = a1(:,t);
P1(:,:,t) = T*P(:,:,t-1)*T' + QQ; %transition according to (6.14) in DK (2012)
P(:,:,t) = P1(:,:,t);
di = data_index{t}';
% store info for lik
if not(isempty(di))
vv = Y(di,t) - Z(di,:)*a(:,t);
F = Z(di,:)*P(:,:,t)*Z(di,:)' + diag(H(di));
sig=sqrt(diag(F));
lik=inf;
end
if isempty(di)
Fi(:,t) = 0;
Ki(:,:,t) = 0;
lik =0;
end
for i=di
Zi = Z(i,:);
......@@ -391,6 +356,11 @@ for i=di
% p. 157, DK (2012)
end
end
if not(isempty(di))
log_dF = log(det(F./(sig*sig')))+2*sum(log(sig));
iF = inv(F./(sig*sig'))./(sig*sig');
lik = log_dF + transpose(vv)*iF*vv + length(di)*log(2*pi);
end
%% do backward pass
ri=zeros(mm,1);
......
matlab/+occbin/kalman_update_engine.m 0 → 100644
View file @ 4a32acc2
function [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V, Fix, Kix, TTx,RRx,CCx] = ...
kalman_update_engine(a0,a1,P0,P1,t,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_,base_regime,d_index,M_,...
dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options, Fi,Ki,kalman_tol,nk)
% [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V, Fix, Kix, TTx,RRx,CCx] = kalman_update_engine(
% a0,a1,P0,P1,t,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_,base_regime,d_index,M_,
% dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options, Fi,Ki,kalman_tol,nk)
% Copyright © 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
is_multivariate = true;
Fix=[];
Kix=[];
TTx=[];
RRx=[];
CCx=[];
V=[];
if nargin>26
is_multivariate = false;
end
if is_multivariate
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,struct(),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, regx, info, M_, likx, alphahat, etahat,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,struct(),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
likvec = likx;
regvec = regx(1);
info0=info;
if info
if ~isequal(regimes_(1:2),[base_regime base_regime])
if is_multivariate
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, regx, info, M_, likx, alphahat, etahat,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
end
likvec = likx;
regvec = regx(1);
info1=info;
else
if ~isequal(regimes_(1:2),[base_regime base_regime])
if is_multivariate
[ax1, a1x1, Px1, P1x1, vx1, Tx1, Rx1, Cx1, regx1, info1, M_1, likx1, etahat1, alphahat1, V1] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax1, a1x1, Px1, P1x1, vx1, Fix1, Kix1, Tx1, Rx1, Cx1, regx1, info1, M_1, likx1, alphahat1, etahat1,TTx1,RRx1,CCx1] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info1==0 && not(isequal(regx1,regx))
likvec = [likvec likx1];
regvec = [regvec; regx1(1)];
end
if info1==0 && likx1<likx
ax=ax1;
a1x=a1x1;
Px=Px1;
P1x=P1x1;
vx=vx1;
Tx=Tx1;
Rx=Rx1;
Cx=Cx1;
regx=regx1;
info=info1;
M_= M_1;
likx=likx1;
etahat=etahat1;
alphahat=alphahat1;
if is_multivariate
V=V1;
else
Fix = Fix1;
Kix = Kix1;
TTx = TTx1;
RRx = RRx1;
CCx = CCx1;
end
end
else
if t>options_.occbin.likelihood.number_of_initial_periods_with_extra_regime_guess
info1=0;
else
% may help in first 2 periods to try some other guess regime, due to
% larger state uncertainty
info1=1;
options_.occbin.likelihood.brute_force_regime_guess = true;
options_.occbin.likelihood.loss_function_regime_guess = true;
end
end
end
diffstart=0;
if info==0
if M_.occbin.constraint_nbr==1
oldstart = regimes_(1).regimestart(end);
newstart = regx(1).regimestart(end);
diffstart = newstart-oldstart;
regname = 'regimestart';
reg_string = 'regime';
else
newstart1 = regx(1).regimestart1(end);
newstart2 = regx(1).regimestart2(end);
oldstart1 = regimes_(1).regimestart1(end);
oldstart2 = regimes_(1).regimestart2(end);
[diffstart, diffregime] = max([newstart1-oldstart1,newstart2-oldstart2]);
switch diffregime
case 1
regname = 'regimestart1';
reg_string = 'regime1';
case 2
regname = 'regimestart2';
reg_string = 'regime2';
end
end
end
if options_.occbin.filter.use_relaxation && diffstart>options_.occbin.filter.use_relaxation
guess_regime = [base_regime base_regime];
options_.occbin.filter.guess_regime = true;
guess_regime(1) = regx(1);
regx2= regx;
while isequal(guess_regime(1),regx2(1))
% we reduce length until the converged regime does not change
guess_regime(1).(regname)(end) = regx2(1).(regname)(end)-1;
if guess_regime(1).(regname)(end)==1
% make sure we enforce base regime
guess_regime(1).(regname)=guess_regime(1).(regname)(end);
guess_regime(1).(reg_string)=0;
else
if guess_regime(1).(regname)(end-1)==guess_regime(1).(regname)(end)
guess_regime(1).(regname)(end-1) = guess_regime(1).(regname)(end-1)-1;
end
end
if is_multivariate
[ax2, a1x2, Px2, P1x2, vx2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, etahat2, alphahat2, V2] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2, a1x2, Px2, P1x2, vx2, Fix2, Kix2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, alphahat2, etahat2,TTx2,RRx2,CCx2] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
isnew=true;
for kr=1:length(likvec)
% make sure likelihood does not differ by rounding issue
% but due to different regimes
if isequal(regx2(1),regvec(kr))
isnew = false;
end
end
if isnew
likvec = [likvec likx2];
regvec = [regvec; regx2(1)];
end
if info2==0 && likx2<likx
ax=ax2;
a1x=a1x2;
Px=Px2;
P1x=P1x2;
vx=vx2;
Tx=Tx2;
Rx=Rx2;
Cx=Cx2;
regx=regx2;
info=info2;
likx=likx2;
M_= M_2;
etahat=etahat2;
alphahat=alphahat2;
if is_multivariate
V=V2;
else
Fix = Fix2;
Kix = Kix2;
TTx = TTx2;
RRx = RRx2;
CCx = CCx2;
end
end
end
options_.occbin.filter.guess_regime = false;
options_.occbin.likelihood.brute_force_regime_guess = false;
options_.occbin.likelihood.loss_function_regime_guess = false;
end
if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (info==0 && ~isequal(regx(1),base_regime))
guess_regime = [base_regime base_regime];
options_.occbin.filter.guess_regime = true;
use_index = 0;
if M_.occbin.constraint_nbr==1
for k=1:5
guess_regime(1).regimestart=[1 5 5+4*k];
guess_regime(1).regime=[0 1 0];
if is_multivariate
[ax2{1}, a1x2{1}, Px2{1}, P1x2{1}, vx2{1}, Tx2{1}, Rx2{1}, Cx2{1}, regx2{1}, info2, M_2{1}, likx2{1}, etahat2{1}, alphahat2{1}, V2{1}] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2{1}, a1x2{1}, Px2{1}, P1x2{1}, vx2{1}, Fix2{1}, Kix2{1}, Tx2{1}, Rx2{1}, Cx2{1}, regx2{1}, info2, M_2{1}, likx2{1}, alphahat2{1}, etahat2{1},TTx2{1},RRx2{1},CCx2{1}] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info2==0
use_index= 1;
if not(info==0 && isequal(regx2{1},regx))
% found a solution, different from previous
break
end
end
guess_regime(1).regimestart=[1 1+4*k];
guess_regime(1).regime=[1 0];
if is_multivariate
[ax2{2}, a1x2{2}, Px2{2}, P1x2{2}, vx2{2}, Tx2{2}, Rx2{2}, Cx2{2}, regx2{2}, info2, M_2{2}, likx2{2}, etahat2{2}, alphahat2{2}, V2{2}] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2{2}, a1x2{2}, Px2{2}, P1x2{2}, vx2{2}, Fix2{2}, Kix2{2}, Tx2{2}, Rx2{2}, Cx2{2}, regx2{2}, info2, M_2{2}, likx2{2}, alphahat2{2}, etahat2{2},TTx2{2},RRx2{2},CCx2{2}] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info2==0
use_index = 2;
end
if use_index
% in case the second guess does not find a solution!
info2=0;
% a solution was found
break
end
end
end
if M_.occbin.constraint_nbr==2
for jk=0:1 % loop over other regime duration. this loop is shorter for parsimony. one may add an option ...
for k=1:5 % loop over current regime duration
gindex = 0;
for jr=1:2 % loop over current regime 1 or 2
if jr==1
regstart1 = 'regimestart1';
reg1 = 'regime1';
regstart2 = 'regimestart2';
reg2 = 'regime2';
else
regstart1 = 'regimestart2';
reg1 = 'regime2';
regstart2 = 'regimestart1';
reg2 = 'regime1';
end
for kk=1:2 % loop over current regime binding in expectation vs binding in current period
if kk==1
guess_regime(1).(regstart1)=[1 5 5+4*k];
guess_regime(1).(reg1)=[0 1 0];
else
guess_regime(1).(regstart1)=[1 1+4*k];
guess_regime(1).(reg1)=[1 0];
end
for kj=1:1+1*(jk>0)
% loop over other regime slack or binding in current period or binding in
% expectation
if jk==0
% other regime is slack
guess_regime(1).(regstart2) = 1;
guess_regime(1).(reg2) = 0;
else % jk>0
if kj==1
% other regime binding in current period
guess_regime(1).(regstart2)=[1 1+4*jk];
guess_regime(1).(reg2) = [1 0];
else
% other regime binding in expectation
guess_regime(1).(regstart2)=[1 5 5+4*jk];
guess_regime(1).(reg2) = [0 1 0];
end
end
gindex = gindex+1;
if is_multivariate
[ax2{gindex}, a1x2{gindex}, Px2{gindex}, P1x2{gindex}, vx2{gindex}, Tx2{gindex}, Rx2{gindex}, Cx2{gindex}, regx2{gindex}, info2, M_2{gindex}, likx2{gindex}, etahat2{gindex}, alphahat2{gindex}, V2{gindex}] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2{gindex}, a1x2{gindex}, Px2{gindex}, P1x2{gindex}, vx2{gindex}, Fix2{gindex}, Kix2{gindex}, Tx2{gindex}, Rx2{gindex}, Cx2{gindex}, regx2{gindex}, info2, M_2{gindex}, likx2{gindex}, alphahat2{gindex}, etahat2{gindex},TTx2{gindex},RRx2{gindex},CCx2{gindex}] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info2==0
use_index= gindex;
if not(info==0 && isequal(regx2{gindex},regx))
% found a solution, different from previous one
break
end
end
end % loop over other regime slack, binding in expectation or binding in current period
if info2==0
if not(info==0 && isequal(regx2{gindex},regx))
% found a solution, different from previous one
break
end
end
end % loop over current regime binding in expectation vs binding in current period
if info2==0
if not(info==0 && isequal(regx2{gindex},regx))
% found a solution, different from previous one
break
end
end
end % loop over current regime 1 or 2
if use_index
info2=0;
break
end
end % loop over current regime duration
if use_index
break
end
end % loop over other regime duration
end % 2 constraints
if info2==0
% so that we DO NOT enter IVF step
info0=0;
info1=0;
isnew=true;
for kr=1:length(likvec)
% make sure likelihood does not differ by rounding issue
% but truly for different regimes
if isequal(regx2{use_index}(1),regvec(kr))
isnew = false;
end
end
if isnew
likvec = [likvec likx2{use_index}];
regvec = [regvec; regx2{use_index}(1)];
end
end
if info2==0 && likx2{use_index}<likx
ax=ax2{use_index};
a1x=a1x2{use_index};
Px=Px2{use_index};
P1x=P1x2{use_index};
vx=vx2{use_index};
Tx=Tx2{use_index};
Rx=Rx2{use_index};
Cx=Cx2{use_index};
regx=regx2{use_index};
info=info2;
M_= M_2{use_index};
likx=likx2{use_index};
etahat=etahat2{use_index};
alphahat=alphahat2{use_index};
if is_multivariate
V=V2{use_index};
else
Fix = Fix2{use_index};
Kix = Kix2{use_index};
TTx = TTx2{use_index};
RRx = RRx2{use_index};
CCx = CCx2{use_index};
end
end
options_.occbin.filter.guess_regime = false;
end
if options_.occbin.likelihood.loss_function_regime_guess && (info0 || info1) %|| (info==0 && ~isequal(regx(1),base_regime))
[~, out] = occbin.findmin(d_index, a0, P1, Qt, Y, Z, occbin_options.opts_simul,M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_);
if out.error_flag==0
options_.occbin.filter.guess_regime = true;
guess_regime=out.regime_history;
guess_regime = [guess_regime base_regime];
if is_multivariate
[ax2, a1x2, Px2, P1x2, vx2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, etahat2, alphahat2, V2] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2, a1x2, Px2, P1x2, vx2, Fix2, Kix2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, alphahat2, etahat2,TTx2,RRx2,CCx2] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
options_.occbin.filter.guess_regime = false;
isnew=true;
for kr=1:length(likvec)
% make sure likelihood does not differ by rounding issue
% but due to different regimes
if isequal(regx2(1),regvec(kr))
isnew = false;
end
end
if isnew
likvec = [likvec likx2];
regvec = [regvec; regx2(1)];
end
if info2==0 && likx2<likx
ax=ax2;
a1x=a1x2;
Px=Px2;
P1x=P1x2;
vx=vx2;
Tx=Tx2;
Rx=Rx2;
Cx=Cx2;
regx=regx2;
info=info2;
likx=likx2;
M_= M_2;
etahat=etahat2;
alphahat=alphahat2;
if is_multivariate
V=V2;
else
Fix = Fix2;
Kix = Kix2;
TTx = TTx2;
RRx = RRx2;
CCx = CCx2;
end
end
end
end
if length(likvec)>1
% sum the likelihood of multiple solutions
likx = -2*log(sum(exp(-likvec./2)));
end
if info(1)==0
if isnan(likx)
info = 323;
elseif any(any(isnan(ax)))
info = 324;
end
end
end