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
Loading items

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
Loading items
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 499 additions and 286 deletions
macOS/welcome_template.html
View file @ 2ddfa086
<html lang="en"> <html lang="en">
<body> <head>
<h3 style="text-align: center;">Welcome to Dynare</h3> <meta charset="UTF-8">
<p style="text-align: center;">Version VERSION_NO_SPACE</p> <meta name="viewport" content="width=device-width, initial-scale=1">
<p style="text-align: center;">DATE</p> <title>Dynare Installer</title>
<style>
<p><b>Just a few things to note</b>. This installation can be customized as you can choose not to install the GNU Compiler Collection (GCC). Installing GCC is necessary if you want to use the <tt>use_dll</tt> option in Dynare, but otherwise unnecessary.</p> .center-text {
text-align: center;
}
.bold-text {
font-weight: bold;
}
code {
font-family: "Courier New", Courier, monospace;
}
</style>
</head>
<p>To install GCC we run a script that first installs the XCode Command Line Tools (an Apple product). The script then installs Homebrew, a package manager for macOS and, finally, GCC itself. Both Homebrew and GCC will be installed in your Dynare installation folder. So, when you delete this folder, they too will be deleted.</p> <body>
<header class="center-text">
<h1>Welcome to Dynare</h1>
<p>Version VERSION_NO_SPACE</p>
<p>DATE</p>
</header>
<p>Installing GCC will require an active internet connection with the ability to connect to the Apple servers and GitHub. The installation will take anywhere from a few minutes to a half an hour during the <i>Running package scripts</i> phase of Installation. The time it takes depends on your internet speed, the speed of your computer, and whether or not you already have XCode Command Line Tools installed. The progress bar will not advance during this phase. Please be patient.</p> <main>
<p>Thank you for choosing Dynare!</p>
<p>This installation does not require administrative privileges. If for some reason admin rights are requested, use 'Change Install Location' and select 'Install for me only'.</p>
<p>By default Dynare will be installed into /Applications/Dynare. To modify the installation path, click <strong>Customize</strong> after accepting the license. Then, under <strong>Location</strong>, select your desired folder.</p>
<p>Installing into /Applications/Dynare might fail if you have older versions of Dynare installed in /Applications/Dynare. To fix this, modify the ownership by executing the following command in Terminal.app:</p>
<code>
sudo chown -R "$USER":staff /Applications/Dynare
</code>
</main>
<p> You can choose not to install GCC by choosing <i>Customize</i> from the <i>Installation Type</i> screen and deselecting <i>GCC compiler</i>. If you already have <tt>GCC_BINARY</tt> installed under <tt>/usr/local</tt>, you can forgo the installation of GCC here as Dynare will find your system compiler when you use <tt>use_dll</tt>.</p>
</body> </body>
</html> </html>
\ No newline at end of file
matlab/+estimate/nls.m
View file @ 2ddfa086
...@@ -136,10 +136,27 @@ else ...@@ -136,10 +136,27 @@ else
end end
end end
% Remove residuals from the equation. Note that a plus or minus will % Remove residuals from the equation. Note that a plus or minus will remain in the equation
% remain in the equation, but this seems to be without consequence.
rhs = regexprep(rhs, rname, ''); rhs = regexprep(rhs, rname, '');
% FIXME The JSON output for rhs (with aux variables substitutions) is not always
% the same regarding to the position of the residual. If the residual appears at
% the end of the equation, after the removal of rname the rhs will end with a +
% symbol which will result in a crash later when evaluating the sum of square
% residuals. If the residual appears at the begining of the equation, after the
% removal of rname the rhs will begin with a + symbol which is just awful (but
% will not cause any trouble).
% Remove trailing + (if any, introduced when removing residual)
if isequal(rhs(end), '+')
rhs = rhs(1:end-1);
end
% Remove leading + (if any, introduced when removing residual)
if isequal(rhs(1), '+')
rhs = rhs(2:end);
end
% %
% Rewrite and print the equation. % Rewrite and print the equation.
% %
...@@ -234,6 +251,13 @@ fprintf(fid, 'r = %s-(%s);\n', lhs, rhs); ...@@ -234,6 +251,13 @@ fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);
fprintf(fid, 's = r''*r;\n'); fprintf(fid, 's = r''*r;\n');
fclose(fid); fclose(fid);
% Workaround for Octave bug https://savannah.gnu.org/bugs/?46282
% Octave will randomly fail to read the ssr_* file generated in the +folder
if isoctave
rename(['+' M_.fname], ['+' M_.fname '-tmp']);
rename(['+' M_.fname '-tmp'], ['+' M_.fname]);
end
% Create a function handle returning the sum of square residuals for a given vector of parameters. % Create a function handle returning the sum of square residuals for a given vector of parameters.
ssrfun = @(p) feval([M_.fname '.ssr_' eqname], p, DATA, M_, oo_); ssrfun = @(p) feval([M_.fname '.ssr_' eqname], p, DATA, M_, oo_);
......
matlab/method_of_moments/method_of_moments_check_plot.m → matlab/+mom/check_plot.m
View file @ 2ddfa086
function method_of_moments_check_plot(fun,xparam,SE_vec,options_,M_,estim_params_,Bounds,bayestopt_,varargin) function check_plot(fun,xparam,SE_vec,options_,M_,estim_params_,Bounds,bayestopt_,varargin)
%function check_plot(fun,xparam,SE_vec,options_,M_,estim_params_,Bounds,bayestopt_,varargin)
% Checks the estimated local minimum of the moment's distance objective % Checks the estimated local minimum of the moment's distance objective
...@@ -44,7 +45,7 @@ if ~exist([M_.dname filesep 'graphs'],'dir') ...@@ -44,7 +45,7 @@ if ~exist([M_.dname filesep 'graphs'],'dir')
end end
if TeX && any(strcmp('eps',cellstr(options_.graph_format))) if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([M_.dname, '/graphs/', M_.fname '_MoMCheckPlots.tex'],'w'); fidTeX = fopen([M_.dname, '/graphs/', M_.fname '_MoMCheckPlots.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by method_of_moments_check_plot.m (Dynare).\n'); fprintf(fidTeX,'%% TeX eps-loader file generated by mom.check_plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']); fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n'); fprintf(fidTeX,' \n');
end end
...@@ -119,7 +120,7 @@ for plt = 1:nbplt ...@@ -119,7 +120,7 @@ for plt = 1:nbplt
else else
y(i,1) = NaN; y(i,1) = NaN;
if options_.debug if options_.debug
fprintf('method_of_moments_check_plot:: could not solve model for parameter %s at value %4.3f, error code: %u\n',name,z(i),info(1)) fprintf('mom.check_plot:: could not solve model for parameter %s at value %4.3f, error code: %u\n',name,z(i),info(1))
end end
end end
if options_.mom.penalized_estimator if options_.mom.penalized_estimator
......
matlab/method_of_moments/method_of_moments_data_moments.m → matlab/+mom/data_moments.m
View file @ 2ddfa086
function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, matched_moments_, options_mom_) function [dataMoments, m_data] = data_moments(data, oo_, matched_moments_, options_mom_)
% [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, matched_moments_, options_mom_) % [dataMoments, m_data] = data_moments(data, oo_, matched_moments_, options_mom_)
% This function computes the user-selected empirical moments from data % This function computes the user-selected empirical moments from data
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
...@@ -13,8 +13,8 @@ function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, match ...@@ -13,8 +13,8 @@ function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, match
% o m_data [T x numMom] selected empirical moments at each point in time % o m_data [T x numMom] selected empirical moments at each point in time
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% o method_of_moments.m % o mom.run.m
% o method_of_moments_objective_function.m % o mom.objective_function.m
% ========================================================================= % =========================================================================
% Copyright (C) 2020-2021 Dynare Team % Copyright (C) 2020-2021 Dynare Team
% %
...@@ -70,6 +70,3 @@ end ...@@ -70,6 +70,3 @@ end
end %function end end %function end
\ No newline at end of file
matlab/method_of_moments/method_of_moments_objective_function.m → matlab/+mom/objective_function.m
View file @ 2ddfa086
function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_) function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
% [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_) % [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function evaluates the objective function for GMM/SMM estimation % This function evaluates the objective function for GMM/SMM estimation
% ========================================================================= % =========================================================================
...@@ -20,11 +20,14 @@ function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_m ...@@ -20,11 +20,14 @@ function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_m
% o oo_: structure containing the results with the following updated fields: % o oo_: structure containing the results with the following updated fields:
% - mom.model_moments [numMom x 1] vector with model moments % - mom.model_moments [numMom x 1] vector with model moments
% - mom.Q value of the quadratic form of the moment difference % - mom.Q value of the quadratic form of the moment difference
% - mom.model_moments_params_derivs
% [numMom x numParams] Jacobian matrix of derivatives of model_moments with respect to estimated parameters
% (only for GMM with analytical derivatives)
% o M_: Matlab's structure describing the model % o M_: Matlab's structure describing the model
% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_) % o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% o method_of_moments.m % o mom.run.m
% o dynare_minimize_objective.m % o dynare_minimize_objective.m
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls % This function calls
...@@ -66,7 +69,7 @@ if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian ...@@ -66,7 +69,7 @@ if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
df = nan(size(oo_.mom.data_moments,1),length(xparam1)); df = nan(size(oo_.mom.data_moments,1),length(xparam1));
end end
else else
df = nan(1,length(xparam1)); df = nan(length(xparam1),1);
end end
else else
df=[]; %required to be empty by e.g. newrat df=[]; %required to be empty by e.g. newrat
...@@ -154,64 +157,54 @@ if strcmp(options_mom_.mom.mom_method,'GMM') ...@@ -154,64 +157,54 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
end end
oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1); oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
offset = 0; for jm = 1:size(M_.matched_moments,1)
% First moments % First moments
if ~options_mom_.prefilter && isfield(options_mom_.mom.index,'E_y') && nnz(options_mom_.mom.index.E_y) > 0 if ~options_mom_.prefilter && (sum(M_.matched_moments{jm,3}) == 1)
E_y = pruned_state_space.E_y; idx1 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}) );
E_y_nbr = nnz(options_mom_.mom.index.E_y); oo_.mom.model_moments(jm,1) = pruned_state_space.E_y(idx1);
oo_.mom.model_moments(offset+1:E_y_nbr,1) = E_y(options_mom_.mom.index.E_y);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(offset+1:E_y_nbr,:) = pruned_state_space.dE_y(options_mom_.mom.index.E_y,:); oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
end end
offset = offset + E_y_nbr;
end end
% Second moments % Second moments
% Contemporaneous covariance if (sum(M_.matched_moments{jm,3}) == 2)
if isfield(options_mom_.mom.index,'E_yy') && nnz(options_mom_.mom.index.E_yy) > 0 idx1 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(1)) );
idx2 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(2)) );
if nnz(M_.matched_moments{jm,2}) == 0
% Covariance
if options_mom_.prefilter if options_mom_.prefilter
E_yy = pruned_state_space.Var_y; oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
dE_yy = pruned_state_space.dVar_y; oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_y(idx1,idx2,:);
end end
else else
E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y'; oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
dE_yy = pruned_state_space.dVar_y;
for jp=1:totparam_nbr for jp=1:totparam_nbr
dE_yy(:,:,jp) = dE_yy(:,:,jp) + pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y' + pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)'; oo_.mom.model_moments_params_derivs(jm,jp) = pruned_state_space.dVar_y(idx1,idx2,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)';
end end
end end
end end
E_yy_nbr = nnz(tril(options_mom_.mom.index.E_yy)); else
oo_.mom.model_moments(offset+(1:E_yy_nbr),1) = E_yy(tril(options_mom_.mom.index.E_yy)); % Autocovariance
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) lag = -M_.matched_moments{jm,2}(2); %note that leads/lags in matched_moments are transformed such that first entry is always 0 and the second is a lag
oo_.mom.model_moments_params_derivs(offset+(1:E_yy_nbr),:) = reshape(dE_yy(repmat(tril(options_mom_.mom.index.E_yy),[1 1 totparam_nbr])),E_yy_nbr,totparam_nbr);
end
offset = offset + E_yy_nbr;
end
% Lead/lags covariance
if isfield(options_mom_.mom.index,'E_yyt') && nnz(options_mom_.mom.index.E_yyt) > 0
if options_mom_.prefilter if options_mom_.prefilter
E_yyt = pruned_state_space.Var_yi; oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
dE_yyt = pruned_state_space.dVar_yi; oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_yi(idx1,idx2,lag,:);
end end
else else
E_yyt = pruned_state_space.Var_yi + repmat(pruned_state_space.E_y*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)]); oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
dE_yyt = pruned_state_space.dVar_yi;
for jp=1:totparam_nbr for jp=1:totparam_nbr
dE_yyt(:,:,:,jp) = dE_yyt(:,:,:,jp) + repmat(pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)])... oo_.mom.model_moments_params_derivs(jm,jp) = vec( pruned_state_space.dVar_yi(idx1,idx2,lag,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)');
+ repmat(pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)',[1 1 size(pruned_state_space.Var_yi,3)]);
end end
end end
end end
E_yyt_nbr = nnz(options_mom_.mom.index.E_yyt);
oo_.mom.model_moments(offset+(1:E_yyt_nbr),1) = E_yyt(options_mom_.mom.index.E_yyt);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(offset+(1:E_yyt_nbr),:) = reshape(dE_yyt(repmat(options_mom_.mom.index.E_yyt,[1 1 1 totparam_nbr])),E_yyt_nbr,totparam_nbr);
end end
end end
end
elseif strcmp(options_mom_.mom.mom_method,'SMM') elseif strcmp(options_mom_.mom.mom_method,'SMM')
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
...@@ -228,7 +221,7 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM') ...@@ -228,7 +221,7 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM')
y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order); y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order);
% provide meaningful penalty if data is nan or inf % provide meaningful penalty if data is nan or inf
if any(any(isnan(y_sim))) || any(any(isinf(y_sim))) if any(any(isnan(y_sim))) || any(any(isinf(y_sim)))
if options_mom_.mode_compute==13 if options_mom_.vector_output == 1 % lsqnonlin requires vector output
fval = Inf(size(oo_.mom.Sw,1),1); fval = Inf(size(oo_.mom.Sw,1),1);
else else
fval = Inf; fval = Inf;
...@@ -236,8 +229,8 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM') ...@@ -236,8 +229,8 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM')
info(1)=180; info(1)=180;
info(4) = 0.1; info(4) = 0.1;
exit_flag = 0; exit_flag = 0;
if options_mom_.mode_compute == 13 if options_mom_.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.dataMoments,1),1)*options_mom_.huge_number; fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
end end
return return
end end
...@@ -257,7 +250,7 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM') ...@@ -257,7 +250,7 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM')
if options_mom_.prefilter if options_mom_.prefilter
y_sim = bsxfun(@minus, y_sim, mean(y_sim,1)); y_sim = bsxfun(@minus, y_sim, mean(y_sim,1));
end end
oo_.mom.model_moments = method_of_moments_data_moments(y_sim, oo_, M_.matched_moments, options_mom_); oo_.mom.model_moments = mom.data_moments(y_sim, oo_, M_.matched_moments, options_mom_);
end end
...@@ -295,9 +288,9 @@ if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian ...@@ -295,9 +288,9 @@ if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
df(:,jp) = dresiduals; df(:,jp) = dresiduals;
end end
else else
df(:,jp) = dresiduals'*residuals + residuals'*dresiduals; df(jp,1) = dresiduals'*residuals + residuals'*dresiduals;
if options_mom_.mom.penalized_estimator if options_mom_.mom.penalized_estimator
df(:,jp)=df(:,jp)+(dxparam1(:,jp))'/oo_.prior.variance*(xparam1-oo_.prior.mean)+(xparam1-oo_.prior.mean)'/oo_.prior.variance*(dxparam1(:,jp)); df(jp,1)=df(jp,1)+(dxparam1(:,jp))'/oo_.prior.variance*(xparam1-oo_.prior.mean)+(xparam1-oo_.prior.mean)'/oo_.prior.variance*(dxparam1(:,jp));
end end
end end
end end
......
matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m → matlab/+mom/optimal_weighting_matrix.m
View file @ 2ddfa086
function W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag) function W_opt = optimal_weighting_matrix(m_data, moments, q_lag)
% W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag) % W_opt = optimal_weighting_matrix(m_data, moments, q_lag)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag q_lag % This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag q_lag
% Adapted from replication codes of % 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. % 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.
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% o m_data [T x numMom] selected data moments at each point in time % o m_data [T x numMom] selected data moments at each point in time
...@@ -14,7 +14,7 @@ function W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_l ...@@ -14,7 +14,7 @@ function W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_l
% o W_opt [numMom x numMom] optimal weighting matrix % o W_opt [numMom x numMom] optimal weighting matrix
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% o method_of_moments.m % o mom.run.m
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls: % This function calls:
% o CorrMatrix (embedded) % o CorrMatrix (embedded)
...@@ -67,6 +67,7 @@ end ...@@ -67,6 +67,7 @@ end
% The estimate of W % The estimate of W
W_opt = S\eye(size(S,1)); W_opt = S\eye(size(S,1));
W_opt=(W_opt+W_opt')/2; %assure symmetry
end end
% The correlation matrix % The correlation matrix
......
matlab/method_of_moments/method_of_moments.m → matlab/+mom/run.m
View file @ 2ddfa086
function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, estim_params_, M_, options_mom_) function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_, M_, options_mom_)
%function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, estim_params_, M_, options_mom_) %function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_, M_, options_mom_)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function performs a method of moments estimation with the following steps: % This function performs a method of moments estimation with the following steps:
% Step 0: Check if required structures and options exist % Step 0: Check if required structures and options exist
...@@ -49,11 +49,11 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, ...@@ -49,11 +49,11 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
% o get_all_parameters.m % o get_all_parameters.m
% o get_matrix_entries_for_psd_check.m % o get_matrix_entries_for_psd_check.m
% o makedataset.m % o makedataset.m
% o method_of_moments_check_plot.m % o mom.check_plot.m
% o method_of_moments_data_moments.m % o mom.data_moments.m
% o method_of_moments_objective_function.m % o mom.objective_function.m
% o method_of_moments_optimal_weighting_matrix % o mom.optimal_weighting_matrix
% o method_of_moments_standard_errors % o mom-standard_errors
% o plot_priors.m % o plot_priors.m
% o print_info.m % o print_info.m
% o prior_bounds.m % o prior_bounds.m
...@@ -63,7 +63,7 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, ...@@ -63,7 +63,7 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
% o set_all_parameters.m % o set_all_parameters.m
% o test_for_deep_parameters_calibration.m % o test_for_deep_parameters_calibration.m
% ========================================================================= % =========================================================================
% Copyright (C) 2020-2021 Dynare Team % Copyright (C) 2020-2022 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -91,8 +91,8 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, ...@@ -91,8 +91,8 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
% - [ ] add option to use autocorrelations (we have useautocorr in identification toolbox already) % - [ ] add option to use autocorrelations (we have useautocorr in identification toolbox already)
% - [ ] SMM with extended path % - [ ] SMM with extended path
% - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox) % - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
% - [ ] improve check for duplicate moments by using the cellfun and unique functions
% - [ ] dirname option to save output to different directory not yet implemented % - [ ] dirname option to save output to different directory not yet implemented
% - [ ] display scaled moments
% The TeX option crashes MATLAB R2014a run with "-nodisplay" option % The TeX option crashes MATLAB R2014a run with "-nodisplay" option
% (as is done from the testsuite). % (as is done from the testsuite).
...@@ -121,7 +121,7 @@ if isempty(estim_params_) % structure storing the info about estimated parameter ...@@ -121,7 +121,7 @@ if isempty(estim_params_) % structure storing the info about estimated parameter
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
end end
if isempty(M_.matched_moments) % structure storing the moments used for the method of moments estimation if ~isfield(M_,'matched_moments') || isempty(M_.matched_moments) % structure storing the moments used for the method of moments estimation
error('method_of_moments: You need to provide a ''matched_moments'' block') error('method_of_moments: You need to provide a ''matched_moments'' block')
end end
if ~isempty(bayestopt_) && any(bayestopt_.pshape==0) && any(bayestopt_.pshape~=0) if ~isempty(bayestopt_) && any(bayestopt_.pshape==0) && any(bayestopt_.pshape~=0)
...@@ -152,7 +152,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth ...@@ -152,7 +152,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
options_mom_.mom = set_default_option(options_mom_.mom,'bartlett_kernel_lag',20); % bandwith in optimal weighting matrix options_mom_.mom = set_default_option(options_mom_.mom,'bartlett_kernel_lag',20); % bandwith in optimal weighting matrix
options_mom_.mom = set_default_option(options_mom_.mom,'penalized_estimator',false); % include deviation from prior mean as additional moment restriction and use prior precision as weight options_mom_.mom = set_default_option(options_mom_.mom,'penalized_estimator',false); % include deviation from prior mean as additional moment restriction and use prior precision as weight
options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false); % display and store intermediate estimation results options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false); % display and store intermediate estimation results
options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix',{'DIAGONAL'; 'DIAGONAL'}); % weighting matrix in moments distance objective function at each iteration of estimation; options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix',{'DIAGONAL'; 'OPTIMAL'}); % weighting matrix in moments distance objective function at each iteration of estimation;
% possible values are 'OPTIMAL', 'IDENTITY_MATRIX' ,'DIAGONAL' or a filename. Size of cell determines stages in iterated estimation. % possible values are 'OPTIMAL', 'IDENTITY_MATRIX' ,'DIAGONAL' or a filename. Size of cell determines stages in iterated estimation.
options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1); % scaling of weighting matrix in objective function options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1); % scaling of weighting matrix in objective function
options_mom_.mom = set_default_option(options_mom_.mom,'se_tolx',1e-5); % step size for numerical computation of standard errors options_mom_.mom = set_default_option(options_mom_.mom,'se_tolx',1e-5); % step size for numerical computation of standard errors
...@@ -167,10 +167,10 @@ if strcmp(options_mom_.mom.mom_method,'SMM') ...@@ -167,10 +167,10 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
options_mom_.mom = set_default_option(options_mom_.mom,'burnin',500); % number of periods dropped at beginning of simulation options_mom_.mom = set_default_option(options_mom_.mom,'burnin',500); % number of periods dropped at beginning of simulation
options_mom_.mom = set_default_option(options_mom_.mom,'bounded_shock_support',false); % trim shocks in simulation to +- 2 stdev options_mom_.mom = set_default_option(options_mom_.mom,'bounded_shock_support',false); % trim shocks in simulation to +- 2 stdev
options_mom_.mom = set_default_option(options_mom_.mom,'seed',24051986); % seed used in simulations options_mom_.mom = set_default_option(options_mom_.mom,'seed',24051986); % seed used in simulations
options_mom_.mom = set_default_option(options_mom_.mom,'simulation_multiple',5); % multiple of the data length used for simulation options_mom_.mom = set_default_option(options_mom_.mom,'simulation_multiple',7); % multiple of the data length used for simulation
if options_mom_.mom.simulation_multiple < 1 if options_mom_.mom.simulation_multiple < 1
fprintf('The simulation horizon is shorter than the data. Dynare resets the simulation_multiple to 5.\n') fprintf('The simulation horizon is shorter than the data. Dynare resets the simulation_multiple to 5.\n')
options_mom_.mom.simulation_multiple = 5; options_mom_.mom.simulation_multiple = 7;
end end
end end
if strcmp(options_mom_.mom.mom_method,'GMM') if strcmp(options_mom_.mom.mom_method,'GMM')
...@@ -197,6 +197,7 @@ options_mom_ = set_default_option(options_mom_,'noprint',false); % do not ...@@ -197,6 +197,7 @@ options_mom_ = set_default_option(options_mom_,'noprint',false); % do not
options_mom_ = set_default_option(options_mom_,'plot_priors',true); % control plotting of priors options_mom_ = set_default_option(options_mom_,'plot_priors',true); % control plotting of priors
options_mom_ = set_default_option(options_mom_,'prior_trunc',1e-10); % probability of extreme values of the prior density that is ignored when computing bounds for the parameters options_mom_ = set_default_option(options_mom_,'prior_trunc',1e-10); % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
options_mom_ = set_default_option(options_mom_,'TeX',false); % print TeX tables and graphics options_mom_ = set_default_option(options_mom_,'TeX',false); % print TeX tables and graphics
options_mom_ = set_default_option(options_mom_,'verbosity',false); % print TeX tables and graphics
% Data and model options that can be set by the user in the mod file, otherwise default values are provided % Data and model options that can be set by the user in the mod file, otherwise default values are provided
options_mom_ = set_default_option(options_mom_,'first_obs',1); % number of first observation options_mom_ = set_default_option(options_mom_,'first_obs',1); % number of first observation
...@@ -205,6 +206,14 @@ options_mom_ = set_default_option(options_mom_,'nobs',NaN); % number of o ...@@ -205,6 +206,14 @@ options_mom_ = set_default_option(options_mom_,'nobs',NaN); % number of o
options_mom_ = set_default_option(options_mom_,'prefilter',false); % demean each data series by its empirical mean and use centered moments options_mom_ = set_default_option(options_mom_,'prefilter',false); % demean each data series by its empirical mean and use centered moments
options_mom_ = set_default_option(options_mom_,'xls_sheet',1); % name of sheet with data in Excel options_mom_ = set_default_option(options_mom_,'xls_sheet',1); % name of sheet with data in Excel
options_mom_ = set_default_option(options_mom_,'xls_range',''); % range of data in Excel sheet options_mom_ = set_default_option(options_mom_,'xls_range',''); % range of data in Excel sheet
% temporary workaround for https://git.dynare.org/Dynare/dseries/-/issues/51
if options_mom_.xls_sheet~=1
evalin('base','options_.xls_sheet=options_mom_.xls_sheet');
end
if ~isempty(options_mom_.xls_range)
evalin('base','options_.xls_range=options_mom_.xls_range');
end
% Recursive estimation and forecast are not supported % Recursive estimation and forecast are not supported
if numel(options_mom_.nobs)>1 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''.'); error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''nobs''.');
...@@ -215,7 +224,11 @@ end ...@@ -215,7 +224,11 @@ end
% Optimization options that can be set by the user in the mod file, otherwise default values are provided % Optimization options that can be set by the user in the mod file, otherwise default values are provided
options_mom_ = set_default_option(options_mom_,'huge_number',1e7); % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons options_mom_ = set_default_option(options_mom_,'huge_number',1e7); % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
options_mom_ = set_default_option(options_mom_,'mode_compute',13); % specifies the optimizer for minimization of moments distance if (isoctave && user_has_octave_forge_package('optim')) || (~isoctave && user_has_matlab_license('optimization_toolbox'))
options_mom_ = set_default_option(options_mom_,'mode_compute',13); % specifies lsqnonlin as default optimizer for minimization of moments distance
else
options_mom_ = set_default_option(options_mom_,'mode_compute',4); % specifies csminwel as fallback default option for minimization of moments distance
end
options_mom_ = set_default_option(options_mom_,'additional_optimizer_steps',[]); % vector of additional mode-finders run after mode_compute options_mom_ = set_default_option(options_mom_,'additional_optimizer_steps',[]); % vector of additional mode-finders run after mode_compute
options_mom_ = set_default_option(options_mom_,'optim_opt',[]); % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute options_mom_ = set_default_option(options_mom_,'optim_opt',[]); % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
options_mom_ = set_default_option(options_mom_,'silent_optimizer',false); % run minimization of moments distance silently without displaying results or saving files in between options_mom_ = set_default_option(options_mom_,'silent_optimizer',false); % run minimization of moments distance silently without displaying results or saving files in between
...@@ -359,7 +372,7 @@ analytic_jacobian_optimizers = [1, 3, 4, 13, 101]; %these are currently supporte ...@@ -359,7 +372,7 @@ analytic_jacobian_optimizers = [1, 3, 4, 13, 101]; %these are currently supporte
options_mom_.initialize_estimated_parameters_with_the_prior_mode = 0; % needed by set_prior.m options_mom_.initialize_estimated_parameters_with_the_prior_mode = 0; % needed by set_prior.m
options_mom_.figures.textwidth = 0.8; %needed by plot_priors.m options_mom_.figures.textwidth = 0.8; %needed by plot_priors.m
options_mom_.ramsey_policy = 0; % needed by evaluate_steady_state options_mom_.ramsey_policy = 0; % needed by evaluate_steady_state
options_mom_.debug = false; %needed by resol.m options_mom_ = set_default_option(options_mom_,'debug',false); %neeeded by e.g. check_plot
options_mom_.risky_steadystate = false; %needed by resol options_mom_.risky_steadystate = false; %needed by resol
options_mom_.threads = options_.threads; %needed by resol options_mom_.threads = options_.threads; %needed by resol
options_mom_.jacobian_flag = true; options_mom_.jacobian_flag = true;
...@@ -384,91 +397,78 @@ for i=1:options_mom_.obs_nbr ...@@ -384,91 +397,78 @@ for i=1:options_mom_.obs_nbr
end end
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% Step 2: Checks and transformations for matched moments structure (preliminary) % Step 2: Checks and transformations for matched_moments structure
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
M_.matched_moments_orig = M_.matched_moments; %save original structure
% Initialize indices
options_mom_.mom.index.E_y = false(options_mom_.obs_nbr,1); %unconditional first order product moments
options_mom_.mom.index.E_yy = false(options_mom_.obs_nbr,options_mom_.obs_nbr); %unconditional second order product moments
options_mom_.mom.index.E_yyt = false(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %unconditional temporal second order product moments
options_mom_.mom.index.E_y_pos = zeros(options_mom_.obs_nbr,1); %position in matched moments block
options_mom_.mom.index.E_yy_pos = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr); %position in matched moments block
options_mom_.mom.index.E_yyt_pos = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %position in matched moments block
for jm=1:size(M_.matched_moments,1)
% higher-order product moments not supported yet for GMM % higher-order product moments not supported yet for GMM
if strcmp(options_mom_.mom.mom_method, 'GMM') && sum(M_.matched_moments{jm,3}) > 2 if strcmp(options_mom_.mom.mom_method, 'GMM') && any(cellfun(@sum,M_.matched_moments(:,3))> 2)
error('method_of_moments: GMM does not yet support product moments higher than 2. Change row %d in ''matched_moments'' block.',jm); error('method_of_moments: GMM does not yet support product moments higher than 2. Change row %d in ''matched_moments'' block.',jm);
end end
% Check if declared variables are also observed (needed as otherwise the dataset variables won't coincide)
if any(~ismember(oo_.dr.inv_order_var(M_.matched_moments{jm,1})', oo_.dr.obs_var))
error('method_of_moments: Variables in row %d in ''matched_moments'' block need to be declared as VAROBS.', jm)
end
if strcmp(options_mom_.mom.mom_method, 'GMM') % check for duplicate moment conditions
% Check (for now) that only lags are declared for jm=1:size(M_.matched_moments,1)
if any(M_.matched_moments{jm,2}>0) % expand powers to vector of ones
error('method_of_moments: Leads in row %d in the ''matched_moments'' block are not supported for GMM, shift the moments and declare only lags.', jm) if any(M_.matched_moments{jm,3}>1)
end tmp1=[]; tmp2=[]; tmp3=[];
% Check (for now) that first declared variable has zero lag for jjm=1:length(M_.matched_moments{jm,3})
if M_.matched_moments{jm,2}(1)~=0 tmp1 = [tmp1 repmat(M_.matched_moments{jm,1}(jjm),[1 M_.matched_moments{jm,3}(jjm)]) ];
error('method_of_moments: The first variable declared in row %d in the ''matched_moments'' block is not allowed to have a lead or lag for GMM;\n reorder the variables in the row such that the first variable has zero lag!',jm) tmp2 = [tmp2 repmat(M_.matched_moments{jm,2}(jjm),[1 M_.matched_moments{jm,3}(jjm)]) ];
end tmp3 = [tmp3 repmat(1,[1 M_.matched_moments{jm,3}(jjm)]) ];
end end
vars = oo_.dr.inv_order_var(M_.matched_moments{jm,1})'; M_.matched_moments{jm,1} = tmp1;
if sum(M_.matched_moments{jm,3}) == 1 M_.matched_moments{jm,2} = tmp2;
% First-order product moment M_.matched_moments{jm,3} = tmp3;
vpos = (oo_.dr.obs_var == vars); end
options_mom_.mom.index.E_y(vpos,1) = true; % shift time structure to focus only on lags
options_mom_.mom.index.E_y_pos(vpos,1) = jm; M_.matched_moments{jm,2} = M_.matched_moments{jm,2} - max(M_.matched_moments{jm,2});
M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}},')']; % sort such that t=0 variable comes first
M_.matched_moments{jm,5}=['$E(',M_.endo_names_tex{M_.matched_moments{jm,1}},')$']; [M_.matched_moments{jm,2},idx_sort] = sort(M_.matched_moments{jm,2},'descend');
elseif sum(M_.matched_moments{jm,3}) == 2 M_.matched_moments{jm,1} = M_.matched_moments{jm,1}(idx_sort);
% Second-order product moment M_.matched_moments{jm,3} = M_.matched_moments{jm,3}(idx_sort);
idx1 = (oo_.dr.obs_var == vars(1)); end
idx2 = (oo_.dr.obs_var == vars(2));
lag1 = M_.matched_moments{jm,2}(1); % find duplicate rows in cell array by making groups according to powers as we can then use cell2mat for the unique function
lag2 = M_.matched_moments{jm,2}(2); powers = cellfun(@sum,M_.matched_moments(:,3))';
if lag1==0 && lag2==0 % contemporaneous covariance matrix UniqueMomIdx = [];
options_mom_.mom.index.E_yy(idx1,idx2) = true; for jpow = unique(powers)
options_mom_.mom.index.E_yy(idx2,idx1) = true; idx1 = find(powers==jpow);
options_mom_.mom.index.E_yy_pos(idx1,idx2) = jm; [~,idx2] = unique(cell2mat(M_.matched_moments(idx1,:)),'rows');
options_mom_.mom.index.E_yy_pos(idx2,idx1) = jm; UniqueMomIdx = [UniqueMomIdx idx1(idx2)];
M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}(1)},',',M_.endo_names{M_.matched_moments{jm,1}(2)},')']; end
M_.matched_moments{jm,5}=['$E({',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t,{',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t)$'];
elseif lag1==0 && lag2 < 0 % remove duplicate elements
options_mom_.mom.index.E_yyt(idx1,idx2,-lag2) = true; DuplicateMoms = setdiff(1:size(M_.matched_moments_orig,1),UniqueMomIdx);
options_mom_.mom.index.E_yyt_pos(idx1,idx2,-lag2) = jm;
M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}(1)},',',M_.endo_names{M_.matched_moments{jm,1}(2)},'(',num2str(lag2),'))'];
M_.matched_moments{jm,5}=['$E({',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t\times{',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'_{t',num2str(lag2) ,'})$'];
end
end
end
%Remove duplicate elements
UniqueMomIdx = [nonzeros(options_mom_.mom.index.E_y_pos); nonzeros(tril(options_mom_.mom.index.E_yy_pos)); nonzeros(options_mom_.mom.index.E_yyt_pos)];
DuplicateMoms = setdiff(1:size(M_.matched_moments,1),UniqueMomIdx);
if ~isempty(DuplicateMoms) if ~isempty(DuplicateMoms)
fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows: %s.\n',num2str(DuplicateMoms)) fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows:\n %s.\n',num2str(DuplicateMoms))
fprintf('Dynare will continue with remaining moment conditions\n');
end end
%reorder M_.matched_moments to be compatible with options_mom_.mom.index
M_.matched_moments = M_.matched_moments(UniqueMomIdx,:);
if strcmp(options_mom_.mom.mom_method, 'SMM') if strcmp(options_mom_.mom.mom_method, 'SMM')
options_mom_.mom=rmfield(options_mom_.mom,'index'); % for SMM we can keep the original structure but get rid of duplicate moments
M_.matched_moments = M_.matched_moments_orig(sort(UniqueMomIdx),:);
elseif strcmp(options_mom_.mom.mom_method, 'GMM')
% for GMM we use the transformed matched_moments structure
M_.matched_moments = M_.matched_moments(sort(UniqueMomIdx),:);
end end
% Check if both prefilter and first moments were specified % Check if both prefilter and first moments were specified
options_mom_.mom.first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,M_.matched_moments(:,3)))'; first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,M_.matched_moments(:,3)))';
if options_mom_.prefilter && ~isempty(options_mom_.mom.first_moment_indicator) if options_mom_.prefilter && ~isempty(first_moment_indicator)
fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block in rows: %u.\n',options_mom_.mom.first_moment_indicator'); fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block.\n');
M_.matched_moments(options_mom_.mom.first_moment_indicator,:)=[]; %remove first moments entries M_.matched_moments(first_moment_indicator,:)=[]; %remove first moments entries
options_mom_.mom.first_moment_indicator = [];
end end
options_mom_.mom.mom_nbr = size(M_.matched_moments,1); options_mom_.mom.mom_nbr = size(M_.matched_moments,1);
% Get maximum lag number for autocovariances/autocorrelations % Get maximum lag number for autocovariances/autocorrelations
options_mom_.ar = max(cellfun(@max,M_.matched_moments(:,2))) - min(cellfun(@min,M_.matched_moments(:,2))); 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_.dr.obs_var);
if ~isempty(not_observed_variables)
error('\nmethod_of_moments: You specified moments involving %s, but it is not a varobs.',M_.endo_names{oo_.dr.order_var(not_observed_variables)})
end
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% Step 3: Checks and transformations for estimated parameters, priors, and bounds % Step 3: Checks and transformations for estimated parameters, priors, and bounds
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
...@@ -497,7 +497,7 @@ end ...@@ -497,7 +497,7 @@ end
bayestopt_laplace=bayestopt_; bayestopt_laplace=bayestopt_;
% Check on specified priors and penalized estimation % Check on specified priors and penalized estimation
if any(bayestopt_laplace.pshape > 0) % prior specified, not ML if any(bayestopt_laplace.pshape > 0) % prior specified
if ~options_mom_.mom.penalized_estimator if ~options_mom_.mom.penalized_estimator
fprintf('\nPriors were specified, but the penalized_estimator-option was not set.\n') fprintf('\nPriors were specified, but the penalized_estimator-option was not set.\n')
fprintf('Dynare sets penalized_estimator to 1. Conducting %s with penalty.\n',options_mom_.mom.mom_method) fprintf('Dynare sets penalized_estimator to 1. Conducting %s with penalty.\n',options_mom_.mom.mom_method)
...@@ -530,8 +530,7 @@ if any(bayestopt_laplace.pshape > 0) % prior specified, not ML ...@@ -530,8 +530,7 @@ if any(bayestopt_laplace.pshape > 0) % prior specified, not ML
end end
end end
% Check for calibrated covariances before updating parameters estim_params_ = check_for_calibrated_covariances(estim_params_,M_);
estim_params_ = check_for_calibrated_covariances(xparam0,estim_params_,M_);
% Checks on parameter calibration and initialization % Checks on parameter calibration and initialization
xparam1_calib = get_all_parameters(estim_params_,M_); %get calibrated parameters xparam1_calib = get_all_parameters(estim_params_,M_); %get calibrated parameters
...@@ -643,7 +642,7 @@ end ...@@ -643,7 +642,7 @@ end
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'); 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 % Get data moments for the method of moments
[oo_.mom.data_moments, oo_.mom.m_data] = method_of_moments_data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_); [oo_.mom.data_moments, oo_.mom.m_data] = mom.data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_);
% Get shock series for SMM and set variance correction factor % Get shock series for SMM and set variance correction factor
if strcmp(options_mom_.mom.mom_method,'SMM') if strcmp(options_mom_.mom.mom_method,'SMM')
...@@ -722,7 +721,7 @@ test_for_deep_parameters_calibration(M_); ...@@ -722,7 +721,7 @@ test_for_deep_parameters_calibration(M_);
% If steady state of observed variables is non zero, set noconstant equal 0 % If steady state of observed variables is non zero, set noconstant equal 0
if all(abs(oo_.steady_state(oo_.dr.order_var(oo_.dr.obs_var)))<1e-9) if all(abs(oo_.steady_state(oo_.dr.order_var(oo_.dr.obs_var)))<1e-9)
options_mom_.noconstant = 1; options_mom_.noconstant = 0; %identifying the constant based on just the initial parameter value is not feasible
else else
options_mom_.noconstant = 0; options_mom_.noconstant = 0;
end end
...@@ -730,7 +729,7 @@ end ...@@ -730,7 +729,7 @@ end
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% Step 6: checks for objective function at initial parameters % Step 6: checks for objective function at initial parameters
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
objective_function = str2func('method_of_moments_objective_function'); objective_function = str2func('mom.objective_function');
try try
% Check for NaN or complex values of moment-distance-funtion evaluated % Check for NaN or complex values of moment-distance-funtion evaluated
...@@ -855,7 +854,6 @@ if size(options_mom_.mom.weighting_matrix,1)>1 && ~(any(strcmpi('diagonal',optio ...@@ -855,7 +854,6 @@ if size(options_mom_.mom.weighting_matrix,1)>1 && ~(any(strcmpi('diagonal',optio
fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n') fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n')
end end
optim_opt0 = options_mom_.optim_opt; % store original options set by user
for stage_iter=1:size(options_mom_.mom.weighting_matrix,1) for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
fprintf('Estimation stage %u\n',stage_iter); fprintf('Estimation stage %u\n',stage_iter);
Woptflag = false; Woptflag = false;
...@@ -867,19 +865,19 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1) ...@@ -867,19 +865,19 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
fprintf(' - diagonal of optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag); fprintf(' - diagonal of optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1 if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n'); fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = diag(diag( method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag) )); weighting_matrix = diag(diag( mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag) ));
else else
fprintf(' and using previous stage estimate of model-moments\n'); fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = diag(diag( method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag) )); weighting_matrix = diag(diag( mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag) ));
end end
case 'optimal' case 'optimal'
fprintf(' - optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag); fprintf(' - optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1 if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n'); fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag); weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag);
else else
fprintf(' and using previous stage estimate of model-moments\n'); fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag); weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
Woptflag = true; Woptflag = true;
end end
otherwise %user specified matrix in file otherwise %user specified matrix in file
...@@ -938,7 +936,7 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1) ...@@ -938,7 +936,7 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
[fval, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_); [fval, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization
SE = method_of_moments_standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag); SE = mom.standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag);
% Store results in output structure % Store results in output structure
oo_.mom = display_estimation_results_table(xparam1,SE,M_,options_mom_,estim_params_,bayestopt_laplace,oo_.mom,prior_dist_names,sprintf('%s (STAGE %u)',options_mom_.mom.mom_method,stage_iter),sprintf('%s_stage_%u',lower(options_mom_.mom.mom_method),stage_iter)); oo_.mom = display_estimation_results_table(xparam1,SE,M_,options_mom_,estim_params_,bayestopt_laplace,oo_.mom,prior_dist_names,sprintf('%s (STAGE %u)',options_mom_.mom.mom_method,stage_iter),sprintf('%s_stage_%u',lower(options_mom_.mom.mom_method),stage_iter));
...@@ -950,7 +948,7 @@ end ...@@ -950,7 +948,7 @@ end
if options_mom_.mom.mom_nbr > length(xparam1) if options_mom_.mom.mom_nbr > length(xparam1)
%get optimal weighting matrix for J test, if necessary %get optimal weighting matrix for J test, if necessary
if ~Woptflag if ~Woptflag
W_opt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag); W_opt = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
oo_j=oo_; oo_j=oo_;
oo_j.mom.Sw = chol(W_opt); oo_j.mom.Sw = chol(W_opt);
[fval] = feval(objective_function, xparam1, Bounds, oo_j, estim_params_, M_, options_mom_); [fval] = feval(objective_function, xparam1, Bounds, oo_j, estim_params_, M_, options_mom_);
...@@ -972,19 +970,43 @@ end ...@@ -972,19 +970,43 @@ end
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% Step 9: Display estimation results % Step 9: Display estimation results
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
title = ['Data moments and model moments (',options_mom_.mom.mom_method,')']; title = ['Comparison of data moments and model moments (',options_mom_.mom.mom_method,')'];
headers = {'Moment','Data','Model','% dev. target'}; headers = {'Moment','Data','Model'};
labels= M_.matched_moments(:,4); for jm = 1:size(M_.matched_moments,1)
data_mat=[oo_.mom.data_moments oo_.mom.model_moments 100*abs((oo_.mom.model_moments-oo_.mom.data_moments)./oo_.mom.data_moments)]; lables_tmp = 'E[';
lables_tmp_tex = 'E \left[ ';
for jvar = 1:length(M_.matched_moments{jm,1})
lables_tmp = [lables_tmp M_.endo_names{M_.matched_moments{jm,1}(jvar)}];
lables_tmp_tex = [lables_tmp_tex, '{', M_.endo_names_tex{M_.matched_moments{jm,1}(jvar)}, '}'];
if M_.matched_moments{jm,2}(jvar) ~= 0
lables_tmp = [lables_tmp, '(', num2str(M_.matched_moments{jm,2}(jvar)), ')'];
lables_tmp_tex = [lables_tmp_tex, '_{t', num2str(M_.matched_moments{jm,2}(jvar)), '}'];
else
lables_tmp_tex = [lables_tmp_tex, '_{t}'];
end
if M_.matched_moments{jm,3}(jvar) > 1
lables_tmp = [lables_tmp, '^', num2str(M_.matched_moments{jm,3}(jvar))];
lables_tmp_tex = [lables_tmp_tex, '^{', num2str(M_.matched_moments{jm,3}(jvar)) '}'];
end
if jvar == length(M_.matched_moments{jm,1})
lables_tmp = [lables_tmp, ']'];
lables_tmp_tex = [lables_tmp_tex, ' \right]'];
else
lables_tmp = [lables_tmp, '*'];
lables_tmp_tex = [lables_tmp_tex, ' \times '];
end
end
labels{jm,1} = lables_tmp;
labels_TeX{jm,1} = lables_tmp_tex;
end
data_mat=[oo_.mom.data_moments oo_.mom.model_moments ];
dyntable(options_mom_, title, headers, labels, data_mat, cellofchararraymaxlength(labels)+2, 10, 7); dyntable(options_mom_, title, headers, labels, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
if options_mom_.TeX if options_mom_.TeX
lh = cellofchararraymaxlength(labels)+2; dyn_latex_table(M_, options_mom_, title, ['comparison_moments_', options_mom_.mom.mom_method], headers, labels_TeX, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
labels_TeX = M_.matched_moments(:,5);
dyn_latex_table(M_, options_mom_, title, 'sim_corr_matrix', headers, labels_TeX, data_mat, lh, 10, 7);
end end
if options_mom_.mode_check.status if options_mom_.mode_check.status
method_of_moments_check_plot(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_laplace,... mom.check_plot(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_laplace,...
Bounds, oo_, estim_params_, M_, options_mom_) Bounds, oo_, estim_params_, M_, options_mom_)
end end
...@@ -994,8 +1016,11 @@ fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mo ...@@ -994,8 +1016,11 @@ fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mo
% Step 9: Clean up % Step 9: Clean up
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
%reset warning state %reset warning state
if isoctave warning_config;
warning('on')
else if isoctave && isfield(options_, 'prior_restrictions') && ...
warning on isfield(options_.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_.prior_restrictions.routine = [];
end end
matlab/method_of_moments/method_of_moments_standard_errors.m → matlab/+mom/standard_errors.m
View file @ 2ddfa086
function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Wopt_flag) function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Wopt_flag)
% [SE_values, Asympt_Var] = method_of_moments_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)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function computes standard errors to the method of moments estimates % This function computes standard errors to the method of moments estimates
% Adapted from replication codes of % Adapted from replication codes of
% o Andreasen, Fernndez-Villaverde, Rubio-Ramrez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49. % 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.
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% o xparam: value of estimated parameters as returned by set_prior() % o xparam: value of estimated parameters as returned by set_prior()
% o objective_function string of objective function, either method_of_moments_GMM.m or method_of_moments_SMM.m % o objective_function string of objective function
% o Bounds: structure containing parameter bounds % o Bounds: structure containing parameter bounds
% o oo_: structure for results % o oo_: structure for results
% o estim_params_: structure describing the estimated_parameters % o estim_params_: structure describing the estimated_parameters
...@@ -20,13 +20,13 @@ function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, obj ...@@ -20,13 +20,13 @@ function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, obj
% o Asympt_Var [nparam x nparam] asymptotic covariance matrix % o Asympt_Var [nparam x nparam] asymptotic covariance matrix
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% o method_of_moments.m % o mom.run.m
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls: % This function calls:
% o get_the_name % o get_the_name
% o get_error_message % o get_error_message
% o method_of_moments_objective_function % o mom.objective_function
% o method_of_moments_optimal_weighting_matrix % o mom.optimal_weighting_matrix
% ========================================================================= % =========================================================================
% Copyright (C) 2020-2021 Dynare Team % Copyright (C) 2020-2021 Dynare Team
% %
...@@ -87,10 +87,19 @@ else ...@@ -87,10 +87,19 @@ else
D(:,i) = (oo__p.mom.model_moments - oo__m.mom.model_moments)/(2*eps_value); D(:,i) = (oo__p.mom.model_moments - oo__m.mom.model_moments)/(2*eps_value);
else else
problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_); problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_);
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)
else
message_p = get_error_message(info_p, options_mom_); 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)
else
message_m = get_error_message(info_m, options_mom_); message_m = get_error_message(info_m, options_mom_);
end
warning('method_of_moments:info','Cannot compute the Jacobian for parameter %s - no standard errors available\n %s %s\nCheck your bounds and/or priors, or use a different optimizer.\n',problpar, message_p, message_m) 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)
end
Asympt_Var = NaN(length(xparam),length(xparam)); Asympt_Var = NaN(length(xparam),length(xparam));
SE_values = NaN(length(xparam),1); SE_values = NaN(length(xparam),1);
return return
...@@ -109,7 +118,7 @@ if Wopt_flag ...@@ -109,7 +118,7 @@ if Wopt_flag
Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params)); Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params));
else else
% We do not have the optimal weighting matrix yet % We do not have the optimal weighting matrix yet
WWopt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag); WWopt = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
S = WWopt\eye(size(WWopt,1)); S = WWopt\eye(size(WWopt,1));
AA = (D'*WW*D)\eye(dim_params); AA = (D'*WW*D)\eye(dim_params);
Asympt_Var = 1/T*AA*D'*WW*S*WW*D*AA; Asympt_Var = 1/T*AA*D'*WW*S*WW*D*AA;
......
matlab/+occbin/DSGE_smoother.m
View file @ 2ddfa086
...@@ -106,6 +106,8 @@ opts_simul.maxit = options_.occbin.smoother.maxit; ...@@ -106,6 +106,8 @@ opts_simul.maxit = options_.occbin.smoother.maxit;
opts_simul.waitbar = options_.occbin.smoother.waitbar; opts_simul.waitbar = options_.occbin.smoother.waitbar;
opts_simul.periods = options_.occbin.smoother.periods; opts_simul.periods = options_.occbin.smoother.periods;
opts_simul.check_ahead_periods = options_.occbin.smoother.check_ahead_periods; opts_simul.check_ahead_periods = options_.occbin.smoother.check_ahead_periods;
opts_simul.max_check_ahead_periods = options_.occbin.smoother.max_check_ahead_periods;
opts_simul.periodic_solution = options_.occbin.smoother.periodic_solution;
opts_simul.full_output = options_.occbin.smoother.full_output; opts_simul.full_output = options_.occbin.smoother.full_output;
opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only; 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 = options_.occbin.smoother.init_mode; % 0 = standard; 1 = unconditional frcsts zero shocks+smoothed states in each period
...@@ -138,6 +140,7 @@ opts_simul.exo_pos = 1:M_.exo_nbr; ...@@ -138,6 +140,7 @@ opts_simul.exo_pos = 1:M_.exo_nbr;
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1); opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues! opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
options_.noprint = true;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, out, ss] = occbin.solver(M_,oo_,options_);
regime_history = out.regime_history; regime_history = out.regime_history;
if options_.smoother_redux if options_.smoother_redux
...@@ -184,7 +187,9 @@ while is_changed && maxiter>iter && ~is_periodic ...@@ -184,7 +187,9 @@ while is_changed && maxiter>iter && ~is_periodic
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT; [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT;
sto_etahat(iter)={etahat}; sto_etahat(iter)={etahat};
regime_history0(iter,:) = regime_history; regime_history0(iter,:) = regime_history;
if occbin_smoother_debug
save('info1','regime_history0'); save('info1','regime_history0');
end
sto_CC = CC; sto_CC = CC;
sto_RR = RR; sto_RR = RR;
...@@ -250,7 +255,7 @@ while is_changed && maxiter>iter && ~is_periodic ...@@ -250,7 +255,7 @@ while is_changed && maxiter>iter && ~is_periodic
err_TT(iter-1) = max(max(max(abs(TT-sto_TT)))); err_TT(iter-1) = max(max(max(abs(TT-sto_TT))));
end end
if occbin_smoother_debug if occbin_smoother_debug || is_periodic
regime_ = cell(0); regime_ = cell(0);
regime_new = regime_; regime_new = regime_;
start_ = regime_; start_ = regime_;
...@@ -267,7 +272,12 @@ while is_changed && maxiter>iter && ~is_periodic ...@@ -267,7 +272,12 @@ while is_changed && maxiter>iter && ~is_periodic
qq={regime_history(indx_init_1).regimestart1}'; qq={regime_history(indx_init_1).regimestart1}';
for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end
disp('Time points where regime 1 differs') disp('Time points where regime 1 differs')
if ~isoctave
disp(table(indx_init_1, regime_, start_, regime_new, start_new)) disp(table(indx_init_1, regime_, start_, regime_new, start_new))
else % The table() function is not implemented in Octave, print something more or less equivalent (though much less readable)
disp(vertcat({'indx_init_1', 'regime_', 'start_', 'regime_new', 'start_new'}, ...
horzcat(num2cell(indx_init_1), regime_, start_, regime_new, start_new)))
end
end end
indx_init_2 = find(isdiff_(:,2)); indx_init_2 = find(isdiff_(:,2));
...@@ -285,7 +295,12 @@ while is_changed && maxiter>iter && ~is_periodic ...@@ -285,7 +295,12 @@ while is_changed && maxiter>iter && ~is_periodic
qq={regime_history(indx_init_2).regimestart2}'; qq={regime_history(indx_init_2).regimestart2}';
for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end
disp('Time points where regime 2 differs ') disp('Time points where regime 2 differs ')
if ~isoctave
disp(table(indx_init_2, regime_, start_, regime_new, start_new)) disp(table(indx_init_2, regime_, start_, regime_new, start_new))
else % The table() function is not implemented in Octave, print something more or less equivalent (though much less readable)
disp(vertcat({'indx_init_2', 'regime_', 'start_', 'regime_new', 'start_new'}, ...
horzcat(num2cell(indx_init_2), regime_, start_, regime_new, start_new)))
end
end end
else else
indx_init_1 = find(isdiff_(:,1)); indx_init_1 = find(isdiff_(:,1));
...@@ -299,7 +314,12 @@ while is_changed && maxiter>iter && ~is_periodic ...@@ -299,7 +314,12 @@ while is_changed && maxiter>iter && ~is_periodic
qq={regime_history(indx_init_1).regimestart}'; qq={regime_history(indx_init_1).regimestart}';
for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end
disp('Time points where regime differs') disp('Time points where regime differs')
if ~isoctave
disp(table(indx_init_1, regime_, start_, regime_new, start_new)) disp(table(indx_init_1, regime_, start_, regime_new, start_new))
else % The table() function is not implemented in Octave, print something more or less equivalent (though much less readable)
disp(vertcat({'indx_init_1', 'regime_', 'start_', 'regime_new', 'start_new'}, ...
horzcat(num2cell(indx_init_1), regime_, start_, regime_new, start_new)))
end
end end
end end
...@@ -312,17 +332,27 @@ end ...@@ -312,17 +332,27 @@ end
regime_history0(max(iter+1,1),:) = regime_history; regime_history0(max(iter+1,1),:) = regime_history;
oo_.occbin.smoother.regime_history=regime_history0(end,:); oo_.occbin.smoother.regime_history=regime_history0(end,:);
oo_.occbin.smoother.regime_history_iter=regime_history0; oo_.occbin.smoother.regime_history_iter=regime_history0;
if occbin_smoother_debug
save('info1','regime_history0') save('info1','regime_history0')
end
if (maxiter==iter && is_changed) || is_periodic if (maxiter==iter && is_changed) || is_periodic
disp(['Occbin smoother did not converge.']) 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')
if is_periodic if is_periodic
disp(['Occbin smoother algo loops between two solutions.']) 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')
else
fprintf('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.\n')
oo_.occbin.smoother=[];
oo_.occbin.smoother.error_flag=1;
end end
else else
disp(['Occbin smoother converged.']) disp('occbin.DSGE_smoother: smoother converged.')
oo_.occbin.smoother.error_flag=0;
if occbin_smoother_fast && is_changed_start if occbin_smoother_fast && is_changed_start
disp('WARNING: fast algo is used, regime(s) duration(s) was not forced to converge') disp('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge')
end end
end end
if (~is_changed || occbin_smoother_debug) && nargin==12 if (~is_changed || occbin_smoother_debug) && nargin==12
......
matlab/+occbin/IVF_core.m
View file @ 2ddfa086
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)
% function [filtered_errs, resids, Emat, stateval] = 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 % Computes thre
% %
% Outputs: % Outputs:
...@@ -19,7 +19,7 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o ...@@ -19,7 +19,7 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o
% - obs [T by N_obs] observed data % - obs [T by N_obs] observed data
% - init_val [N by 1] initial value of endogenous variables % - init_val [N by 1] initial value of endogenous variables
% Original authors: Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello and Molin Zhong % Original authors: Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello, and Molin Zhong
% Original file downloaded from: % Original file downloaded from:
% http://www.lguerrieri.com/jae-replication.zip % http://www.lguerrieri.com/jae-replication.zip
% Adapted for Dynare by Dynare Team. % Adapted for Dynare by Dynare Team.
...@@ -28,7 +28,7 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o ...@@ -28,7 +28,7 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o
% However the authors would appreciate acknowledgement of the source by % However the authors would appreciate acknowledgement of the source by
% citation of any of the following papers: % citation of any of the following papers:
% %
% Pablo Cuba-Borda, Luca Guerrieri, and Matteo Iacoviello (2019): "Likelihood evaluation of models % Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello, and Molin Zhong (2019): "Likelihood evaluation of models
% with occasionally binding constraints", Journal of Applied Econometrics, % with occasionally binding constraints", Journal of Applied Econometrics,
% 34(7), 1073-1085 % 34(7), 1073-1085
...@@ -60,6 +60,7 @@ opts_simul.maxit = options_.occbin.likelihood.maxit; ...@@ -60,6 +60,7 @@ opts_simul.maxit = options_.occbin.likelihood.maxit;
opts_simul.waitbar = false; opts_simul.waitbar = false;
opts_simul.periods = options_.occbin.likelihood.periods; opts_simul.periods = options_.occbin.likelihood.periods;
opts_simul.check_ahead_periods = options_.occbin.likelihood.check_ahead_periods; opts_simul.check_ahead_periods = options_.occbin.likelihood.check_ahead_periods;
opts_simul.max_check_ahead_periods = options_.occbin.likelihood.max_check_ahead_periods;
opts_simul.periodic_solution = options_.occbin.likelihood.periodic_solution; opts_simul.periodic_solution = options_.occbin.likelihood.periodic_solution;
opts_simul.restrict_state_space = options_.occbin.likelihood.restrict_state_space; opts_simul.restrict_state_space = options_.occbin.likelihood.restrict_state_space;
opts_simul.piecewise_only = 1; opts_simul.piecewise_only = 1;
......
matlab/+occbin/IVF_posterior.m
View file @ 2ddfa086
...@@ -50,7 +50,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti ...@@ -50,7 +50,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti
DLIK=[]; DLIK=[];
Hess=[]; Hess=[];
trend_coeff = []; trend_coeff = [];
obs = obs_info.rawdata; obs = dataset_.data;
obs_list = DynareOptions.varobs(:); obs_list = DynareOptions.varobs(:);
exit_flag = 1; exit_flag = 1;
...@@ -106,7 +106,13 @@ filtered_errs_init = zeros(sample_length,sum(err_index)); ...@@ -106,7 +106,13 @@ filtered_errs_init = zeros(sample_length,sum(err_index));
if info(1) if info(1)
fval = Inf; fval = Inf;
exit_flag = 0; exit_flag = 0;
atT=NaN(size(stateval(:,DynareResults.dr.order_var)'));
innov=NaN(Model.exo_nbr,sample_length);
return return
else
atT = stateval(:,DynareResults.dr.order_var)';
innov = zeros(Model.exo_nbr,sample_length);
innov(diag(Model.Sigma_e)~=0,:)=filtered_errs';
end end
nobs=size(filtered_errs,1); nobs=size(filtered_errs,1);
...@@ -185,14 +191,3 @@ end ...@@ -185,14 +191,3 @@ end
% remember that the likelihood has already been multiplied by -1 % remember that the likelihood has already been multiplied by -1
% hence, posterior is -1 times the log of the prior % hence, posterior is -1 times the log of the prior
fval = like+prior; fval = like+prior;
\ No newline at end of file
atT = stateval(:,DynareResults.dr.order_var)';
innov = zeros(Model.exo_nbr,sample_length);
innov(diag(Model.Sigma_e)~=0,:)=filtered_errs';
updated_variables = atT*nan;
BayesInfo.mf = BayesInfo.smoother_var_list(BayesInfo.smoother_mf);
initDynareOptions=DynareOptions;
DynareOptions.nk=[]; %unset options_.nk and reset it later
[DynareResults]=store_smoother_results(Model,DynareResults,DynareOptions,BayesInfo,dataset_,obs_info,atT,innov,[],updated_variables,DynareResults.dr.ys,zeros(length(DynareOptions.varobs_id),1));
DynareOptions=initDynareOptions;
matlab/+occbin/graph.m
View file @ 2ddfa086
...@@ -42,9 +42,9 @@ end ...@@ -42,9 +42,9 @@ end
var_list_plots=var_list(index_uniques); var_list_plots=var_list(index_uniques);
var_list_TeX = M_.endo_names_tex(i_var); var_list_TeX = M_.endo_names_tex(i_var);
data_to_plot(:,:,1)=oo_.occbin.piecewise(:,i_var); data_to_plot(:,:,1)=oo_.occbin.simul.piecewise(:,i_var);
if isfield(oo_.occbin,'linear') if isfield(oo_.occbin,'simul') && isfield(oo_.occbin.simul,'linear')
data_to_plot(:,:,2)=oo_.occbin.linear(:,i_var); data_to_plot(:,:,2)=oo_.occbin.simul.linear(:,i_var);
legend_list = {'Piecewise Linear','Linear'}; legend_list = {'Piecewise Linear','Linear'};
else else
legend_list = {'Piecewise Linear'}; legend_list = {'Piecewise Linear'};
...@@ -54,7 +54,7 @@ nperiods=size(data_to_plot,1); ...@@ -54,7 +54,7 @@ nperiods=size(data_to_plot,1);
ndim=size(data_to_plot,3); ndim=size(data_to_plot,3);
if ~options_.occbin.graph.steady_state if ~options_.occbin.graph.steady_state
data_to_plot=data_to_plot-repmat(oo_.occbin.ys(i_var)',nperiods,1,ndim); data_to_plot=data_to_plot-repmat(oo_.occbin.simul.ys(i_var)',nperiods,1,ndim);
end end
%get exogenous variables %get exogenous variables
...@@ -65,14 +65,14 @@ var_list_TeX = [var_list_TeX; M_.exo_names_tex(i_var_exo)]; ...@@ -65,14 +65,14 @@ var_list_TeX = [var_list_TeX; M_.exo_names_tex(i_var_exo)];
if number_of_plots_to_draw_exo>0 if number_of_plots_to_draw_exo>0
exo_index=NaN(number_of_plots_to_draw_exo); exo_index=NaN(number_of_plots_to_draw_exo);
for ii=1:length(i_var_exo) for ii=1:length(i_var_exo)
temp_index=find(oo_.occbin.exo_pos==i_var_exo(ii)); temp_index=find(oo_.occbin.simul.exo_pos==i_var_exo(ii));
if ~isempty(temp_index) if ~isempty(temp_index)
exo_index(ii)=temp_index; exo_index(ii)=temp_index;
else 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
end end
data_to_plot(:,end+1:end+number_of_plots_to_draw_exo,1)=[oo_.occbin.shocks_sequence(:,exo_index); zeros(nperiods-size(oo_.occbin.shocks_sequence,1),number_of_plots_to_draw_exo)]; 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)];
data_to_plot(:,end+1:end+number_of_plots_to_draw_exo,2)=NaN; data_to_plot(:,end+1:end+number_of_plots_to_draw_exo,2)=NaN;
end end
......
matlab/+occbin/kalman_update_algo_1.m
View file @ 2ddfa086
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_, info, 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)
% INPUTS % INPUTS
% - a [N by 1] t-1's state estimate % - a [N by 1] t-1's state estimate
% - a1 [N by 2] state predictions made at t-1:t % - a1 [N by 2] state predictions made at t-1:t
...@@ -42,7 +42,7 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal ...@@ -42,7 +42,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 % Philipp Pfeiffer, Marco Ratto (2021), Efficient and robust inference of models with occasionally binding
% constraints, Working Papers 2021-03, Joint Research Centre, European Commission % constraints, Working Papers 2021-03, Joint Research Centre, European Commission
% Copyright (C) 2021 Dynare Team % Copyright (C) 2021-2022 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -61,6 +61,12 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal ...@@ -61,6 +61,12 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal
warning off warning off
options_.noprint = true;
R=NaN(size(RR));
C=NaN(size(CC));
T=NaN(size(TT));
lik=Inf;
sto.a=a; sto.a=a;
sto.a1=a1; sto.a1=a1;
sto.P=P; sto.P=P;
...@@ -76,6 +82,7 @@ else ...@@ -76,6 +82,7 @@ else
base_regime.regime2 = 0; base_regime.regime2 = 0;
base_regime.regimestart2 = 1; base_regime.regimestart2 = 1;
end end
regimes_ = [base_regime base_regime base_regime];
mm=size(a,1); mm=size(a,1);
%% store info in t=1 %% store info in t=1
...@@ -94,17 +101,20 @@ PZI = P1(:,:,t)*ZZ'*iF(di,di,t); ...@@ -94,17 +101,20 @@ PZI = P1(:,:,t)*ZZ'*iF(di,di,t);
L(:,:,t) = eye(mm)-PZI*ZZ; L(:,:,t) = eye(mm)-PZI*ZZ;
if ~options_.occbin.filter.use_relaxation if ~options_.occbin.filter.use_relaxation
[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, 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);
else else
[~,~,~,~,~,~, TTx, RRx, CCx] ... [~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, base_regime,'reduced_state_space',T0,R0); = occbin.dynare_resolve(M_,options_,oo_, base_regime,'reduced_state_space',T0,R0);
regimes0(1)=base_regime;
TT(:,:,2) = TTx(:,:,end); TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end); RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,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); [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);
regimes0(1)=base_regime; end
if error_flag
etahat=NaN(size(QQQ,1),1);
return;
end end
%% run here the occbin simul %% run here the occbin simul
opts_simul = occbin_options.opts_simul; opts_simul = occbin_options.opts_simul;
...@@ -122,6 +132,11 @@ else ...@@ -122,6 +132,11 @@ else
end end
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, out, ss] = occbin.solver(M_,oo_,options_);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
return;
end
regimes_ = out.regime_history; regimes_ = out.regime_history;
if M_.occbin.constraint_nbr==1 if M_.occbin.constraint_nbr==1
...@@ -205,6 +220,11 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) ...@@ -205,6 +220,11 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.periods = max(opts_simul.periods,max(myregimestart)); opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, out, ss] = occbin.solver(M_,oo_,options_);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
return;
end
regimes0=regimes_; regimes0=regimes_;
regimes_ = out.regime_history; regimes_ = out.regime_history;
if niter>1 if niter>1
...@@ -250,6 +270,11 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) ...@@ -250,6 +270,11 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.maxit=1; opts_simul.maxit=1;
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, out, ss] = occbin.solver(M_,oo_,options_);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
return;
end
end end
end end
end end
...@@ -303,6 +328,11 @@ if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~ ...@@ -303,6 +328,11 @@ if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~
opts_simul.periods = max(opts_simul.periods,max(myregimestart)); opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, out, ss] = occbin.solver(M_,oo_,options_);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
return;
end
if isequal(out.regime_history(1),regimes_(1)) if isequal(out.regime_history(1),regimes_(1))
error_flag=0; error_flag=0;
break break
...@@ -312,26 +342,33 @@ if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~ ...@@ -312,26 +342,33 @@ if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~
end end
end end
if ~error_flag
a = out.piecewise(1:2,my_order_var)' - repmat(out.ys(my_order_var),1,2); a = out.piecewise(1:2,my_order_var)' - repmat(out.ys(my_order_var),1,2);
regimes_=regimes_(1:3);
end
T = ss.T(my_order_var,my_order_var,1:2); T = ss.T(my_order_var,my_order_var,1:2);
R = ss.R(my_order_var,:,1:2); R = ss.R(my_order_var,:,1:2);
C = ss.C(my_order_var,1:2); C = ss.C(my_order_var,1:2);
QQ = R(:,:,2)*QQQ(:,:,3)*transpose(R(:,:,2)); QQ = R(:,:,2)*QQQ(:,:,3)*transpose(R(:,:,2));
P(:,:,1) = P(:,:,2); P(:,:,1) = P(:,:,2);
P(:,:,2) = T(:,:,2)*P(:,:,1)*transpose(T(:,:,2))+QQ; P(:,:,2) = T(:,:,2)*P(:,:,1)*transpose(T(:,:,2))+QQ;
% P = cat(3,P(:,:,2),P2);
regimes_=regimes_(1:3);
etahat=etahat(:,2); etahat=etahat(:,2);
warning on warning_config;
end end
function [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, rescale_prediction_error_covariance, IF_likelihood) 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)
alphahat=NaN(size(a));
etahat=NaN(size(QQQ,1),2);
lik=Inf;
error_flag=0;
warning off warning off
if nargin<18 if nargin<18
IF_likelihood=0; IF_likelihood=0;
end end
t=2; t=2;
lik=0;
%% forward pass %% forward pass
% given updated variables and covarnace in t=1, we make the step to t=2 % given updated variables and covarnace in t=1, we make the step to t=2
T = TT(:,:,t); T = TT(:,:,t);
...@@ -354,6 +391,10 @@ else ...@@ -354,6 +391,10 @@ else
v(di,t) = Y(di,t) - ZZ*a(:,t); v(di,t) = Y(di,t) - ZZ*a(:,t);
F = ZZ*P(:,:,t)*ZZ' + H(di,di); F = ZZ*P(:,:,t)*ZZ' + H(di,di);
sig=sqrt(diag(F)); sig=sqrt(diag(F));
if any(any(isnan(F)))
error_flag=1;
return;
end
if rank(F)<size(F,1) if rank(F)<size(F,1)
% here we trap cases when some OBC regime triggers singularity % here we trap cases when some OBC regime triggers singularity
% e.g. no shock to interest rate at ZLB % e.g. no shock to interest rate at ZLB
...@@ -424,5 +465,5 @@ while t > 1 ...@@ -424,5 +465,5 @@ while t > 1
end end
end end
warning on warning_config;
end end
matlab/+occbin/kalman_update_algo_3.m
View file @ 2ddfa086
...@@ -69,12 +69,21 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, ...@@ -69,12 +69,21 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat,
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
warning off
options_.noprint = true;
T=[];
R=[];
C=[];
regimes_ = struct();
if isempty(nk) if isempty(nk)
nk=1; nk=1;
end end
nk=max(nk,1); nk=max(nk,1);
opts_simul = occbin_options.opts_regime;base_regime = struct(); opts_simul = occbin_options.opts_regime;
base_regime = struct();
if M_.occbin.constraint_nbr==1 if M_.occbin.constraint_nbr==1
base_regime.regime = 0; base_regime.regime = 0;
base_regime.regimestart = 1; base_regime.regimestart = 1;
...@@ -121,6 +130,10 @@ end ...@@ -121,6 +130,10 @@ end
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, out, ss] = occbin.solver(M_,oo_,options_);
if out.error_flag
error_flag = out.error_flag;
return;
end
regimes_ = out.regime_history; regimes_ = out.regime_history;
if M_.occbin.constraint_nbr==1 if M_.occbin.constraint_nbr==1
...@@ -210,6 +223,10 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) ...@@ -210,6 +223,10 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.periods = max(opts_simul.periods,max(myregimestart)); opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, out, ss] = occbin.solver(M_,oo_,options_);
if out.error_flag
error_flag = out.error_flag;
return;
end
regimes0=regimes_; regimes0=regimes_;
regimes_ = out.regime_history; regimes_ = out.regime_history;
if niter>1 if niter>1
...@@ -290,7 +307,7 @@ if error_flag==0 && niter>options_.occbin.likelihood.max_number_of_iterations && ...@@ -290,7 +307,7 @@ if error_flag==0 && niter>options_.occbin.likelihood.max_number_of_iterations &&
CC(:,2) = CCx(:,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] = 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)'; opts_simul.SHOCKS(1,:) = etahat(:,2)';
if occbin_options.opts_algo.restrict_state_space if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1); tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1); tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1); opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
......
matlab/+occbin/match_function.m
View file @ 2ddfa086
...@@ -28,7 +28,7 @@ function [resids, grad, state_out, E, M_, out] = match_function(err_0, obs_list, ...@@ -28,7 +28,7 @@ function [resids, grad, state_out, E, M_, out] = match_function(err_0, obs_list,
% However the authors would appreciate acknowledgement of the source by % However the authors would appreciate acknowledgement of the source by
% citation of any of the following papers: % citation of any of the following papers:
% %
% Pablo Cuba-Borda, Luca Guerrieri, and Matteo Iacoviello (2019): "Likelihood evaluation of models % Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello, and Molin Zhong (2019): "Likelihood evaluation of models
% with occasionally binding constraints", Journal of Applied Econometrics, % with occasionally binding constraints", Journal of Applied Econometrics,
% 34(7), 1073-1085 % 34(7), 1073-1085
...@@ -37,19 +37,18 @@ options_.occbin.simul=opts_simul; ...@@ -37,19 +37,18 @@ options_.occbin.simul=opts_simul;
options_.occbin.simul.full_output=1; options_.occbin.simul.full_output=1;
options_.noprint = 1; options_.noprint = 1;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, out, ss] = occbin.solver(M_,oo_,options_);
state_out= out.piecewise(1,:)' - out.ys;
E = ss.R(:,opts_simul.exo_pos);
grad = ss.R(opts_simul.varobs_id,opts_simul.exo_pos);
nobs = size(obs_list,1); nobs = size(obs_list,1);
resids = zeros(nobs,1); resids = zeros(nobs,1);
if ~out.error_flag if ~out.error_flag
% -- add observation block in model ---% state_out= out.piecewise(1,:)' - out.ys;
% % put in model file
E = ss.R(:,opts_simul.exo_pos);
grad = ss.R(opts_simul.varobs_id,opts_simul.exo_pos);
resids = (out.piecewise(1,opts_simul.varobs_id)-current_obs)'; %-out.endo_ss.(obs_list{this_obs}); resids = (out.piecewise(1,opts_simul.varobs_id)-current_obs)'; %-out.endo_ss.(obs_list{this_obs});
else else
grad = NaN(length(opts_simul.varobs_id),length(opts_simul.exo_pos));
resids = resids+100; resids = resids+100;
end end
......
matlab/+occbin/mkdata.m
View file @ 2ddfa086
...@@ -30,7 +30,7 @@ function [data_mat]=mkdata(n_periods,dr_A,dr_B,endo_names,exo_names,wish_list,sh ...@@ -30,7 +30,7 @@ function [data_mat]=mkdata(n_periods,dr_A,dr_B,endo_names,exo_names,wish_list,sh
% given decision rule % given decision rule
neqs = size(endo_names,1); neqs = size(endo_names,1);
if nargin<9 if nargin<9 || isempty(var_init)
var_init = zeros(neqs,1); var_init = zeros(neqs,1);
end end
......
matlab/+occbin/mkdatap_anticipated_2constraints_dyn.m
View file @ 2ddfa086
...@@ -45,7 +45,7 @@ T = DM.decrulea; ...@@ -45,7 +45,7 @@ T = DM.decrulea;
CONST = zeros(n_vars,1); CONST = zeros(n_vars,1);
R = DM.decruleb; R = DM.decruleb;
if nargin<7 if nargin<7 || isempty(init)
init=zeros(n_vars,1); init=zeros(n_vars,1);
end end
...@@ -90,15 +90,55 @@ if T_max > 0 ...@@ -90,15 +90,55 @@ if T_max > 0
dictionary.ss(1).R = temp(:,n_vars+1:n_vars+n_exo); dictionary.ss(1).R = temp(:,n_vars+1:n_vars+n_exo);
dictionary.ss(1).C = temp(:,n_vars+n_exo+1:end); dictionary.ss(1).C = temp(:,n_vars+n_exo+1:end);
end end
end % we do not know what is the last binding regime between 10 01 and 11!\
ireg(T_max)=1; ireg(T_max)=1;
icount = 1;
else
icount=length(dictionary.ss); icount=length(dictionary.ss);
% check if last binding regime was already stored
tmp = 0*binding_indicator;
tmp(1:end-T_max+1,:) = binding_indicator(T_max:end,:);
if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
itmp = find(~any(bsxfun(@minus, dictionary.binding_indicator(1:length(tmp)*2,:), tmp(:))));
else
itmp = find(~any(dictionary.binding_indicator(1:length(tmp)*2,:)-tmp(:)));
end
if ~isempty(itmp)
ireg(T_max) = itmp;
else
icount=icount+1;
ireg(T_max) = icount;
tmp = [binding_indicator(T_max,:); zeros(size(binding_indicator,1)-1,2)];
dictionary.binding_indicator(:,icount) = tmp(:);
if (binding_indicator(T_max,1) && ~binding_indicator(T_max,2))
temp = -(DM.Abarmat10*DM.decrulea+DM.Bbarmat10)\[DM.Cbarmat10 DM.Jbarmat10 DM.Dbarmat10];
dictionary.ss(icount).T = temp(:,1:n_vars);
dictionary.ss(icount).R = temp(:,n_vars+1:n_vars+n_exo);
dictionary.ss(icount).C = temp(:,n_vars+n_exo+1:end);
elseif (binding_indicator(T_max,1) && binding_indicator(T_max,2))
temp = -(DM.Abarmat11*DM.decrulea+DM.Bbarmat11)\[DM.Cbarmat11 DM.Jbarmat11 DM.Dbarmat11];
dictionary.ss(icount).T = temp(:,1:n_vars);
dictionary.ss(icount).R = temp(:,n_vars+1:n_vars+n_exo);
dictionary.ss(icount).C = temp(:,n_vars+n_exo+1:end);
else
temp = -(DM.Abarmat01*DM.decrulea+DM.Bbarmat01)\[DM.Cbarmat01 DM.Jbarmat01 DM.Dbarmat01];
dictionary.ss(icount).T = temp(:,1:n_vars);
dictionary.ss(icount).R = temp(:,n_vars+1:n_vars+n_exo);
dictionary.ss(icount).C = temp(:,n_vars+n_exo+1:end);
end
end
end
for i = T_max-1:-1:1 for i = T_max-1:-1:1
tmp = 0*binding_indicator; tmp = 0*binding_indicator;
tmp(1:end-i+1,:) = binding_indicator(i:end,:); tmp(1:end-i+1,:) = binding_indicator(i:end,:);
if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
itmp = find(~any(bsxfun(@minus, dictionary.binding_indicator(1:length(tmp)*2,:), tmp(:))));
else
itmp = find(~any(dictionary.binding_indicator(1:length(tmp)*2,:)-tmp(:))); itmp = find(~any(dictionary.binding_indicator(1:length(tmp)*2,:)-tmp(:)));
end
if ~isempty(itmp) if ~isempty(itmp)
ireg(i) = itmp; ireg(i) = itmp;
else else
......
matlab/+occbin/mkdatap_anticipated_dyn.m
View file @ 2ddfa086
...@@ -48,7 +48,7 @@ T = DM.decrulea; ...@@ -48,7 +48,7 @@ T = DM.decrulea;
CONST = zeros(n_vars,1); CONST = zeros(n_vars,1);
R = DM.decruleb; R = DM.decruleb;
if nargin<7 if nargin<7 || isempty(init)
init=zeros(n_vars,1); init=zeros(n_vars,1);
end end
...@@ -89,7 +89,11 @@ if T_max > 0 ...@@ -89,7 +89,11 @@ if T_max > 0
tmp = 0*binding_indicator; tmp = 0*binding_indicator;
tmp(1:end-i+1) = binding_indicator(i:end); tmp(1:end-i+1) = binding_indicator(i:end);
if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
itmp = find(~any(bsxfun(@minus, dictionary.binding_indicator, tmp)));
else
itmp = find(~any(dictionary.binding_indicator-tmp)); itmp = find(~any(dictionary.binding_indicator-tmp));
end
if ~isempty(itmp) if ~isempty(itmp)
ireg(i) = itmp; ireg(i) = itmp;
else else
......
matlab/+occbin/set_default_options.m
View file @ 2ddfa086
...@@ -37,6 +37,7 @@ if ismember(flag,{'all'}) ...@@ -37,6 +37,7 @@ if ismember(flag,{'all'})
options_occbin_.solver.solve_tolf=1e-5; options_occbin_.solver.solve_tolf=1e-5;
options_occbin_.solver.maxit=10; options_occbin_.solver.maxit=10;
options_occbin_.write_regimes.periods=[]; options_occbin_.write_regimes.periods=[];
options_occbin_.write_regimes.type='simul';
options_occbin_.write_regimes.filename=[M_.fname '_occbin_regimes']; options_occbin_.write_regimes.filename=[M_.fname '_occbin_regimes'];
end end
...@@ -66,21 +67,22 @@ end ...@@ -66,21 +67,22 @@ end
if ismember(flag,{'likelihood','all'}) if ismember(flag,{'likelihood','all'})
options_occbin_.likelihood.curb_retrench = false; options_occbin_.likelihood.curb_retrench = false;
options_occbin_.likelihood.first_period_occbin_update = true; options_occbin_.likelihood.first_period_occbin_update = 1;
options_occbin_.likelihood.full_output = false; options_occbin_.likelihood.full_output = false;
options_occbin_.likelihood.IF_likelihood = false; options_occbin_.likelihood.IF_likelihood = false;
options_occbin_.likelihood.init_regime_history = []; options_occbin_.likelihood.init_regime_history = [];
options_occbin_.likelihood.init_binding_indicator = false(0); options_occbin_.likelihood.init_binding_indicator = false(0);
options_occbin_.likelihood.inversion_filter = false; options_occbin_.likelihood.inversion_filter = false;
options_occbin_.likelihood.IVF_shock_observable_mapping = []; options_occbin_.likelihood.IVF_shock_observable_mapping = [];
options_occbin_.likelihood.maxit = 50; % this is for occbin solver algo options_occbin_.likelihood.maxit = 30; % this is for occbin solver algo
options_occbin_.likelihood.max_number_of_iterations = 10; % this is for occbin_kalman_update loop options_occbin_.likelihood.max_number_of_iterations = 10; % this is for occbin_kalman_update loop
options_occbin_.likelihood.max_check_ahead_periods=inf;
options_occbin_.likelihood.periods = 100; options_occbin_.likelihood.periods = 100;
options_occbin_.likelihood.check_ahead_periods=200; options_occbin_.likelihood.check_ahead_periods=200;
options_occbin_.likelihood.periodic_solution=true; options_occbin_.likelihood.periodic_solution=false;
options_occbin_.likelihood.piecewise_only = true; options_occbin_.likelihood.piecewise_only = true;
options_occbin_.likelihood.restrict_state_space = true; options_occbin_.likelihood.restrict_state_space = true;
options_occbin_.likelihood.status=true; options_occbin_.likelihood.status=true; %initialized to false in default_option_values
options_occbin_.likelihood.use_updated_regime = true; options_occbin_.likelihood.use_updated_regime = true;
options_occbin_.likelihood.waitbar=false; options_occbin_.likelihood.waitbar=false;
end end
...@@ -168,13 +170,13 @@ if ismember(flag,{'simul','all'}) ...@@ -168,13 +170,13 @@ if ismember(flag,{'simul','all'})
options_occbin_.simul.init_regime=[]; options_occbin_.simul.init_regime=[];
options_occbin_.simul.init_binding_indicator=false(0); options_occbin_.simul.init_binding_indicator=false(0);
options_occbin_.simul.exo_pos=1:M_.exo_nbr; options_occbin_.simul.exo_pos=1:M_.exo_nbr;
options_occbin_.simul.local=true; options_occbin_.simul.maxit=30;
options_occbin_.simul.maxit=10; options_occbin_.simul.max_check_ahead_periods=inf;
options_occbin_.simul.max_periods=inf; options_occbin_.simul.periods=100;
options_occbin_.simul.periods=30;
options_occbin_.simul.check_ahead_periods=200; options_occbin_.simul.check_ahead_periods=200;
options_occbin_.simul.periodic_solution=true; options_occbin_.simul.periodic_solution=false;
options_occbin_.simul.piecewise_only = false; options_occbin_.simul.piecewise_only = false;
options_occbin_.simul.reset_check_ahead_periods_in_new_period = false;
options_occbin_.simul.reset_regime_in_new_period = false; options_occbin_.simul.reset_regime_in_new_period = false;
options_occbin_.simul.restrict_state_space=false; options_occbin_.simul.restrict_state_space=false;
options_occbin_.simul.SHOCKS=zeros(1,M_.exo_nbr); options_occbin_.simul.SHOCKS=zeros(1,M_.exo_nbr);
...@@ -185,7 +187,7 @@ if ismember(flag,{'smoother','all'}) ...@@ -185,7 +187,7 @@ if ismember(flag,{'smoother','all'})
options_occbin_.smoother.curb_retrench = false; options_occbin_.smoother.curb_retrench = false;
options_occbin_.smoother.debug = false; options_occbin_.smoother.debug = false;
options_occbin_.smoother.fast = false; options_occbin_.smoother.fast = false;
options_occbin_.smoother.first_period_occbin_update = true; options_occbin_.smoother.first_period_occbin_update = 1;
options_occbin_.smoother.full_output = false; options_occbin_.smoother.full_output = false;
% options.occbin.smoother.init_mode = 1; % 0 = standard; 1 = unconditional frcsts zero shocks+smoothed states in each period % options.occbin.smoother.init_mode = 1; % 0 = standard; 1 = unconditional frcsts zero shocks+smoothed states in each period
options_occbin_.smoother.init_regime_history = []; options_occbin_.smoother.init_regime_history = [];
...@@ -193,9 +195,11 @@ if ismember(flag,{'smoother','all'}) ...@@ -193,9 +195,11 @@ if ismember(flag,{'smoother','all'})
options_occbin_.smoother.inversion_filter = false; options_occbin_.smoother.inversion_filter = false;
options_occbin_.smoother.linear_smoother = true; options_occbin_.smoother.linear_smoother = true;
options_occbin_.smoother.maxit = 30; % this is for occbin solver algo options_occbin_.smoother.maxit = 30; % this is for occbin solver algo
options_occbin_.smoother.max_check_ahead_periods=inf;
options_occbin_.smoother.max_number_of_iterations = 10; % this is for smoother loop options_occbin_.smoother.max_number_of_iterations = 10; % this is for smoother loop
options_occbin_.smoother.periods = 100; options_occbin_.smoother.periods = 100;
options_occbin_.smoother.check_ahead_periods=200; options_occbin_.smoother.check_ahead_periods=200;
options_occbin_.smoother.periodic_solution=false;
options_occbin_.smoother.piecewise_only = true; options_occbin_.smoother.piecewise_only = true;
options_occbin_.smoother.plot = true; options_occbin_.smoother.plot = true;
options_occbin_.smoother.status=true; options_occbin_.smoother.status=true;
......
matlab/+occbin/setup.m
View file @ 2ddfa086
...@@ -32,13 +32,23 @@ options_ = occbin.set_option(options_,options_occbin_,'simul.periods'); ...@@ -32,13 +32,23 @@ options_ = occbin.set_option(options_,options_occbin_,'simul.periods');
options_ = occbin.set_option(options_,options_occbin_,'simul.curb_retrench'); options_ = occbin.set_option(options_,options_occbin_,'simul.curb_retrench');
options_ = occbin.set_option(options_,options_occbin_,'simul.maxit'); options_ = occbin.set_option(options_,options_occbin_,'simul.maxit');
options_ = occbin.set_option(options_,options_occbin_,'simul.check_ahead_periods'); options_ = occbin.set_option(options_,options_occbin_,'simul.check_ahead_periods');
options_ = occbin.set_option(options_,options_occbin_,'simul.debug');
options_ = occbin.set_option(options_,options_occbin_,'simul.periodic_solution');
options_ = occbin.set_option(options_,options_occbin_,'smoother.periods'); options_ = occbin.set_option(options_,options_occbin_,'smoother.periods');
options_ = occbin.set_option(options_,options_occbin_,'smoother.curb_retrench'); options_ = occbin.set_option(options_,options_occbin_,'smoother.curb_retrench');
options_ = occbin.set_option(options_,options_occbin_,'smoother.maxit'); options_ = occbin.set_option(options_,options_occbin_,'smoother.maxit');
options_ = occbin.set_option(options_,options_occbin_,'smoother.check_ahead_periods'); options_ = occbin.set_option(options_,options_occbin_,'smoother.check_ahead_periods');
options_ = occbin.set_option(options_,options_occbin_,'smoother.debug');
options_ = occbin.set_option(options_,options_occbin_,'smoother.periodic_solution');
options_ = occbin.set_option(options_,options_occbin_,'smoother.inversion_filter');
options_ = occbin.set_option(options_,options_occbin_,'likelihood.periods');
options_ = occbin.set_option(options_,options_occbin_,'likelihood.curb_retrench');
options_ = occbin.set_option(options_,options_occbin_,'likelihood.maxit');
options_ = occbin.set_option(options_,options_occbin_,'likelihood.check_ahead_periods');
options_ = occbin.set_option(options_,options_occbin_,'likelihood.periodic_solution');
options_ = occbin.set_option(options_,options_occbin_,'likelihood.max_number_of_iterations');
options_ = occbin.set_option(options_,options_occbin_,'filter.use_relaxation'); options_ = occbin.set_option(options_,options_occbin_,'filter.use_relaxation');
options_ = occbin.set_option(options_,options_occbin_,'likelihood.inversion_filter'); options_ = occbin.set_option(options_,options_occbin_,'likelihood.inversion_filter');
options_ = occbin.set_option(options_,options_occbin_,'smoother.inversion_filter');
if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks) if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks)
temp=zeros(max(cat(2,M_.surprise_shocks.periods)),M_.exo_nbr); temp=zeros(max(cat(2,M_.surprise_shocks.periods)),M_.exo_nbr);
...@@ -46,7 +56,7 @@ if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks) ...@@ -46,7 +56,7 @@ if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks)
ivar = M_.surprise_shocks(ii).exo_id; ivar = M_.surprise_shocks(ii).exo_id;
temp(M_.surprise_shocks(ii).periods,ivar) = M_.surprise_shocks(ii).value; temp(M_.surprise_shocks(ii).periods,ivar) = M_.surprise_shocks(ii).value;
end end
shock_index=~all(temp==0); shock_index=~all(temp==0,1);
options_.occbin.simul.SHOCKS=temp(:,shock_index); options_.occbin.simul.SHOCKS=temp(:,shock_index);
options_.occbin.simul.exo_pos=find(shock_index); options_.occbin.simul.exo_pos=find(shock_index);
end end