diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index b4fcff64563ec0664fe44308bae67d6f5f1ee701..062a623012a096c6eae59e8b4fd0f218a224082a 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -6387,10 +6387,6 @@ observed variables.
initially far from the steady state. This option is not compatible with
``analytical_derivation``.
- .. option:: lik_algo = INTEGER
-
- For internal use and testing only.
-
.. option:: conf_sig = DOUBLE
Level of significance of the confidence interval used for classical forecasting after
diff --git a/matlab/+mom/default_option_mom_values.m b/matlab/+mom/default_option_mom_values.m
index 9a27c76134200e06da184e80a547ea21b124bf43..38ad0cf201ad29036910e0d32564f524c6542311 100644
--- a/matlab/+mom/default_option_mom_values.m
+++ b/matlab/+mom/default_option_mom_values.m
@@ -205,7 +205,6 @@ options_mom_.solve_tolx = options_.solve_tolx;
options_mom_.steady = options_.steady;
options_mom_.steadystate = options_.steadystate;
options_mom_.steadystate_flag = options_.steadystate_flag;
-%options_mom_.steadystate_partial
options_mom_.threads = options_.threads; % needed by resol
options_mom_.debug = options_.debug; % debug option needed by some functions, e.g. check_plot
diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m
index 181e27f1ea2a73db38ca082cce1161413f6d6bfc..b9fe38ba666d78956165e48f78352f7d8ded05ce 100644
--- a/matlab/+mom/objective_function.m
+++ b/matlab/+mom/objective_function.m
@@ -101,7 +101,7 @@ end
%--------------------------------------------------------------------------
% Compute linear approximation around the deterministic steady state
-[dr, info, M_, oo_] = resol(0, M_, options_mom_, oo_);
+[oo_.dr, info, M_.params] = resol(0, M_, options_mom_, oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
% Return, with endogenous penalty when possible, if resol issues an error code
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
@@ -151,11 +151,11 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
stderrparam_nbr = estim_params_.nvx; % number of stderr parameters
corrparam_nbr = estim_params_.ncx; % number of corr parameters
totparam_nbr = stderrparam_nbr+corrparam_nbr+modparam_nbr;
- dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices
+ oo_.dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices
oo_.mom.model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr);
- pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.mom.obs_var, options_mom_.ar, 0, 1);
+ pruned_state_space = pruned_state_space_system(M_, options_mom_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 1);
else
- pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.mom.obs_var, options_mom_.ar, 0, 0);
+ pruned_state_space = pruned_state_space_system(M_, options_mom_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 0);
end
oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
for jm = 1:size(M_.matched_moments,1)
@@ -218,7 +218,7 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
scaled_shock_series = zeros(size(options_mom_.mom.shock_series)); % initialize
scaled_shock_series(:,i_exo_var) = options_mom_.mom.shock_series(:,i_exo_var)*chol_S; % set non-zero entries
% simulate series
- y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order);
+ y_sim = simult_(M_, options_mom_, oo_.dr.ys, oo_.dr, scaled_shock_series, options_mom_.order);
% provide meaningful penalty if data is nan or inf
if any(any(isnan(y_sim))) || any(any(isinf(y_sim)))
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m
index 79319b3868e28efd8fd8b9ab57e3d65fa2ec93ae..903df87394e698542c6350a7fecd43cfbb93c79b 100644
--- a/matlab/+occbin/DSGE_smoother.m
+++ b/matlab/+occbin/DSGE_smoother.m
@@ -1,4 +1,4 @@
-function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
+function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
%function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
% Runs a DSGE smoother with occasionally binding constraints
%
@@ -65,7 +65,7 @@ regime_history=[];
if options_.occbin.smoother.linear_smoother && nargin==12
%% linear smoother
options_.occbin.smoother.status=false;
- [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
+ [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
tmp_smoother=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
for jf=1:length(smoother_field_list)
oo_.occbin.linear_smoother.(smoother_field_list{jf}) = tmp_smoother.(smoother_field_list{jf});
@@ -120,7 +120,7 @@ occbin_options.first_period_occbin_update = options_.occbin.smoother.first_perio
occbin_options.opts_regime = opts_simul; % this builds the opts_simul options field needed by occbin.solver
occbin_options.opts_regime.binding_indicator = options_.occbin.likelihood.init_binding_indicator;
occbin_options.opts_regime.regime_history=options_.occbin.likelihood.init_regime_history;
-[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
+[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
oo_.occbin.smoother.realtime_regime_history = oo_.occbin.smoother.regime_history;
regime_history = oo_.occbin.smoother.regime_history;
@@ -143,7 +143,7 @@ opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
options_.occbin.simul=opts_simul;
options_.noprint = true;
-[~, out, ss] = occbin.solver(M_,oo_,options_);
+[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
print_info(out.error_flag, options_.noprint, options_)
@@ -153,7 +153,7 @@ end
regime_history = out.regime_history;
if options_.smoother_redux
occbin_options.opts_simul.restrict_state_space =1;
- [T0,R0] = dynare_resolve(M_,options_,oo_);
+ [T0,R0] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
oo_.occbin.linear_smoother.T0=T0;
oo_.occbin.linear_smoother.R0=R0;
% oo_.occbin.linear_smoother.T0=ss.T(oo_.dr.order_var,oo_.dr.order_var,1);
@@ -168,7 +168,7 @@ CC = ss.C(oo_.dr.order_var,:);
opts_regime.regime_history = regime_history;
opts_regime.binding_indicator = [];
-[TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, oo_, options_);
+[TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, options_ , oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
is_changed = ~isequal(regime_history0,regime_history);
if isempty(regime_history0)
regime_history0 = regime_history;
@@ -195,7 +195,7 @@ while is_changed && maxiter>iter && ~is_periodic
iter=iter+1;
fprintf('Occbin smoother iteration %u.\n', iter)
occbin_options.opts_regime.regime_history=regime_history;
- [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT;
+ [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT;
sto_etahat(iter)={etahat};
regime_history0(iter,:) = regime_history;
if occbin_smoother_debug
@@ -209,7 +209,7 @@ while is_changed && maxiter>iter && ~is_periodic
opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
options_.occbin.simul=opts_simul;
- [~, out, ss] = occbin.solver(M_,oo_,options_);
+ [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
print_info(out.error_flag, false, options_)
@@ -225,7 +225,7 @@ while is_changed && maxiter>iter && ~is_periodic
% CC = cat(2,CC(:,1),CC);
opts_regime.regime_history = regime_history;
- [TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, oo_, options_);
+ [TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, options_ , oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
is_changed = ~isequal(regime_history0(iter,:),regime_history);
if M_.occbin.constraint_nbr==2
for k=1:size(regime_history0,2)
diff --git a/matlab/+occbin/IVF_posterior.m b/matlab/+occbin/IVF_posterior.m
index db4d51ec06fd178c9e4620e5e079814266042795..726f082562674cfaa6ddb2a950baafdaf145f146 100644
--- a/matlab/+occbin/IVF_posterior.m
+++ b/matlab/+occbin/IVF_posterior.m
@@ -76,7 +76,7 @@ err_index=DynareOptions.occbin.likelihood.IVF_shock_observable_mapping; % err_in
COVMAT1 = Model.Sigma_e(err_index,err_index);
% Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
-[T,R,SteadyState,info,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+[T,R,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
diff --git a/matlab/+occbin/check_regimes.m b/matlab/+occbin/check_regimes.m
index a0dcf62c600f3aa96ad9462ee1a063f6d62dac48..c6cb51c159c8fd45d37dfecaf8060583d3da2a5e 100644
--- a/matlab/+occbin/check_regimes.m
+++ b/matlab/+occbin/check_regimes.m
@@ -1,5 +1,5 @@
-function [TT, RR, CC, regime_history] = check_regimes(TT, RR, CC, opts_regime, M_, oo_, options_)
-%-function function [TT, RR, CC, regime_history] = check_regimes(TT, RR, CC, opts_regime, M_, oo_, options_)
+function [TT, RR, CC, regime_history] = check_regimes(TT, RR, CC, opts_regime, M_, options_, dr ,steady_state, exo_steady_state, exo_det_steady_state)
+%function [TT, RR, CC, regime_history] = check_regimes(TT, RR, CC, opts_regime, M_, options_ dr ,steady_state, exo_steady_state, exo_det_steady_state)
%
% INPUTS
% - TT [N by N] transition matrix of state space
@@ -7,8 +7,11 @@ function [TT, RR, CC, regime_history] = check_regimes(TT, RR, CC, opts_regime, M
% - CC [N by 1] constant of state space
% - opts_regime_ [structure] structure describing the regime
% - M_ [structure] Matlab's structure describing the model
-% - oo_ [structure] Matlab's structure containing the results
% - options_ [structure] Matlab's structure describing the current options
+% - dr [structure] Reduced form model.
+% - endo_steady_state [vector] steady state value for endogenous variables
+% - exo_steady_state [vector] steady state value for exogenous variables
+% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - TT [N by N] transition matrix of state space for each period
@@ -135,10 +138,10 @@ for tp=1:gend-1
opts_simul.maxit=1;
opts_simul.periods=nperi;
options_.occbin.simul=opts_simul;
- [~, ~, ss] = occbin.solver(M_,oo_,options_);
- TT(:,:,tp+1) = ss.T(oo_.dr.order_var,oo_.dr.order_var);
- RR(:,:,tp+1) = ss.R(oo_.dr.order_var,:);
- CC(:,tp+1) = ss.C(oo_.dr.order_var);
+ [~, ~, ss] = occbin.solver(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
+ TT(:,:,tp+1) = ss.T(dr.order_var,dr.order_var);
+ RR(:,:,tp+1) = ss.R(dr.order_var,:);
+ CC(:,tp+1) = ss.C(dr.order_var);
end
end
diff --git a/matlab/+occbin/dynare_resolve.m b/matlab/+occbin/dynare_resolve.m
index 6402853c0d29a6d16584acb67d0f4f1ebf8e2e8a..6c704854220ee64c21196f09171f7a121d99e6e0 100644
--- a/matlab/+occbin/dynare_resolve.m
+++ b/matlab/+occbin/dynare_resolve.m
@@ -1,4 +1,4 @@
-function [A,B,ys,info,M_,oo_,TT, RR, CC, A0, B0] ...
+function [A,B,ys,info,dr,params,TT, RR, CC, A0, B0] ...
= dynare_resolve(M_,options_,oo_,regime_history, reduced_state_space, A, B)
% function [A,B,ys,info,M_,options_,oo_,TT, RR, CC, A0, B0] ...
% = dynare_resolve(M_,options_,oo_,regime_history, reduced_state_space, A, B)
@@ -12,14 +12,14 @@ function [A,B,ys,info,M_,oo_,TT, RR, CC, A0, B0] ...
% - reduced_state_space [string]
% - A [double] State transition matrix
% - B [double] shock impact matrix
+%
% Outputs:
% - A [double] State transition matrix (potentially for restricted state space)
% - B [double] shock impact matrix (potentially for restricted state space)
% - ys [double] vector of steady state values
% - info [double] 4 by 1 vector with exit flag and information
-% - M_ [structure] Matlab's structure describing the model
-% - options_ [structure] Matlab's structure containing the options
-% - oo_ [structure] Matlab's structure containing the results
+% - dr [structure] Reduced form model.
+% - params [double] vector of potentially updated parameters
% - TT [N by N] transition matrix of state space for each period
% - RR [N by N_exo by T] shock impact matrix of state space for each period
% - CC [N by N_exo by T] constant of state space for each period
@@ -44,16 +44,18 @@ function [A,B,ys,info,M_,oo_,TT, RR, CC, A0, B0] ...
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if nargin<6
- [A,B,ys,info,M_,oo_] = dynare_resolve(M_,options_,oo_);
+ [A,B,ys,info,dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
ys = oo_.dr.ys;
+ dr = oo_.dr;
info = 0;
end
+params=M_.params;
if ~info(1) && nargin>4 && ~isempty(regime_history)
opts_regime.regime_history=regime_history;
opts_regime.binding_indicator=[];
[TT, RR, CC] = ...
- occbin.check_regimes(A, B, [], opts_regime, M_,oo_,options_);
+ occbin.check_regimes(A, B, [], opts_regime, M_, options_, dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
else
TT=A;
RR=B;
@@ -63,7 +65,7 @@ end
A0=A;
B0=B;
if ~info(1) && nargin>4 && ischar(reduced_state_space) && ~isempty(reduced_state_space)
- iv = oo_.dr.restrict_var_list;
+ iv = dr.restrict_var_list;
A=A(iv,iv);
B=B(iv,:);
TT=TT(iv,iv,:);
diff --git a/matlab/+occbin/kalman_update_algo_1.m b/matlab/+occbin/kalman_update_algo_1.m
index 8381e701110bbedcd76120ca777f495e19012f1b..2f304308665663fafb6afd81fe0cb7b7b5a237f3 100644
--- a/matlab/+occbin/kalman_update_algo_1.m
+++ b/matlab/+occbin/kalman_update_algo_1.m
@@ -131,10 +131,10 @@ else
my_order_var = oo_.dr.order_var;
end
options_.occbin.simul=opts_simul;
-[~, out, ss] = occbin.solver(M_,oo_,options_);
+[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
options_.occbin.simul.init_regime=regimes0;
- [~, out, ss] = occbin.solver(M_,oo_,options_);
+ [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
if out.error_flag
error_flag = out.error_flag;
@@ -223,7 +223,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
- [~, out, ss] = occbin.solver(M_,oo_,options_);
+ [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
@@ -273,7 +273,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
opts_simul.maxit=1;
options_.occbin.simul=opts_simul;
- [~, out, ss] = occbin.solver(M_,oo_,options_);
+ [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
@@ -331,7 +331,7 @@ if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
- [~, out, ss] = occbin.solver(M_,oo_,options_);
+ [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
diff --git a/matlab/+occbin/kalman_update_algo_3.m b/matlab/+occbin/kalman_update_algo_3.m
index 9894de15dc1ce46c55977060165d926311717271..c370920e239e756ac0ddbe980692cda9a1c5a1a7 100644
--- a/matlab/+occbin/kalman_update_algo_3.m
+++ b/matlab/+occbin/kalman_update_algo_3.m
@@ -129,7 +129,7 @@ else
end
options_.occbin.simul=opts_simul;
-[~, out, ss] = occbin.solver(M_,oo_,options_);
+[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
return;
@@ -222,7 +222,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
- [~, out, ss] = occbin.solver(M_,oo_,options_);
+ [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
return;
@@ -271,7 +271,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
opts_simul.maxit=1;
options_.occbin.simul=opts_simul;
- [~, out, ss] = occbin.solver(M_,oo_,options_);
+ [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
end
end
diff --git a/matlab/+occbin/match_function.m b/matlab/+occbin/match_function.m
index ae26236201fa04ca4d6f3ce21b6f129421cb5064..c1059470be7152ee1ac4e081dd26d6bceb3c89ba 100644
--- a/matlab/+occbin/match_function.m
+++ b/matlab/+occbin/match_function.m
@@ -36,7 +36,7 @@ opts_simul.SHOCKS = err_0';
options_.occbin.simul=opts_simul;
options_.occbin.simul.full_output=1;
options_.noprint = 1;
-[~, out, ss] = occbin.solver(M_,oo_,options_);
+[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
nobs = size(obs_list,1);
resids = zeros(nobs,1);
diff --git a/matlab/+occbin/solver.m b/matlab/+occbin/solver.m
index 3dc465ba3cfeb102ed811a6e78cd4fade92d275e..33b26e0289b8f81df16bafdf269427408e851262 100644
--- a/matlab/+occbin/solver.m
+++ b/matlab/+occbin/solver.m
@@ -1,5 +1,5 @@
-function [oo_, out, ss] = solver(M_,oo_,options_)
-% function [oo_, out, ss] = solver(M_,oo_,options_)
+function [dr, out, ss] = solver(M_, options_, dr ,steady_state, exo_steady_state, exo_det_steady_state)
+% function [oo_, out, ss] = solver(M_,oo_,options_, dr ,steady_state, exo_steady_state, exo_det_steady_state
% Solves the model with an OBC and produces simulations/IRFs
%
% INPUT:
@@ -58,17 +58,16 @@ if solve_dr
options_.order = 1;
end
- [dr,error_flag,M_,oo_] = resol(0,M_,options_,oo_);
+ [dr,error_flag,M_.params] = resol(0,M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
out.error_flag=error_flag;
if error_flag
print_info(error_flag, options_.noprint, options_)
return;
end
- oo_.dr = dr;
sto_dr=dr;
sto_M=M_;
else
- oo_.dr=sto_dr;
+ dr=sto_dr;
end
if options_.occbin.simul.check_ahead_periods>options_.occbin.simul.max_check_ahead_periods
@@ -78,9 +77,9 @@ if options_.occbin.simul.check_ahead_periods>options_.occbin.simul.max_check_ahe
end
if M_.occbin.constraint_nbr==1
- [out, ss, error_flag ] = occbin.solve_one_constraint(M_,oo_.dr,options_.occbin.simul,solve_dr);
+ [out, ss, error_flag ] = occbin.solve_one_constraint(M_,dr,options_.occbin.simul,solve_dr);
elseif M_.occbin.constraint_nbr==2
- [out, ss, error_flag ] = occbin.solve_two_constraints(M_,oo_.dr,options_.occbin.simul,solve_dr);
+ [out, ss, error_flag ] = occbin.solve_two_constraints(M_,dr,options_.occbin.simul,solve_dr);
end
out.error_flag=error_flag;
@@ -102,6 +101,4 @@ if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introdu
else
out.piecewise = out.piecewise + out.ys';
end
-out.exo_pos = options_.occbin.simul.exo_pos;
-
-oo_.occbin.simul=out;
+out.exo_pos = options_.occbin.simul.exo_pos;
\ No newline at end of file
diff --git a/matlab/+osr/objective.m b/matlab/+osr/objective.m
index 758d17738068bbdf46af65ce98036e4a265290a2..5905912411f3b5d5a700fdc42974fe46998de5c2 100644
--- a/matlab/+osr/objective.m
+++ b/matlab/+osr/objective.m
@@ -43,7 +43,7 @@ df=NaN(length(i_params),1);
% set parameters of the policy rule
M_.params(i_params) = x;
-[oo_.dr,info] = resol(0,M_,options_,oo_);
+[oo_.dr,info,M_.params] = resol(0,M_,options_,oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index f2a70113488e965b4059dd8ea4a837a33d5912cd..109ad465820990bfb38ebe4abd516b8011fdd374 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -1,4 +1,4 @@
-function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0,d] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin)
+function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0,d] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin)
% Estimation of the smoothed variables and innovations.
%
% INPUTS
@@ -32,7 +32,6 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
% o trend_addition [double] (n*T) pure trend component; stored in options_.varobs order
% o state_uncertainty [double] (K,K,T) array, storing the uncertainty
% about the smoothed state (decision-rule order)
-% o M_ [structure] decribing the model
% o oo_ [structure] storing the results
% o bayestopt_ [structure] describing the priors
%
@@ -105,7 +104,7 @@ if ~options_.smoother_redux
oo_.dr.restrict_var_list = bayestopt_.smoother_var_list;
oo_.dr.restrict_columns = bayestopt_.smoother_restrict_columns;
- [T,R,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_);
+ [T,R,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
%get location of observed variables and requested smoothed variables in
%decision rules
@@ -113,9 +112,9 @@ if ~options_.smoother_redux
else
if ~options_.occbin.smoother.status
- [T,R,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
+ [T,R,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
else
- [T,R,SteadyState,info,M_,oo_,~,~,~, T0, R0] = ...
+ [T,R,SteadyState,info,oo_.dr, M_.params,~,~,~, T0, R0] = ...
occbin.dynare_resolve(M_,options_,oo_,[],'restrict');
varargin{length_varargin+1}=T0;
varargin{length_varargin+2}=R0;
@@ -382,7 +381,7 @@ else
if size(TT,3)<(smpl+1)
TT=repmat(T,1,1,smpl+1);
RR=repmat(R,1,1,smpl+1);
- CC=repmat(zeros(mm,1),1,smpl+1);
+ CC=repmat(zeros(size(T,1),1),1,smpl+1);
end
end
if isoccbin==0
@@ -471,7 +470,7 @@ else
opts_simul.init_regime = []; %regimes_(k);
opts_simul.waitbar=0;
options_.occbin.simul=opts_simul;
- [~, out] = occbin.solver(M_,oo_,options_);
+ [~, out] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
% regime in out should be identical to regimes_(k-2) moved one
% period ahead (so if regimestart was [1 5] it should be [1 4]
% in out
@@ -521,7 +520,7 @@ else
opts_simul.init_regime = []; %regimes_(k);
opts_simul.waitbar=0;
options_.occbin.simul=opts_simul;
- [~, out] = occbin.solver(M_,oo_,options_);
+ [~, out] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
% regime in out should be identical to regimes_(k-2) moved one
% period ahead (so if regimestart was [1 5] it should be [1 4]
% in out
diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m
index 61f7361d973116943fd5eb7cdf196ede38fc71a3..c07fc92313968cc3fe6f81d6fb2c5b42b714ed23 100644
--- a/matlab/PosteriorIRF.m
+++ b/matlab/PosteriorIRF.m
@@ -1,10 +1,18 @@
-function PosteriorIRF(type,dispString)
+function oo_=PosteriorIRF(type,options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString)
+% PosteriorIRF(type,options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString)
% Builds posterior IRFs after the MH algorithm.
%
% INPUTS
-% o type [char] string specifying the joint density of the
-% deep parameters ('prior','posterior').
-% o dispString [char] string to display in the console.
+% o type [char] string specifying the joint density of the
+% deep parameters ('prior','posterior').
+% o options_ [structure] storing the options
+% o estim_params_ [structure] storing information about estimated parameters
+% o oo_ [structure] storing the results
+% o M_ [structure] storing the model information
+% o bayestopt_ [structure] storing information about priors
+% o dataset_ [structure] storing the dataset
+% o dataset_info [structure] Various information about the dataset
+% o dispString [char] string to display in the console.
%
% OUTPUTS
% None (oo_ and plots).
@@ -34,9 +42,6 @@ function PosteriorIRF(type,dispString)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-
-global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
-
% Set the number of periods
if isempty(options_.irf) || ~options_.irf
options_.irf = 40;
@@ -62,13 +67,7 @@ end
irf_shocks_indx = getIrfShocksIndx(M_, options_);
% Set various parameters & Check or create directories
-nvx = estim_params_.nvx;
-nvn = estim_params_.nvn;
-ncx = estim_params_.ncx;
-ncn = estim_params_.ncn;
-np = estim_params_.np ;
-npar = nvx+nvn+ncx+ncn+np;
-offset = npar-np; clear('nvx','nvn','ncx','ncn','np');
+npar = estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np ;
nvobs = dataset_.vobs;
gend = dataset_.nobs;
@@ -118,7 +117,8 @@ elseif strcmpi(type,'gsa')
end
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat istable
- B=size(x,1); options_.B = B;
+ B=size(x,1);
+ options_.B = B;
else% type = 'prior'
B = options_.prior_draws;
options_.B = B;
@@ -130,23 +130,7 @@ irun2 = 0;
NumberOfIRFfiles_dsge = 1;
NumberOfIRFfiles_dsgevar = 1;
ifil2 = 1;
-% Create arrays
-if B <= MAX_nruns
- stock_param = zeros(B, npar);
-else
- stock_param = zeros(MAX_nruns, npar);
-end
-if B >= MAX_nirfs_dsge
- stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
-else
- stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B);
-end
if MAX_nirfs_dsgevar
- if B >= MAX_nirfs_dsgevar
- stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
- else
- stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,B);
- end
NumberOfLags = options_.dsge_varlag;
NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
if options_.noconstant
@@ -154,10 +138,9 @@ if MAX_nirfs_dsgevar
else
NumberOfParametersPerEquation = NumberOfLagsTimesNvobs+1;
end
- Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
end
-% First block of code executed in parallel. The function devoted to do it is PosteriorIRF_core1.m
+% First block of code executed in parallel. The function devoted to do it is PosteriorIRF_core1.m
% function.
b = 0;
@@ -183,7 +166,6 @@ if ~strcmpi(type,'prior')
localVars.x=x;
end
-b=0;
if options_.dsge_var
localVars.gend = gend;
localVars.nvobs = nvobs;
@@ -202,6 +184,16 @@ localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
localVars.ifil2=ifil2;
localVars.MhDirectoryName=MhDirectoryName;
+%store main structures
+localVars.options_=options_;
+localVars.estim_params_= estim_params_;
+localVars.M_= M_;
+localVars.oo_= oo_;
+localVars.bayestopt_= bayestopt_;
+localVars.dataset_= dataset_;
+localVars.dataset_info= dataset_info;
+
+
% Like sequential execution!
if isnumeric(options_.parallel)
[fout] = PosteriorIRF_core1(localVars,1,B,0);
@@ -225,14 +217,6 @@ else
localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
localVars.ifil2=ifil2;
- globalVars = struct('M_',M_, ...
- 'options_', options_, ...
- 'bayestopt_', bayestopt_, ...
- 'estim_params_', estim_params_, ...
- 'oo_', oo_, ...
- 'dataset_',dataset_, ...
- 'dataset_info',dataset_info);
-
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '.static.m']};
NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']};
@@ -246,7 +230,7 @@ else
NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '.steadystate.m']};
end
end
- [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, globalVars, options_.parallel_info);
+ [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, [], options_.parallel_info);
nosaddle=0;
for j=1:length(fout)
nosaddle = nosaddle + fout(j).nosaddle;
@@ -260,9 +244,9 @@ if nosaddle
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
end
-ReshapeMatFiles('irf_dsge',type)
+ReshapeMatFiles(M_.fname,M_.dname,M_.exo_nbr,M_.endo_nbr,options_,'irf_dsge',type)
if MAX_nirfs_dsgevar
- ReshapeMatFiles('irf_bvardsge')
+ ReshapeMatFiles(M_.fname,M_.dname,M_.exo_nbr,M_.endo_nbr,options_,'irf_bvardsge')
end
if strcmpi(type,'gsa')
@@ -293,21 +277,21 @@ tit = M_.exo_names;
kdx = 0;
for file = 1:NumberOfIRFfiles_dsge
- load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']);
+ temp=load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']);
for i = irf_shocks_indx
for j = 1:nvar
- for k = 1:size(STOCK_IRF_DSGE,1)
+ for k = 1:size(temp.STOCK_IRF_DSGE,1)
kk = k+kdx;
[MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),...
- DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig);
+ DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(temp.STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig);
end
end
end
- kdx = kdx + size(STOCK_IRF_DSGE,1);
+ kdx = kdx + size(temp.STOCK_IRF_DSGE,1);
end
-clear STOCK_IRF_DSGE;
+clear temp;
for i = irf_shocks_indx
for j = 1:nvar
@@ -332,20 +316,20 @@ if MAX_nirfs_dsgevar
tit = M_.exo_names;
kdx = 0;
for file = 1:NumberOfIRFfiles_dsgevar
- load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']);
+ temp=load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']);
for i = irf_shocks_indx
for j = 1:nvar
- for k = 1:size(STOCK_IRF_BVARDSGE,1)
+ for k = 1:size(temp.STOCK_IRF_BVARDSGE,1)
kk = k+kdx;
[MeanIRFdsgevar(kk,j,i),MedianIRFdsgevar(kk,j,i),VarIRFdsgevar(kk,j,i),...
HPDIRFdsgevar(kk,:,j,i),DistribIRFdsgevar(kk,:,j,i)] = ...
- posterior_moments(squeeze(STOCK_IRF_BVARDSGE(k,j,i,:)),0,options_.mh_conf_sig);
+ posterior_moments(squeeze(temp.STOCK_IRF_BVARDSGE(k,j,i,:)),0,options_.mh_conf_sig);
end
end
end
- kdx = kdx + size(STOCK_IRF_BVARDSGE,1);
+ kdx = kdx + size(temp.STOCK_IRF_BVARDSGE,1);
end
- clear STOCK_IRF_BVARDSGE;
+ clear temp;
for i = irf_shocks_indx
for j = 1:nvar
name = sprintf('%s_%s', M_.endo_names{IndxVariables(j)}, tit{i});
@@ -358,25 +342,24 @@ if MAX_nirfs_dsgevar
end
end
end
-%
-% Finally I build the plots.
-%
+%% Finally I build the plots.
% Second block of code executed in parallel, with the exception of file
-% .tex generation always run in sequentially. This portion of code is execute in parallel by
+% .tex generation always run in sequentially. This portion of code is executed in parallel by
% PosteriorIRF_core2.m function.
if ~options_.nograph && ~options_.no_graph.posterior
% Save the local variables.
localVars=[];
-
+ localVars.dname=M_.dname;
+ localVars.fname=M_.fname;
+ localVars.options_=options_;
Check=options_.TeX;
if (Check)
localVars.varlist_TeX=varlist_TeX;
end
-
localVars.nvar=nvar;
localVars.MeanIRF=MeanIRF;
localVars.tit=tit;
@@ -445,9 +428,7 @@ if ~options_.nograph && ~options_.no_graph.posterior
if isRemoteOctave
[fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
else
- globalVars = struct('M_',M_, ...
- 'options_', options_);
-
+ globalVars = [];
[fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info);
end
end
@@ -455,7 +436,6 @@ if ~options_.nograph && ~options_.no_graph.posterior
[fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
end
% END parallel code!
-
end
-fprintf('%s: Posterior IRFs, done!\n',dispString);
+fprintf('%s: Posterior IRFs, done!\n',dispString);
\ No newline at end of file
diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m
index a0a47dd5ab32de3b92bbaf0e4bd66190cfb00a84..91d5b24ce13980015ed763bc766af6c20ee22e1c 100644
--- a/matlab/PosteriorIRF_core1.m
+++ b/matlab/PosteriorIRF_core1.m
@@ -1,4 +1,5 @@
function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
+%myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
% Generates and stores Posterior IRFs
% PARALLEL CONTEXT
% This function perfoms in parallel execution a portion of the PosteriorIRF.m code.
@@ -40,9 +41,6 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-
-global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
-
if nargin<4
whoiam=0;
end
@@ -50,6 +48,14 @@ end
% Reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
+options_= myinputs.options_;
+estim_params_= myinputs.estim_params_;
+M_= myinputs.M_;
+oo_= myinputs.oo_;
+bayestopt_= myinputs.bayestopt_;
+dataset_= myinputs.dataset_;
+dataset_info= myinputs.dataset_info;
+
IRUN = myinputs.IRUN;
irun =myinputs.irun;
irun2=myinputs.irun2;
@@ -78,13 +84,10 @@ if options_.dsge_var
bounds = prior_bounds(bayestopt_,options_.prior_trunc);
end
-
if whoiam
Parallel=myinputs.Parallel;
end
-
-% MhDirectoryName = myinputs.MhDirectoryName;
if strcmpi(type,'posterior')
MhDirectoryName = CheckPath('metropolis',M_.dname);
elseif strcmpi(type,'gsa')
@@ -115,12 +118,10 @@ else
h = dyn_waitbar(prct0,'Bayesian (prior) IRFs...');
end
-
OutputFileName_bvardsge = {};
OutputFileName_dsge = {};
OutputFileName_param = {};
-
fpar = fpar-1;
fpar0=fpar;
nosaddle=0;
@@ -150,8 +151,7 @@ while fpar<B
end
stock_param(irun2,:) = deep;
M_ = set_parameters_locally(M_, deep);
- [dr,info,M_,oo_] =compute_decision_rules(M_,options_,oo_);
- oo_.dr = dr;
+ [oo_.dr,info,M_.params] =compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if info(1)
nosaddle = nosaddle + 1;
fpar = fpar - 1;
@@ -182,9 +182,9 @@ while fpar<B
for i=irf_shocks_indx
if SS(i,i) > 5e-7
if options_.order>1 && options_.relative_irf % normalize shock to 0.01 before IRF generation for GIRFs; multiply with 100 later
- y=irf(M_,options_,dr,SS(:,i)./SS(i,i)/100, options_.irf, options_.drop,options_.replic,options_.order);
+ y=irf(M_,options_,oo_.dr,SS(:,i)./SS(i,i)/100, options_.irf, options_.drop,options_.replic,options_.order);
else
- y=irf(M_,options_,dr,SS(:,i), options_.irf, options_.drop,options_.replic,options_.order);
+ y=irf(M_,options_,oo_.dr,SS(:,i), options_.irf, options_.drop,options_.replic,options_.order);
end
if options_.relative_irf && options_.order==1 %multiply with 100 for backward compatibility
y = 100*y/SS(i,i);
diff --git a/matlab/PosteriorIRF_core2.m b/matlab/PosteriorIRF_core2.m
index 73ea4c493ecb330552ae884d2c5e82e19473b4f0..09b8650d9bc011bccb9b799a028cee73778ef1b2 100644
--- a/matlab/PosteriorIRF_core2.m
+++ b/matlab/PosteriorIRF_core2.m
@@ -1,5 +1,5 @@
function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam,ThisMatlab)
-% function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
+% myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% Generates the Posterior IRFs plot from the IRFs generated in
% PosteriorIRF_core1
%
@@ -47,8 +47,6 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam,ThisMatlab)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global options_ M_
-
if nargin<4
whoiam=0;
end
@@ -56,6 +54,7 @@ end
% Reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
+options_=myinputs.options_;
Check=options_.TeX;
if (Check)
varlist_TeX=myinputs.varlist_TeX;
@@ -67,6 +66,9 @@ tit=myinputs.tit;
nn=myinputs.nn;
MAX_nirfs_dsgevar=myinputs.MAX_nirfs_dsgevar;
HPDIRF=myinputs.HPDIRF;
+dname=myinputs.dname;
+fname=myinputs.fname;
+
if options_.dsge_var
HPDIRFdsgevar=myinputs.HPDIRFdsgevar;
MeanIRFdsgevar=myinputs.MeanIRFdsgevar;
@@ -82,7 +84,7 @@ end
% To save the figures where the function is computed!
-DirectoryName = CheckPath('Output',M_.dname);
+DirectoryName = CheckPath('Output',dname);
RemoteFlag = 0;
if whoiam
@@ -117,7 +119,6 @@ for i=fpar:npar
h2 = area(1:options_.irf,HPDIRF(:,1,j,i),'FaceColor',[1 1 1],'BaseValue',min(HPDIRF(:,1,j,i))); %white below HPDIinf and minimum of HPDIinf
end
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
- % plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
box on
axis tight
xlim([1 options_.irf]);
@@ -135,7 +136,6 @@ for i=fpar:npar
plot(1:options_.irf,HPDIRFdsgevar(:,1,j,i),'--k','linewidth',1)
plot(1:options_.irf,HPDIRFdsgevar(:,2,j,i),'--k','linewidth',1)
end
- % plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
box on
axis tight
xlim([1 options_.irf]);
@@ -155,9 +155,9 @@ for i=fpar:npar
if subplotnum == MaxNumberOfPlotPerFigure || (j == nvar && subplotnum> 0)
figunumber = figunumber+1;
- dyn_saveas(hh_fig,[DirectoryName '/' M_.fname '_Bayesian_IRF_' tit{i} '_' int2str(figunumber)],options_.nodisplay,options_.graph_format);
+ dyn_saveas(hh_fig,[DirectoryName '/' fname '_Bayesian_IRF_' tit{i} '_' int2str(figunumber)],options_.nodisplay,options_.graph_format);
if RemoteFlag==1
- OutputFileName = [OutputFileName; {[DirectoryName,filesep], [M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.*']}];
+ OutputFileName = [OutputFileName; {[DirectoryName,filesep], [fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.*']}];
end
subplotnum = 0;
end
@@ -165,7 +165,6 @@ for i=fpar:npar
if whoiam
fprintf('Done! \n');
waitbarString = [ 'Exog. shocks ' int2str(i) '/' int2str(npar) ' done.'];
- % fMessageStatus((i-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
dyn_waitbar((i-fpar+1)/(npar-fpar+1),[],waitbarString);
end
end % loop over exo_var
diff --git a/matlab/ReshapeMatFiles.m b/matlab/ReshapeMatFiles.m
index 38e529e59ae1ebc65271d967cf374d05a16f80ae..58107602eb15542bb828e7cbb244dc057b818580 100644
--- a/matlab/ReshapeMatFiles.m
+++ b/matlab/ReshapeMatFiles.m
@@ -1,10 +1,15 @@
-function ReshapeMatFiles(type, type2)
-% function ReshapeMatFiles(type, type2)
+function ReshapeMatFiles(fname, dname, exo_nbr, endo_nbr, options_, type, type2)
+% function ReshapeMatFiles(fname, dname, exo_nbr, endo_nbr, options_, type, type2)
% Reshapes and sorts (along the mcmc simulations) the mat files generated by DYNARE.
% 4D-arrays are splitted along the first dimension.
% 3D-arrays are splitted along the second dimension.
%
% INPUTS:
+% fname: [string] filename
+% dname: [string] directory name
+% exo_nbr: [integer] number of exogenous variables
+% endo_nbr: [integer] number of endogenous variables
+% options_: [struct] options structure
% type: statistics type in the repertory:
% dgse
% irf_bvardsge
@@ -25,7 +30,7 @@ function ReshapeMatFiles(type, type2)
% SPECIAL REQUIREMENTS
% none
-% Copyright © 2003-2017 Dynare Team
+% Copyright © 2003-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -42,43 +47,41 @@ function ReshapeMatFiles(type, type2)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global M_ options_
-
-if nargin==1
- MhDirectoryName = [ CheckPath('metropolis',M_.dname) filesep ];
+if nargin==6
+ MhDirectoryName = [ CheckPath('metropolis',dname) filesep ];
else
if strcmpi(type2,'posterior')
- MhDirectoryName = [CheckPath('metropolis',M_.dname) filesep ];
+ MhDirectoryName = [CheckPath('metropolis',dname) filesep ];
elseif strcmpi(type2,'gsa')
if options_.opt_gsa.morris==1
- MhDirectoryName = [CheckPath('gsa/screen',M_.dname) filesep ];
+ MhDirectoryName = [CheckPath('gsa/screen',dname) filesep ];
elseif options_.opt_gsa.morris==2
- MhDirectoryName = [CheckPath('gsa/identif',M_.dname) filesep ];
+ MhDirectoryName = [CheckPath('gsa/identif',dname) filesep ];
elseif options_.opt_gsa.pprior
- MhDirectoryName = [CheckPath(['gsa' filesep 'prior'],M_.dname) filesep ];
+ MhDirectoryName = [CheckPath(['gsa' filesep 'prior'],dname) filesep ];
else
- MhDirectoryName = [CheckPath(['gsa' filesep 'mc'],M_.dname) filesep ];
+ MhDirectoryName = [CheckPath(['gsa' filesep 'mc'],dname) filesep ];
end
else
- MhDirectoryName = [CheckPath('prior',M_.dname) filesep ];
+ MhDirectoryName = [CheckPath('prior',dname) filesep ];
end
end
switch type
case 'irf_dsge'
CAPtype = 'IRF_DSGE';
- TYPEsize = [ options_.irf , size(options_.varlist,1) , M_.exo_nbr ];
+ TYPEsize = [ options_.irf , size(options_.varlist,1) , exo_nbr ];
TYPEarray = 4;
case 'irf_bvardsge'
CAPtype = 'IRF_BVARDSGE';
- TYPEsize = [ options_.irf , length(options_.varobs) , M_.exo_nbr ];
+ TYPEsize = [ options_.irf , length(options_.varobs) , exo_nbr ];
TYPEarray = 4;
case 'smooth'
CAPtype = 'SMOOTH';
- TYPEsize = [ M_.endo_nbr , options_.nobs ];
+ TYPEsize = [ endo_nbr , options_.nobs ];
TYPEarray = 3;
case 'filter'
CAPtype = 'FILTER';
- TYPEsize = [ M_.endo_nbr , options_.nobs + 1 ];% TO BE CHECKED!
+ TYPEsize = [ endo_nbr , options_.nobs + 1 ];% TO BE CHECKED!
TYPEarray = 3;
case 'error'
CAPtype = 'ERROR';
@@ -86,22 +89,22 @@ switch type
TYPEarray = 3;
case 'innov'
CAPtype = 'INNOV';
- TYPEsize = [ M_.exo_nbr , options_.nobs ];
+ TYPEsize = [ exo_nbr , options_.nobs ];
TYPEarray = 3;
case 'forcst'
CAPtype = 'FORCST';
- TYPEsize = [ M_.endo_nbr , options_.forecast ];
+ TYPEsize = [ endo_nbr , options_.forecast ];
TYPEarray = 3;
case 'forcst1'
CAPtype = 'FORCST1';
- TYPEsize = [ M_.endo_nbr , options_.forecast ];
+ TYPEsize = [ endo_nbr , options_.forecast ];
TYPEarray = 3;
otherwise
disp('ReshapeMatFiles :: Unknown argument!')
return
end
-TYPEfiles = dir([MhDirectoryName M_.fname '_' type '*.mat']);
+TYPEfiles = dir([MhDirectoryName fname '_' type '*.mat']);
NumberOfTYPEfiles = length(TYPEfiles);
B = options_.B;
@@ -116,36 +119,29 @@ switch TYPEarray
for f1=1:NumberOfTYPEfiles-foffset
eval(['STOCK_' CAPtype ' = zeros(NumberOfPeriodsPerTYPEfiles,TYPEsize(2),TYPEsize(3),B);'])
for f2 = 1:NumberOfTYPEfiles
- load([MhDirectoryName M_.fname '_' type int2str(f2) '.mat']);
+ load([MhDirectoryName fname '_' type int2str(f2) '.mat']);
eval(['STOCK_' CAPtype '(:,:,1:+size(stock_' type ',3),idx+1:idx+size(stock_' type ',4))=stock_' ...
type '(jdx+1:jdx+NumberOfPeriodsPerTYPEfiles,:,:,:);'])
eval(['idx = idx + size(stock_' type ',4);'])
end
- %eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',4);'])
- save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(f1) '.mat'],['STOCK_' CAPtype]);
+ save([MhDirectoryName fname '_' CAPtype 's' int2str(f1) '.mat'],['STOCK_' CAPtype]);
jdx = jdx + NumberOfPeriodsPerTYPEfiles;
idx = 0;
end
if reste
eval(['STOCK_' CAPtype ' = zeros(reste,TYPEsize(2),TYPEsize(3),B);'])
for f2 = 1:NumberOfTYPEfiles
- load([MhDirectoryName M_.fname '_' type int2str(f2) '.mat']);
+ load([MhDirectoryName fname '_' type int2str(f2) '.mat']);
eval(['STOCK_' CAPtype '(:,:,:,idx+1:idx+size(stock_' type ',4))=stock_' type '(jdx+1:jdx+reste,:,:,:);'])
eval(['idx = idx + size(stock_' type ',4);'])
end
- %eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',4);'])
- save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(NumberOfTYPEfiles-foffset+1) '.mat'],['STOCK_' CAPtype]);
+ save([MhDirectoryName fname '_' CAPtype 's' int2str(NumberOfTYPEfiles-foffset+1) '.mat'],['STOCK_' CAPtype]);
end
else
- load([MhDirectoryName M_.fname '_' type '1.mat']);
- %eval(['STOCK_' CAPtype ' = sort(stock_' type ',4);'])
+ load([MhDirectoryName fname '_' type '1.mat']);
eval(['STOCK_' CAPtype ' = stock_' type ';'])
- save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]);
+ save([MhDirectoryName fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]);
end
- % Original file format may be useful in some cases...
- % for file = 1:NumberOfTYPEfiles
- % delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat'])
- % end
case 3
if NumberOfTYPEfiles>1
NumberOfPeriodsPerTYPEfiles = ceil( TYPEsize(2)/NumberOfTYPEfiles );
@@ -155,31 +151,24 @@ switch TYPEarray
for f1=1:NumberOfTYPEfiles-1
eval(['STOCK_' CAPtype ' = zeros(TYPEsize(1),NumberOfPeriodsPerTYPEfiles,B);'])
for f2 = 1:NumberOfTYPEfiles
- load([MhDirectoryName M_.fname '_' type int2str(f2) '.mat']);
+ load([MhDirectoryName fname '_' type int2str(f2) '.mat']);
eval(['STOCK_' CAPtype '(:,:,idx+1:idx+size(stock_ ' type ',3))=stock_' type '(:,jdx+1:jdx+NumberOfPeriodsPerTYPEfiles,:);'])
eval(['idx = idx + size(stock_' type ',3);'])
end
- %eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',3);'])
- save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(f1) '.mat'],['STOCK_' CAPtype]);
+ save([MhDirectoryName fname '_' CAPtype 's' int2str(f1) '.mat'],['STOCK_' CAPtype]);
jdx = jdx + NumberOfPeriodsPerTYPEfiles;
idx = 0;
end
eval(['STOCK_' CAPtype ' = zeros(TYPEsize(1),reste,B);'])
for f2 = 1:NumberOfTYPEfiles
- load([MhDirectoryName M_.fname '_' type int2str(f2) '.mat']);
+ load([MhDirectoryName fname '_' type int2str(f2) '.mat']);
eval(['STOCK_' CAPtype '(:,:,idx+1:idx+size(stock_' type ',3))=stock_' type '(:,jdx+1:jdx+reste,:);'])
eval(['idx = idx + size(stock_' type ',3);'])
end
- %eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',3);'])
- save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(NumberOfTYPEfiles) '.mat'],['STOCK_' CAPtype]);
+ save([MhDirectoryName fname '_' CAPtype 's' int2str(NumberOfTYPEfiles) '.mat'],['STOCK_' CAPtype]);
else
- load([MhDirectoryName M_.fname '_' type '1.mat']);
- %eval(['STOCK_' CAPtype ' = sort(stock_' type ',3);'])
+ load([MhDirectoryName fname '_' type '1.mat']);
eval(['STOCK_' CAPtype ' = stock_' type ';'])
- save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]);
+ save([MhDirectoryName fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]);
end
- % Original file format may be useful in some cases...
- % for file = 1:NumberOfTYPEfiles
- % delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat'])
- % end
end
\ No newline at end of file
diff --git a/matlab/backward/simul_backward_model_init.m b/matlab/backward/simul_backward_model_init.m
index 447a2eee458fa42fc21c30d183fbe10cc977cc64..504593fb618b6ca8f1b4aa81bbc18b99309bffd1 100644
--- a/matlab/backward/simul_backward_model_init.m
+++ b/matlab/backward/simul_backward_model_init.m
@@ -162,7 +162,7 @@ if nargin<6 || isempty(innovations)
case 'gaussian'
DynareOutput.bnlms.shocks = randn(samplesize,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
otherwise
- error(['simul_backward_nonlinear_model:: ' DynareOption.bnlms.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
+ error(['simul_backward_nonlinear_model:: ' DynareOptions.bnlms.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
% Put the simulated innovations in DynareOutput.exo_simul.
DynareOutput.exo_simul = zeros(samplesize,number_of_shocks);
diff --git a/matlab/check.m b/matlab/check.m
index ec9ec65136bc93e71309c656321b868b8df5eb70..57775f91833ba8bafef2602203ee17385e0656a7 100644
--- a/matlab/check.m
+++ b/matlab/check.m
@@ -42,7 +42,7 @@ end
oo.dr=set_state_space(oo.dr,M,options);
-[dr,info,M,~] = resol(1,M,options,oo);
+[dr,info] = resol(1,M,options,oo.dr ,oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state);
if info(1) ~= 0 && info(1) ~= 3 && info(1) ~= 4
print_info(info, 0, options);
diff --git a/matlab/cli/prior.m b/matlab/cli/prior.m
index a0956e56297a7238862e1fe668fb463606968c8d..dd0bb52b3ed0a49008825a85e48663ed865761a7 100644
--- a/matlab/cli/prior.m
+++ b/matlab/cli/prior.m
@@ -137,7 +137,7 @@ if ismember('moments', varargin) % Prior simulations (2nd order moments).
oo__ = oo_;
oo__.dr = set_state_space(oo__.dr, Model, options_);
% Solve model
- [T,R,~,info,Model , oo__] = dynare_resolve(Model , options_ ,oo__,'restrict');
+ [T,R,~,info,oo__.dr, Model.params] = dynare_resolve(Model , options_ , oo__.dr, oo__.steady_state, oo__.exo_steady_state, oo__.exo_det_steady_state,'restrict');
if ~info(1)
info=endogenous_prior_restrictions(T,R,Model , options__ , oo__);
end
diff --git a/matlab/compute_decision_rules.m b/matlab/compute_decision_rules.m
index 9c7b561637eb06a029f692b0c33ffbbceeba0293..b8dc6c470631530fecc7e3b10d7bd73b5f98e58f 100644
--- a/matlab/compute_decision_rules.m
+++ b/matlab/compute_decision_rules.m
@@ -1,17 +1,19 @@
-function [dr,info,M_,oo_] =compute_decision_rules(M_,options_,oo_)
-% function [dr,info,M_,oo_] =compute_decision_rules(M_,options_,oo_)
+function [dr,info,params] =compute_decision_rules(M_,options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
+% function [dr,info,params] =compute_decision_rules(M_,options_,oo_)
% INPUTS
% - M_ [structure] Matlab's structure describing the model (M_).
% - options_ [structure] Matlab's structure describing the current options (options_).
-% - oo_ [structure] Matlab's structure containing the results (oo_).
+% - dr [structure] Reduced form model.
+% - endo_steady_state [vector] steady state value for endogenous variables
+% - exo_steady_state [vector] steady state value for exogenous variables
+% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - dr [structure] Reduced form model.
% - info [integer] scalar or vector, error code.
-% - M_ [structure] Matlab's structure describing the model (M_).
-% - oo_ [structure] Matlab's structure containing the results (oo_).
+% - params [double] vector of potentially updated parameters
-% Copyright © 2020 Dynare Team
+% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -29,7 +31,7 @@ function [dr,info,M_,oo_] =compute_decision_rules(M_,options_,oo_)
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if options_.discretionary_policy
- [dr,info,M_,oo_] = discretionary_policy_1(options_.instruments,M_,options_,oo_);
+ [dr,info,params] = discretionary_policy_1(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
else
- [dr,info,M_,oo_] = resol(0,M_,options_,oo_);
+ [dr,info,params] = resol(0,M_,options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
end
diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m
index 8b9e8aae13eab750746be7a8d0162799fbd134df..e50f1adae8845c2e2dccdd31e15757ea5a6e9227 100644
--- a/matlab/compute_moments_varendo.m
+++ b/matlab/compute_moments_varendo.m
@@ -1,4 +1,4 @@
-function oo_ = compute_moments_varendo(type, options_, M_, oo_, var_list_)
+function oo_ = compute_moments_varendo(type, options_, M_, oo_, estim_params_, var_list_)
% Computes the second order moments (autocorrelation function, covariance
% matrix and variance decomposition) distributions for all the endogenous variables selected in
@@ -55,7 +55,7 @@ end
if strcmpi(type,'posterior')
posterior = 1;
- if nargin==4
+ if nargin==5
var_list_ = options_.varobs;
end
if isfield(oo_,'PosteriorTheoreticalMoments')
@@ -63,7 +63,7 @@ if strcmpi(type,'posterior')
end
elseif strcmpi(type,'prior')
posterior = 0;
- if nargin==4
+ if nargin==5
var_list_ = options_.prior_analysis_endo_var_list;
if isempty(var_list_)
options_.prior_analysis_var_list = options_.varobs;
@@ -97,13 +97,13 @@ end
if posterior
for i=1:NumberOfEndogenousVariables
for j=i:NumberOfEndogenousVariables
- oo_ = posterior_analysis('variance', var_list_{i}, var_list_{j}, NumberOfLags, options_, M_, oo_);
+ oo_ = posterior_analysis('variance', var_list_{i}, var_list_{j}, NumberOfLags, options_, M_, oo_, estim_params_);
end
end
else
for i=1:NumberOfEndogenousVariables
for j=i:NumberOfEndogenousVariables
- oo_ = prior_analysis('variance', var_list_{i}, var_list_{j}, [], options_, M_, oo_);
+ oo_ = prior_analysis('variance', var_list_{i}, var_list_{j}, [], options_, M_, oo_, estim_params_);
end
end
end
@@ -113,7 +113,7 @@ if posterior
for h=NumberOfLags:-1:1
for i=1:NumberOfEndogenousVariables
for j=1:NumberOfEndogenousVariables
- oo_ = posterior_analysis('correlation', var_list_{i}, var_list_{j}, h, options_, M_, oo_);
+ oo_ = posterior_analysis('correlation', var_list_{i}, var_list_{j}, h, options_, M_, oo_, estim_params_);
end
end
end
@@ -121,7 +121,7 @@ else
for h=NumberOfLags:-1:1
for i=1:NumberOfEndogenousVariables
for j=1:NumberOfEndogenousVariables
- oo_ = prior_analysis('correlation', var_list_{i}, var_list_{j}, h, options_, M_, oo_);
+ oo_ = prior_analysis('correlation', var_list_{i}, var_list_{j}, h, options_, M_, oo_, estim_params_);
end
end
end
@@ -135,7 +135,7 @@ if options_.order==1
if posterior
for i=1:NumberOfEndogenousVariables
for j=1:NumberOfExogenousVariables
- oo_ = posterior_analysis('decomposition', var_list_{i}, M_.exo_names{j}, [], options_, M_, oo_);
+ oo_ = posterior_analysis('decomposition', var_list_{i}, M_.exo_names{j}, [], options_, M_, oo_, estim_params_);
temp(i,j) = oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.Mean.(var_list_{i}).(M_.exo_names{j});
end
end
@@ -144,7 +144,7 @@ if options_.order==1
else
for i=1:NumberOfEndogenousVariables
for j=1:NumberOfExogenousVariables
- oo_ = prior_analysis('decomposition', var_list_{i}, M_.exo_names{j}, [], options_, M_, oo_);
+ oo_ = prior_analysis('decomposition', var_list_{i}, M_.exo_names{j}, [], options_, M_, oo_, estim_params_);
temp(i,j)=oo_.PriorTheoreticalMoments.dsge.VarianceDecomposition.Mean.(var_list_{i}).(M_.exo_names{j});
end
end
@@ -181,7 +181,7 @@ if options_.order==1
temp(i,j,:) = oo_.PosteriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(observable_name_requested_vars{i}).(M_.exo_names{j});
end
endo_index_varlist = strmatch(observable_name_requested_vars{i}, var_list_, 'exact');
- oo_ = posterior_analysis('decomposition', var_list_{endo_index_varlist}, 'ME', [], options_, M_, oo_);
+ oo_ = posterior_analysis('decomposition', var_list_{endo_index_varlist}, 'ME', [], options_, M_, oo_, estim_params_);
temp(i,j+1,:) = oo_.PosteriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(observable_name_requested_vars{i}).('ME');
end
title='Posterior mean variance decomposition (in percent) with measurement error';
@@ -219,7 +219,7 @@ if options_.order==1
if posterior
for i=1:NumberOfEndogenousVariables
for j=1:NumberOfExogenousVariables
- oo_ = posterior_analysis('conditional decomposition', var_list_{i}, M_.exo_names{j}, Steps, options_, M_, oo_);
+ oo_ = posterior_analysis('conditional decomposition', var_list_{i}, M_.exo_names{j}, Steps, options_, M_, oo_, estim_params_);
temp(i,j,:) = oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecomposition.Mean.(var_list_{i}).(M_.exo_names{j});
end
end
@@ -261,7 +261,7 @@ if options_.order==1
temp(i,j,:) = oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(observable_name_requested_vars{i}).(M_.exo_names{j});
end
endo_index_varlist = strmatch(observable_name_requested_vars{i}, var_list_, 'exact');
- oo_ = posterior_analysis('conditional decomposition', var_list_{endo_index_varlist}, 'ME', Steps, options_, M_, oo_);
+ oo_ = posterior_analysis('conditional decomposition', var_list_{endo_index_varlist}, 'ME', Steps, options_, M_, oo_, estim_params_);
temp(i,j+1,:) = oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(observable_name_requested_vars{i}).('ME');
end
title = 'Posterior mean conditional variance decomposition (in percent) with measurement error';
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 77ba20d6c4855664ed4259ed7a8245e5e7ba81b4..c9aa64c4e4922b7759d32c692356230890590d49 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -89,7 +89,6 @@ elseif exist([M_.fname '_steadystate.m'],'file')
else
options_.steadystate_flag = 0;
end
-options_.steadystate_partial = [];
options_.steadystate.nocheck = false;
% subset of the estimated deep parameters
@@ -400,7 +399,6 @@ options_.kalman.keep_kalman_algo_if_singularity_is_detected = false;
options_.diffuse_kalman_tol = 1e-6;
options_.use_univariate_filters_if_singularity_is_detected = 1;
options_.riccati_tol = 1e-6;
-options_.lik_algo = 1;
options_.lik_init = 1;
options_.load_mh_file = false;
options_.load_results_after_load_mh = false;
diff --git a/matlab/discretionary_policy/discretionary_policy_1.m b/matlab/discretionary_policy/discretionary_policy_1.m
index c2ddc000ec5db5617d80f19a88eacda93af3baa9..fee2055ad0a4030f27af472258cf3e97fece980c 100644
--- a/matlab/discretionary_policy/discretionary_policy_1.m
+++ b/matlab/discretionary_policy/discretionary_policy_1.m
@@ -1,16 +1,17 @@
-function [dr, info, M_, oo_]=discretionary_policy_1(Instruments, M_, options_, oo_)
+function [dr, info, params]=discretionary_policy_1(M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% Higher-level function for solving discretionary optimal policy
% INPUTS
-% - Instruments [cell] array containing instrument names
% - M_ [structure] Matlab's structure describing the model (M_).
% - options_ [structure] Matlab's structure describing the current options (options_).
-% - oo_ [structure] Matlab's structure containing the results (oo_).
+% - dr [struct] Decision rules for stochastic simulations.
+% - endo_steady_state [vector] steady state value for endogenous variables
+% - exo_steady_state [vector] steady state value for exogenous variables
+% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - dr [structure] Reduced form model.
% - info [integer] scalar or vector, error code.
-% - M_ [structure] Matlab's structure describing the model (M_).
-% - oo_ [structure] Matlab's structure containing the results (oo_).
+% - params [double] vector of potentially updated parameters
% Copyright © 2007-2020 Dynare Team
%
@@ -33,14 +34,12 @@ persistent Hold
info = 0;
-dr=oo_.dr; %initialize output argument
-
beta = get_optimal_policy_discount_factor(M_.params, M_.param_names);
%call steady_state_file if present to update parameters
if options_.steadystate_flag
% explicit steady state file
- [ys,M_.params,info] = evaluate_steady_state_file(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_, ...
+ [ys,M_.params,info] = evaluate_steady_state_file(endo_steady_state,[exo_steady_state; exo_det_steady_state],M_, ...
options_,false);
if info(1)
return;
@@ -48,6 +47,8 @@ if options_.steadystate_flag
else
ys=zeros(M_.endo_nbr,1);
end
+params=M_.params;
+
[U,Uy,W] = feval([M_.fname,'.objective.static'],zeros(M_.endo_nbr,1),[], M_.params);
if any(any(isnan(Uy)))
info = 64 ; %the derivatives of the objective function contain NaN
@@ -80,12 +81,7 @@ iyr0 = find(iyv(:)) ;
z = z(iyr0);
it_ = M_.maximum_lag + 1 ;
-if M_.exo_nbr == 0
- oo_.exo_steady_state = [] ;
-end
-
-[junk,jacobia_] = feval([M_.fname '.dynamic'],z, [zeros(size(oo_.exo_simul)) ...
- oo_.exo_det_simul], M_.params, ys, it_);
+[junk,jacobia_] = feval([M_.fname '.dynamic'],z,zeros(M_.exo_nbr+M_.exo_det_nbr,klen), M_.params, ys, it_);
if max(abs(junk))>options_.solve_tolf
info = 65; %the model must be written in deviation form and not have constant terms or have a steady state provided
return;
@@ -126,5 +122,4 @@ dr.ghu=G(dr.order_var,:);
if M_.maximum_endo_lag
Selection=M_.lead_lag_incidence(1,dr.order_var)>0;%select state variables
end
-dr.ghx=T(:,Selection);
-oo_.dr = dr;
+dr.ghx=T(:,Selection);
\ No newline at end of file
diff --git a/matlab/dsge_conditional_likelihood_1.m b/matlab/dsge_conditional_likelihood_1.m
index 55eefe99f54539cdc07e4588cae3b572c28df521..9057ff0afd2dbab5421737e09b172d02c89c5255 100644
--- a/matlab/dsge_conditional_likelihood_1.m
+++ b/matlab/dsge_conditional_likelihood_1.m
@@ -61,9 +61,9 @@ end
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
-% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
-[T, R, SteadyState, info, Model, DynareResults] = ...
- dynare_resolve(Model, DynareOptions, DynareResults, 'restrict');
+% Linearize the model around the deterministic steadystate and extract the matrices of the state equation (T and R).
+[T, R, SteadyState, info,DynareResults.dr, Model.params] = ...
+ dynare_resolve(Model, DynareOptions, DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state, 'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 7c5bbac7c626ff2f91d41c96b7db15c520a0ccf0..3110a02f959cedbc9b1ffbb1ba3444fc113904c7 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -192,7 +192,7 @@ is_restrict_state_space = true;
if DynareOptions.occbin.likelihood.status
occbin_options = set_occbin_options(DynareOptions, Model);
if occbin_options.opts_simul.restrict_state_space
- [T,R,SteadyState,info,Model,DynareResults,TTx,RRx,CCx, T0, R0] = ...
+ [T,R,SteadyState,info,DynareResults.dr, Model.params,TTx,RRx,CCx, T0, R0] = ...
occbin.dynare_resolve(Model,DynareOptions,DynareResults,[],'restrict');
else
is_restrict_state_space = false;
@@ -202,7 +202,7 @@ if DynareOptions.occbin.likelihood.status
DynareResults.dr.restrict_columns = BayesInfo.smoother_restrict_columns;
% Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
- [T,R,SteadyState,info,Model,DynareOptions,DynareResults,TTx,RRx,CCx, T0, R0] = ...
+ [T,R,SteadyState,info,Model,DynareResults.dr, Model.params,TTx,RRx,CCx, T0, R0] = ...
occbin.dynare_resolve(Model,DynareOptions,DynareResults);
DynareResults.dr.restrict_var_list = oldoo.restrict_var_list;
@@ -213,7 +213,7 @@ if DynareOptions.occbin.likelihood.status
occbin_.info= {DynareOptions, DynareResults, Model, occbin_options, TTx, RRx, CCx,T0,R0};
else
% Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
- [T,R,SteadyState,info,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+ [T,R,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict');
occbin_.status = false;
end
@@ -542,7 +542,7 @@ if analytic_derivation
AHess = [];
iv = DynareResults.dr.restrict_var_list;
if nargin<10 || isempty(derivatives_info)
- [A,B,nou,nou,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
+ [A,B,nou,nou,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
if ~isempty(EstimatedParameters.var_exo)
indexo=EstimatedParameters.var_exo(:,1);
else
diff --git a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
index 5918b75e475f62a37f8b7490a3a6993b8e8194e8..e70663ac4abc92f5876362d552c6fc09ad99f30d 100644
--- a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
+++ b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
@@ -132,7 +132,7 @@ for file = 1:NumberOfDrawsFiles
dr = temp.pdraws{linee,2};
else
M_=set_parameters_locally(M_,temp.pdraws{linee,1});
- [dr,info,M_,oo_] = compute_decision_rules(M_,options_,oo_);
+ [dr,info,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
end
M_ = set_measurement_errors(temp.pdraws{linee,1},temp.estim_params_,M_);
diff --git a/matlab/dsge_simulated_theoretical_covariance.m b/matlab/dsge_simulated_theoretical_covariance.m
index e8c4747ad2e5da10a15b8f4b97f6ca24884355c2..3fb978ada543470923f2108f91c0c6825592193d 100644
--- a/matlab/dsge_simulated_theoretical_covariance.m
+++ b/matlab/dsge_simulated_theoretical_covariance.m
@@ -127,7 +127,7 @@ for file = 1:NumberOfDrawsFiles
dr = temp.pdraws{linee,2};
else
M_=set_parameters_locally(M_,temp.pdraws{linee,1});
- [dr,info,M_,oo_] = compute_decision_rules(M_,options_,oo_);
+ [dr,info,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
end
if ~options_.pruning
tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
diff --git a/matlab/dsge_simulated_theoretical_variance_decomposition.m b/matlab/dsge_simulated_theoretical_variance_decomposition.m
index 305dc76d227d3772a6bb01a809c729c1a68ad134..38863c0ac038aad84e3588d1d6e01cd80ed978c2 100644
--- a/matlab/dsge_simulated_theoretical_variance_decomposition.m
+++ b/matlab/dsge_simulated_theoretical_variance_decomposition.m
@@ -137,7 +137,7 @@ for file = 1:NumberOfDrawsFiles
dr = temp.pdraws{linee,2};
else
M_=set_parameters_locally(M_,temp.pdraws{linee,1});
- [dr,info,M_,oo_] = compute_decision_rules(M_,options_,oo_);
+ [dr,info,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
end
if file==1 && linee==1
[tmp, stationary_vars] = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
diff --git a/matlab/dsge_var_likelihood.m b/matlab/dsge_var_likelihood.m
index 2481c45cb587e0dacc630792d49b15f6bb7b61b8..910ef5ad373d523573cb551c1a59f90c6b673a3b 100644
--- a/matlab/dsge_var_likelihood.m
+++ b/matlab/dsge_var_likelihood.m
@@ -136,7 +136,7 @@ end
% Solve the Dsge model and get the matrices of the reduced form solution. T and R are the matrices of the
% state equation
-[T,R,SteadyState,info] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+[T,R,SteadyState,info] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 704c2b9c65f3a2786cbeeda406978c182d94707c..1e216266b40382372c5a57a0af884683b46538b0 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -93,6 +93,7 @@ if options_.order > 1
end
end
+%% set objective function
if ~options_.dsge_var
if options_.particle.status
objective_function = str2func('non_linear_dsge_likelihood');
@@ -133,16 +134,12 @@ if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(esti
end
data = dataset_.data;
-rawdata = dataset_info.rawdata;
data_index = dataset_info.missing.aindex;
missing_value = dataset_info.missing.state;
% Set number of observations
gend = dataset_.nobs;
-% Set the number of observed variables.
-n_varobs = length(options_.varobs);
-
% Get the number of parameters to be estimated.
nvx = estim_params_.nvx; % Variance of the structural innovations (number of parameters).
nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters).
@@ -151,13 +148,10 @@ ncn = estim_params_.ncn; % Covariance of the measurement innovations (number of
np = estim_params_.np ; % Number of deep parameters.
nx = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated.
-dr = oo_.dr;
-
if ~isempty(estim_params_)
M_ = set_all_parameters(xparam1,estim_params_,M_);
end
-
%% perform initial estimation checks;
try
oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);
@@ -175,6 +169,7 @@ catch % if check fails, provide info on using calibration if present
rethrow(e);
end
+%% Run smoother if no estimation or mode-finding are requested
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
if options_.order==1 && ~options_.particle.status
if options_.smoother
@@ -195,14 +190,14 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.
end
else
if options_.occbin.smoother.status
- [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info);
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info);
if oo_.occbin.smoother.error_flag(1)==0
[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
else
fprintf('\nOccbin: smoother did not succeed. No results will be written to oo_.\n')
end
else
- [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
end
end
@@ -229,71 +224,71 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
for optim_iter = 1:length(optimizer_vec)
current_optimizer = optimizer_vec{optim_iter};
- [xparam1, fval, exitflag, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
- fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
+ [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
+ fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
if isnumeric(current_optimizer)
if current_optimizer==5
- newratflag = new_rat_hess_info.newratflag;
- new_rat_hess_info = new_rat_hess_info.new_rat_hess_info;
+ newratflag = new_rat_hess_info.newratflag;
+ new_rat_hess_info = new_rat_hess_info.new_rat_hess_info;
elseif current_optimizer==6 %save scaling factor
- save([M_.dname filesep 'Output' filesep M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
- options_.mh_jscale = Scale;
- bayestopt_.jscale(:) = options_.mh_jscale;
- end
+ save([M_.dname filesep 'Output' filesep M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
+ options_.mh_jscale = Scale;
+ bayestopt_.jscale(:) = options_.mh_jscale;
+ end
end
if ~isnumeric(current_optimizer) || ~isequal(current_optimizer,6) %always already computes covariance matrix
- if options_.cova_compute == 1 %user did not request covariance not to be computed
- if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood')
- ana_deriv_old = options_.analytic_derivation;
- options_.analytic_derivation = 2;
- [~,~,~,~,hh] = feval(objective_function,xparam1, ...
- dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
- options_.analytic_derivation = ana_deriv_old;
+ if options_.cova_compute == 1 %user did not request covariance not to be computed
+ if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood')
+ ana_deriv_old = options_.analytic_derivation;
+ options_.analytic_derivation = 2;
+ [~,~,~,~,hh] = feval(objective_function,xparam1, ...
+ dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
+ options_.analytic_derivation = ana_deriv_old;
elseif ~isnumeric(current_optimizer) || ~(isequal(current_optimizer,5) && newratflag~=1 && strcmp(func2str(objective_function),'dsge_likelihood'))
- % enter here if i) not mode_compute_5, ii) if mode_compute_5 and newratflag==1;
- % with flag==0 or 2 and dsge_likelihood, we force to use
- % the hessian from outer product gradient of optimizer 5 below
- if options_.hessian.use_penalized_objective
- penalized_objective_function = str2func('penalty_objective_function');
- hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
- else
- hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
- end
- hh = reshape(hh, nx, nx);
+ % enter here if i) not mode_compute_5, ii) if mode_compute_5 and newratflag==1;
+ % with flag==0 or 2 and dsge_likelihood, we force to use
+ % the hessian from outer product gradient of optimizer 5 below
+ if options_.hessian.use_penalized_objective
+ penalized_objective_function = str2func('penalty_objective_function');
+ hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
+ else
+ hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
+ end
+ hh = reshape(hh, nx, nx);
elseif isnumeric(current_optimizer) && isequal(current_optimizer,5)
- % other numerical hessian options available with optimizer
- % 5 and dsge_likelihood
- %
- % if newratflag == 0
- % compute outer product gradient of optimizer 5
- %
- % if newratflag == 2
- % compute 'mixed' outer product gradient of optimizer 5
- % with diagonal elements computed with numerical second order derivatives
- %
- % uses univariate filters, so to get max # of available
- % densities for outer product gradient
- kalman_algo0 = options_.kalman_algo;
- compute_hessian = 1;
- if ~((options_.kalman_algo == 2) || (options_.kalman_algo == 4))
- options_.kalman_algo=2;
- if options_.lik_init == 3
- options_.kalman_algo=4;
- end
- elseif newratflag==0 % hh already contains outer product gradient with univariate filter
- compute_hessian = 0;
- end
- if compute_hessian
- crit = options_.newrat.tolerance.f;
- newratflag = newratflag>0;
- hh = reshape(mr_hessian(xparam1,objective_function,fval,newratflag,crit,new_rat_hess_info,[bounds.lb bounds.ub],bayestopt_.p2,0,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx);
+ % other numerical hessian options available with optimizer
+ % 5 and dsge_likelihood
+ %
+ % if newratflag == 0
+ % compute outer product gradient of optimizer 5
+ %
+ % if newratflag == 2
+ % compute 'mixed' outer product gradient of optimizer 5
+ % with diagonal elements computed with numerical second order derivatives
+ %
+ % uses univariate filters, so to get max # of available
+ % densities for outer product gradient
+ kalman_algo0 = options_.kalman_algo;
+ compute_hessian = 1;
+ if ~((options_.kalman_algo == 2) || (options_.kalman_algo == 4))
+ options_.kalman_algo=2;
+ if options_.lik_init == 3
+ options_.kalman_algo=4;
end
- options_.kalman_algo = kalman_algo0;
+ elseif newratflag==0 % hh already contains outer product gradient with univariate filter
+ compute_hessian = 0;
end
+ if compute_hessian
+ crit = options_.newrat.tolerance.f;
+ newratflag = newratflag>0;
+ hh = reshape(mr_hessian(xparam1,objective_function,fval,newratflag,crit,new_rat_hess_info,[bounds.lb bounds.ub],bayestopt_.p2,0,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx);
+ end
+ options_.kalman_algo = kalman_algo0;
end
end
- parameter_names = bayestopt_.name;
+ end
+ parameter_names = bayestopt_.name;
end
if options_.cova_compute || current_optimizer==5 || current_optimizer==6
save([M_.dname filesep 'Output' filesep M_.fname '_mode.mat'],'xparam1','hh','parameter_names','fval');
@@ -306,6 +301,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute
check_hessian_at_the_mode(hh, xparam1, M_, estim_params_, options_, bounds);
end
+%% create mode_check_plots
if options_.mode_check.status && ~options_.mh_posterior_mode_estimation
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0;
@@ -440,7 +436,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
end
%% Estimation of the marginal density from the Mh draws:
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
- [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
+ [~,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
% Store posterior statistics by parameter name
oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, prior_dist_names);
if ~options_.nograph
@@ -476,7 +472,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if error_flag
error('%s: I cannot compute the posterior IRFs!',dispString)
end
- PosteriorIRF('posterior',dispString);
+ oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString);
end
if options_.moments_varendo
if error_flag
@@ -501,14 +497,14 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
end
end
end
- oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
+ oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
end
if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
if error_flag
error('%s: I cannot compute the posterior statistics!',dispString)
end
if options_.order==1 && ~options_.particle.status
- prior_posterior_statistics('posterior',dataset_,dataset_info,dispString); %get smoothed and filtered objects and forecasts
+ oo_=prior_posterior_statistics('posterior',dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_,dispString); %get smoothed and filtered objects and forecasts
else
error('%s: Particle Smoothers are not yet implemented.',dispString)
end
@@ -519,6 +515,8 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
M_ = set_all_parameters(xparam1,estim_params_,M_);
end
+else
+ M_ = set_all_parameters(xparam1,estim_params_,M_); %set to posterior mode
end
if options_.particle.status
@@ -527,238 +525,16 @@ if options_.particle.status
return
end
+%Run and store classical smoother if needed
if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file)) ...
|| ~options_.smoother ) && ~options_.partial_information % to be fixed
%% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
- if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
- [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_);
- if ismember(info(1),[303,304,306])
- fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n')
- else
- updated_variables = atT*nan;
- measurement_error=[];
- ys = oo_.dr.ys;
- trend_coeff = zeros(length(options_.varobs_id),1);
- bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf);
- options_nk=options_.nk;
- options_.nk=[]; %unset options_.nk and reset it later
- [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff);
- options_.nk=options_nk;
- end
- else
- if options_.occbin.smoother.status
- [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info);
- if oo_.occbin.smoother.error_flag(1)==0
- [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
- else
- fprintf('\nOccbin: smoother did not succeed. No results will be written to oo_.\n')
- end
- else
- [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_);
- [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
- end
- [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
- end
- if ~options_.nograph
- [nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);
- if ~exist([M_.dname '/graphs'],'dir')
- mkdir(M_.dname,'graphs');
- end
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_SmoothedShocks.tex'],'w');
- fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
- fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
- fprintf(fidTeX,' \n');
- end
- for plt = 1:nbplt
- fh = dyn_figure(options_.nodisplay,'Name','Smoothed shocks');
- NAMES = [];
- if options_.TeX, TeXNAMES = []; end
- nstar0=min(nstar,M_.exo_nbr-(plt-1)*nstar);
- if gend==1
- marker_string{1,1}='-ro';
- marker_string{2,1}='-ko';
- else
- marker_string{1,1}='-r';
- marker_string{2,1}='-k';
- end
- for i=1:nstar0
- k = (plt-1)*nstar+i;
- subplot(nr,nc,i);
- plot([1 gend],[0 0],marker_string{1,1},'linewidth',.5)
- hold on
- plot(1:gend,innov(k,:),marker_string{2,1},'linewidth',1)
- hold off
- name = M_.exo_names{k};
- if ~isempty(options_.XTick)
- set(gca,'XTick',options_.XTick)
- set(gca,'XTickLabel',options_.XTickLabel)
- end
- if gend>1
- xlim([1 gend])
- end
- if options_.TeX
- title(['$' M_.exo_names_tex{k} '$'],'Interpreter','latex')
- else
- title(name,'Interpreter','none')
- end
- end
- dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_SmoothedShocks' int2str(plt)],options_.nodisplay,options_.graph_format);
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fprintf(fidTeX,'\\begin{figure}[H]\n');
- fprintf(fidTeX,'\\centering \n');
- fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_SmoothedShocks%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
- fprintf(fidTeX,'\\caption{Smoothed shocks.}');
- fprintf(fidTeX,'\\label{Fig:SmoothedShocks:%s}\n',int2str(plt));
- fprintf(fidTeX,'\\end{figure}\n');
- fprintf(fidTeX,'\n');
- end
- end
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fprintf(fidTeX,'\n');
- fprintf(fidTeX,'%% End of TeX file.\n');
- fclose(fidTeX);
- end
- end
- if nvn
- number_of_plots_to_draw = 0;
- index = [];
- for obs_iter=1:n_varobs
- if max(abs(measurement_error(obs_iter,:))) > options_.ME_plot_tol;
- number_of_plots_to_draw = number_of_plots_to_draw + 1;
- index = cat(1,index,obs_iter);
- end
- end
- if ~options_.nograph
- [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw);
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_SmoothedObservationErrors.tex'],'w');
- fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
- fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
- fprintf(fidTeX,' \n');
- end
- for plt = 1:nbplt
- fh = dyn_figure(options_.nodisplay,'Name','Smoothed observation errors');
- nstar0=min(nstar,number_of_plots_to_draw-(plt-1)*nstar);
- if gend==1
- marker_string{1,1}='-ro';
- marker_string{2,1}='-ko';
- else
- marker_string{1,1}='-r';
- marker_string{2,1}='-k';
- end
- for i=1:nstar0
- k = (plt-1)*nstar+i;
- subplot(nr,nc,i);
- plot([1 gend],[0 0],marker_string{1,1},'linewidth',.5)
- hold on
- plot(1:gend,measurement_error(index(k),:),marker_string{2,1},'linewidth',1)
- hold off
- name = options_.varobs{index(k)};
- if gend>1
- xlim([1 gend])
- end
- if ~isempty(options_.XTick)
- set(gca,'XTick',options_.XTick)
- set(gca,'XTickLabel',options_.XTickLabel)
- end
- if options_.TeX
- idx = strmatch(options_.varobs{index(k)}, M_.endo_names, 'exact');
- title(['$' M_.endo_names_tex{idx} '$'],'Interpreter','latex')
- else
- title(name,'Interpreter','none')
- end
- end
- dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_SmoothedObservationErrors' int2str(plt)],options_.nodisplay,options_.graph_format);
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fprintf(fidTeX,'\\begin{figure}[H]\n');
- fprintf(fidTeX,'\\centering \n');
- fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_SmoothedObservationErrors%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
- fprintf(fidTeX,'\\caption{Smoothed observation errors.}');
- fprintf(fidTeX,'\\label{Fig:SmoothedObservationErrors:%s}\n',int2str(plt));
- fprintf(fidTeX,'\\end{figure}\n');
- fprintf(fidTeX,'\n');
- end
- end
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fprintf(fidTeX,'\n');
- fprintf(fidTeX,'%% End of TeX file.\n');
- fclose(fidTeX);
- end
- end
- end
- %%
- %% Historical and smoothed variabes
- %%
- if ~options_.nograph
- [nbplt,nr,nc,lr,lc,nstar] = pltorg(n_varobs);
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_HistoricalAndSmoothedVariables.tex'],'w');
- fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
- fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
- fprintf(fidTeX,' \n');
- end
- for plt = 1:nbplt
- fh = dyn_figure(options_.nodisplay,'Name','Historical and smoothed variables');
- NAMES = [];
- nstar0=min(nstar,n_varobs-(plt-1)*nstar);
- if gend==1
- marker_string{1,1}='-ro';
- marker_string{2,1}='--ko';
- else
- marker_string{1,1}='-r';
- marker_string{2,1}='--k';
- end
- for i=1:nstar0
- k = (plt-1)*nstar+i;
- subplot(nr,nc,i);
- plot(1:gend,yf(k,:),marker_string{1,1},'linewidth',1)
- hold on
- plot(1:gend,rawdata(:,k),marker_string{2,1},'linewidth',1)
- hold off
- name = options_.varobs{k};
- if ~isempty(options_.XTick)
- set(gca,'XTick',options_.XTick)
- set(gca,'XTickLabel',options_.XTickLabel)
- end
- if gend>1
- xlim([1 gend])
- end
- if options_.TeX
- idx = strmatch(options_.varobs{k}, M_.endo_names,'exact');
- title(['$' M_.endo_names_tex{idx} '$'],'Interpreter','latex')
- else
- title(name,'Interpreter','none')
- end
- end
- dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_HistoricalAndSmoothedVariables' int2str(plt)],options_.nodisplay,options_.graph_format);
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fprintf(fidTeX,'\\begin{figure}[H]\n');
- fprintf(fidTeX,'\\centering \n');
- fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_HistoricalAndSmoothedVariables%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
- fprintf(fidTeX,'\\caption{Historical and smoothed variables.}');
- fprintf(fidTeX,'\\label{Fig:HistoricalAndSmoothedVariables:%s}\n',int2str(plt));
- fprintf(fidTeX,'\\end{figure}\n');
- fprintf(fidTeX,'\n');
- end
- end
- if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- fprintf(fidTeX,'\n');
- fprintf(fidTeX,'%% End of TeX file.\n');
- fclose(fidTeX);
- end
- end
+ oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
end
-
if options_.forecast > 0 && options_.mh_replic == 0 && ~options_.load_mh_file
oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info);
end
-if np > 0
- pindx = estim_params_.param_vals(:,1);
- save([M_.dname filesep 'Output' filesep M_.fname '_pindx.mat'] ,'pindx');
-end
-
%reset qz_criterium
options_.qz_criterium=qz_criterium_old;
diff --git a/matlab/dynare_resolve.m b/matlab/dynare_resolve.m
index b3bf6ed1a5cb69f06cf09dacffd13f49b68fd306..8249293470b2e24e88f1db588d4cdb10db527635 100644
--- a/matlab/dynare_resolve.m
+++ b/matlab/dynare_resolve.m
@@ -1,11 +1,14 @@
-function [A,B,ys,info,M_,oo_] = dynare_resolve(M_,options_,oo_,mode)
-% function [A,B,ys,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,mode)
+function [A,B,ys,info,dr,params] = dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,mode)
+% function [A,B,ys,info,dr,params] = dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,mode)
% Computes the linear approximation and the matrices A and B of the transition equation.
%
% Inputs:
% - M_ [structure] Matlab's structure describing the model
% - options_ [structure] Matlab's structure containing the options
-% - oo_ [structure] Matlab's structure containing the results
+% - dr [structure] Reduced form model.
+% - endo_steady_state [vector] steady state value for endogenous variables
+% - exo_steady_state [vector] steady state value for exogenous variables
+% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - mode [string] if provided, use restricted state space
%
% Outputs:
@@ -13,8 +16,8 @@ function [A,B,ys,info,M_,oo_] = dynare_resolve(M_,options_,oo_,mode)
% - B [double] shock impact matrix (potentially for restricted state space)
% - ys [double] vector of steady state values
% - info [double] 4 by 1 vector with exit flag and information
-% - M_ [structure] Matlab's structure describing the model
-% - oo_ [structure] Matlab's structure containing the results
+% - dr [structure] Reduced form model.
+% - params [double] vector of potentially updated parameters
% Copyright © 2001-2023 Dynare Team
%
@@ -33,7 +36,7 @@ function [A,B,ys,info,M_,oo_] = dynare_resolve(M_,options_,oo_,mode)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-[dr,info,M_,oo_] =compute_decision_rules(M_,options_,oo_);
+[dr,info,params] =compute_decision_rules(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
if info(1) > 0
A = [];
@@ -47,23 +50,23 @@ if info(1) > 0
end
switch nargin
- case 3
+ case 6
endo_nbr = M_.endo_nbr;
nstatic = M_.nstatic;
nspred = M_.nspred;
iv = (1:endo_nbr)';
- ic = [ nstatic+(1:nspred) endo_nbr+(1:size(oo_.dr.ghx,2)-nspred) ]';
- case 4
- iv = oo_.dr.restrict_var_list;
- ic = oo_.dr.restrict_columns;
+ ic = [ nstatic+(1:nspred) endo_nbr+(1:size(dr.ghx,2)-nspred) ]';
+ case 7
+ iv = dr.restrict_var_list;
+ ic = dr.restrict_columns;
otherwise
error('dynare_resolve:: Error in the calling sequence!')
end
if nargout==1
- A = kalman_transition_matrix(oo_.dr,iv,ic);
+ A = kalman_transition_matrix(dr,iv,ic);
return
end
-[A,B] = kalman_transition_matrix(oo_.dr,iv,ic);
-ys = oo_.dr.ys;
+[A,B] = kalman_transition_matrix(dr,iv,ic);
+ys = dr.ys;
diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m
index 4a3389d225f08bc92122f6ed692843e306aa7fc2..1b1cfe30bbd3fcf9682478a9f1c64fe933624505 100644
--- a/matlab/dynare_sensitivity.m
+++ b/matlab/dynare_sensitivity.m
@@ -146,7 +146,7 @@ if M_.exo_nbr==0
error('dynare_sensitivity does not support having no varexo in the model. As a workaround you could define a dummy exogenous variable.')
end
-[make,my,day,punk,M_,oo_] = dynare_resolve(M_,options_,oo_);
+[make,my,day,punk,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
options_gsa = set_default_option(options_gsa,'identification',0);
if options_gsa.identification
@@ -410,9 +410,9 @@ if options_gsa.rmse
end
end
- prior_posterior_statistics('gsa',dataset_, dataset_info,'gsa::mcmc');
+ oo_=prior_posterior_statistics('gsa',dataset_, dataset_info,M_,oo_,options_,estim_params_,bayestopt_,'gsa::mcmc');
if options_.bayesian_irf
- PosteriorIRF('gsa','gsa::mcmc');
+ oo_=PosteriorIRF('gsa',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,'gsa::mcmc');
end
options_gsa.load_rmse=0;
% else
diff --git a/matlab/endogenous_prior_restrictions.m b/matlab/endogenous_prior_restrictions.m
index 8549fd63008da8441d0e02e82f5ae2dbff264a5b..71877a30af8ead7ac9cf94649aca54aa2df4de0e 100644
--- a/matlab/endogenous_prior_restrictions.m
+++ b/matlab/endogenous_prior_restrictions.m
@@ -49,7 +49,7 @@ if ~isempty(endo_prior_restrictions.irf)
end
varlist = Model.endo_names(DynareResults.dr.order_var);
if isempty(T)
- [T,R,SteadyState,infox,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
+ [T,R,SteadyState,infox,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
else % check if T and R are given in the restricted form!!!
if size(T,1)<length(varlist)
varlist = varlist(DynareResults.dr.restrict_var_list);
@@ -63,7 +63,7 @@ if ~isempty(endo_prior_restrictions.irf)
end
if ~varlistok
varlist = Model.endo_names(DynareResults.dr.order_var);
- [T,R,SteadyState,infox,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
+ [T,R,SteadyState,infox,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
end
end
NT=1;
diff --git a/matlab/ep/extended_path_initialization.m b/matlab/ep/extended_path_initialization.m
index 522f2c664e4afa8639f4f64b665f2075f2196eec..6f41c4d484eb82cc920ee87b8f03ee6b00fa1d32 100644
--- a/matlab/ep/extended_path_initialization.m
+++ b/matlab/ep/extended_path_initialization.m
@@ -78,7 +78,7 @@ dr = struct();
if ep.init
DynareOptions.order = 1;
DynareResults.dr=set_state_space(dr,DynareModel,DynareOptions);
- [dr,Info,DynareResults] = resol(0,DynareModel,DynareOptions,DynareResults);
+ [DynareResults.dr,Info,DynareModel.params] = resol(0,DynareModel,DynareOptions,DynareResults.dr,DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
end
% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
@@ -106,7 +106,7 @@ if pfm.hybrid_order
DynareResults.dr = set_state_space(DynareResults.dr, DynareModel, DynareOptions);
options = DynareOptions;
options.order = pfm.hybrid_order;
- pfm.dr = resol(0, DynareModel, options, DynareResults);
+ [pfm.dr, DynareModel.params] = resol(0, DynareModel, options, DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
else
pfm.dr = [];
end
diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m
index b0fb55501e60bfa4232f91a7e9f29d1db4739748..6162a3b1716e09760ab4f74170adae773c3ad695 100644
--- a/matlab/evaluate_smoother.m
+++ b/matlab/evaluate_smoother.m
@@ -117,11 +117,11 @@ if options_.occbin.smoother.status
bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf);
end
else
- [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = ...
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = ...
occbin.DSGE_smoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info);
end
else
- [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = ...
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = ...
DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_);
end
if ~(options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter)
diff --git a/matlab/evaluate_steady_state_file.m b/matlab/evaluate_steady_state_file.m
index ee312d881530fee13559b1a3321192152a388a8d..d13ac7091ca84930c6ed2bce8b9b3504750220d2 100644
--- a/matlab/evaluate_steady_state_file.m
+++ b/matlab/evaluate_steady_state_file.m
@@ -114,21 +114,4 @@ if steady_state_checkflag
info(1) = 22;
return
end
-elseif ~isempty(options.steadystate_partial)
- ssvar = options.steadystate_partial.ssvar;
- nov = length(ssvar);
- indv = zeros(nov,1);
- for i = 1:nov
- indv(i) = strmatch(ssvar(i), M.endo_names, 'exact');
- end
- ys = dynare_solve(@restricted_steadystate, ys(indv), options.steady.maxit, options.dynatol.f, options.dynatol.x, options, exo_ss, indv);
-end
-
-function [sR,sG] = restricted_steadystate(y,x,indx)
-global options_ M_ oo_
-inde = options_.steadystate_partial.sseqn;
-ss = oo_.steady_state;
-ss(indx) = y;
-[R,G] = feval([M_.fname '.static'], ss, x, M_.params);
-sR = R(inde);
-sG = G(inde,indx);
+end
\ No newline at end of file
diff --git a/matlab/forcst2.m b/matlab/forcst2.m
deleted file mode 100644
index 4c2c5b5b25826484bd2cc5fc12be7e4cd9724400..0000000000000000000000000000000000000000
--- a/matlab/forcst2.m
+++ /dev/null
@@ -1,67 +0,0 @@
-function yf=forcst2(y0,horizon,dr,n)
-% function yf=forcst2(y0,horizon,dr,n)
-%
-% computes forecasts based on first order model solution, given shocks
-% drawn from the shock distribution, but not including measurement error
-% Inputs:
-% - y0 [endo_nbr by maximum_endo_lag] matrix of starting values
-% - dr [structure] structure with Dynare decision rules
-% - e [horizon by exo_nbr] matrix with shocks
-% - n [scalar] number of repetitions
-%
-% Outputs:
-% - yf [horizon+ykmin_ by endo_nbr by n] array of forecasts
-%
-% Copyright © 2008-2017 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-
-global M_ options_
-
-Sigma_e_ = M_.Sigma_e;
-endo_nbr = M_.endo_nbr;
-exo_nbr = M_.exo_nbr;
-ykmin_ = M_.maximum_endo_lag;
-
-k1 = [ykmin_:-1:1];
-k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
-k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
-
-% eliminate shocks with 0 variance
-i_exo_var = setdiff([1:exo_nbr],find(diag(Sigma_e_) == 0));
-nxs = length(i_exo_var);
-
-chol_S = chol(Sigma_e_(i_exo_var,i_exo_var));
-
-if ~isempty(Sigma_e_)
- e = randn(nxs,n,horizon);
-end
-
-B1 = dr.ghu(:,i_exo_var)*chol_S';
-
-yf = zeros(endo_nbr,horizon+ykmin_,n);
-yf(:,1:ykmin_,:,:) = repmat(y0,[1,1,n]);
-
-j = ykmin_*endo_nbr;
-for i=ykmin_+(1:horizon)
- tempx1 = reshape(yf(:,k1,:),[j,n]);
- tempx = tempx1(k2,:);
- yf(:,i,:) = dr.ghx*tempx+B1*squeeze(e(:,:,i-ykmin_));
- k1 = k1+1;
-end
-
-yf(dr.order_var,:,:) = yf;
-yf=permute(yf,[2 1 3]);
diff --git a/matlab/forcst2a.m b/matlab/forcst2a.m
deleted file mode 100644
index 2b602f30758ba045f2c1fd3652525ed0bb3fce12..0000000000000000000000000000000000000000
--- a/matlab/forcst2a.m
+++ /dev/null
@@ -1,50 +0,0 @@
-function yf=forcst2a(y0,dr,e)
-% function yf=forcst2a(y0,dr,e)
-% computes forecasts based on first order model solution, assuming the absence of shocks
-% Inputs:
-% - y0 [endo_nbr by maximum_endo_lag] matrix of starting values
-% - dr [structure] structure with Dynare decision rules
-% - e [horizon by exo_nbr] matrix with shocks
-%
-% Outputs:
-% - yf [horizon+maximum_endo_lag,endo_nbr] matrix of forecasts
-%
-% Copyright © 2008-2017 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-
-global M_ options_
-
-endo_nbr = M_.endo_nbr;
-ykmin_ = M_.maximum_endo_lag;
-
-horizon = size(e,1);
-
-k1 = [ykmin_:-1:1];
-k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
-k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
-
-yf = zeros(horizon+ykmin_,endo_nbr);
-yf(1:ykmin_,:) = y0';
-
-j = ykmin_*endo_nbr;
-for i=ykmin_+(1:horizon)
- tempx = yf(k1,:)';
- yf(i,:) = tempx(k2)'*dr.ghx';
- k1 = k1+1;
-end
-
-yf(:,dr.order_var) = yf;
diff --git a/matlab/get_perturbation_params_derivs_numerical_objective.m b/matlab/get_perturbation_params_derivs_numerical_objective.m
index 8b8c9591f9a01dab89a1504f27232301388fb603..d9a4bbcd3560be0f892a13857789a4c6e52b40e4 100644
--- a/matlab/get_perturbation_params_derivs_numerical_objective.m
+++ b/matlab/get_perturbation_params_derivs_numerical_objective.m
@@ -51,7 +51,7 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
%% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters
M = set_all_parameters(params,estim_params,M);
-[~,info,M,oo] = compute_decision_rules(M,options,oo);
+[oo.dr,info,M.params] = compute_decision_rules(M,options,oo.dr, oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state);
Sigma_e = M.Sigma_e;
if info(1) > 0
diff --git a/matlab/gsa/map_calibration.m b/matlab/gsa/map_calibration.m
index 2244d5d2fb2c7021e448cc9a92a11ee099584e6d..c1755a4742e5d7a7064efd7dd08b4468a05f297b 100644
--- a/matlab/gsa/map_calibration.m
+++ b/matlab/gsa/map_calibration.m
@@ -102,9 +102,9 @@ if init
for j=1:Nsam
Model = set_all_parameters(lpmat(j,:)',EstimatedParameters,Model);
if nbr_moment_restrictions
- [Tt,Rr,SteadyState,info,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
+ [Tt,Rr,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
else
- [Tt,Rr,SteadyState,info,Model,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+ [Tt,Rr,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict');
end
if info(1)==0
[info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults);
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index 4a3efcfe8cfc52e80820a66b086d0151db2c23c3..ea2a5de62cc99373177fe134a7f4a68ea0e20557 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -242,9 +242,9 @@ if fload==0
%try stoch_simul([]);
try
if ~ isempty(options_.endogenous_prior_restrictions.moment)
- [Tt,Rr,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_);
+ [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
- [Tt,Rr,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
+ [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
end
infox(j,1)=info(1);
if infox(j,1)==0 && ~exist('T','var')
@@ -395,7 +395,7 @@ else
for j=1:ntrans
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
%stoch_simul([]);
- [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
+ [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
if ~exist('T','var')
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans);
end
@@ -570,7 +570,7 @@ if length(iunstable)>0 || length(iwrong)>0
end
M_ = set_all_parameters(x0,estim_params_,M_);
- [oo_.dr,info,M_,oo_] = resol(0,M_,options_,oo_);
+ [oo_.dr,info,M_.params] = resol(0,M_,options_,oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
else
disp('All parameter values in the specified ranges are not acceptable!')
x0=[];
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index 308be40b6a274aec4e75790e97efcd1648a420d6..1feeb6552b77d076fa72f52e71a83e455addbc54 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -134,7 +134,7 @@ error_indicator.identification_minimal=0;
error_indicator.identification_spectrum=0;
%Compute linear approximation and fill dr structure
-[oo_.dr,info,M_,oo_] = compute_decision_rules(M_,options_,oo_);
+[oo_.dr,info,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if info(1) == 0 %no errors in solution
% Compute parameter Jacobians for identification analysis
diff --git a/matlab/identification_numerical_objective.m b/matlab/identification_numerical_objective.m
index c703d2b84ad8119ebb7bc38507ec291e2966e252..172f1f6f2e021bf92ebbd9576bac5c57d5ae54a8 100644
--- a/matlab/identification_numerical_objective.m
+++ b/matlab/identification_numerical_objective.m
@@ -76,7 +76,7 @@ else
end
%% compute Kalman transition matrices and steady state with updated parameters
-[~,info,M,oo] = compute_decision_rules(M,options,oo);
+[oo.dr,info,M.params] = compute_decision_rules(M,options,oo.dr, oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state);
options = rmfield(options,'options_ident');
pruned = pruned_state_space_system(M, options, oo.dr, indvar, nlags, useautocorr, 0);
diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m
index 07ed86af5583ddf0f1e5e2c3b837e8cbd9c7aec3..e202943dc475d8bd487907b267a076db7f3820b0 100644
--- a/matlab/imcforecast.m
+++ b/matlab/imcforecast.m
@@ -126,7 +126,7 @@ if estimated_model
qz_criterium_old=options_.qz_criterium;
options_=select_qz_criterium_value(options_);
options_smoothed_state_uncertainty_old = options_.smoothed_state_uncertainty;
- [atT, ~, ~, ~,ys, ~, ~, ~, ~, ~, ~, ~, ~, ~,M_,oo_,bayestopt_] = ...
+ [atT, ~, ~, ~,ys, ~, ~, ~, ~, ~, ~, ~, ~, ~,oo_,bayestopt_] = ...
DsgeSmoother(xparam, gend, data, data_index, missing_value, M_, oo_, options_, bayestopt_, estim_params_);
options_.smoothed_state_uncertainty = options_smoothed_state_uncertainty_old;
%get constant part
@@ -171,7 +171,7 @@ if options_.logged_steady_state %if steady state was previously logged, undo thi
options_.logged_steady_state=0;
end
-[T, R, ys, ~, M_, oo_] = dynare_resolve(M_, options_, oo_);
+[T, R, ys, ~, oo_.dr,M_.params] = dynare_resolve(M_, options_, oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if options_.loglinear && isfield(oo_.dr,'ys') && options_.logged_steady_state==0 %log steady state
oo_.dr.ys=log_variable(1:M_.endo_nbr,oo_.dr.ys,M_);
diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
index e966ffaa3b8b137e20c72d7042a9830eb568fe8f..8f738a7349b86b4f6ee42d1eca946b2e21b01b29 100644
--- a/matlab/kalman/likelihood/missing_observations_kalman_filter.m
+++ b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
@@ -113,7 +113,7 @@ if occbin_.status
if isempty(opts_regime.binding_indicator) && isempty(opts_regime.regime_history)
opts_regime.binding_indicator=zeros(last+2,M_.occbin.constraint_nbr);
end
- [~, ~, ~, regimes_] = occbin.check_regimes([], [], [], opts_regime, M_, oo_, options_);
+ [~, ~, ~, regimes_] = occbin.check_regimes([], [], [], opts_regime, M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if length(occbin_.info)>4
TT=occbin_.info{5};
RR=occbin_.info{6};
diff --git a/matlab/minus_logged_prior_density.m b/matlab/minus_logged_prior_density.m
index 19203297a65164034d2b3214070d5b1036432f81..3dd49834c82b7aaaec4eaf672deeedff1d2a038b 100644
--- a/matlab/minus_logged_prior_density.m
+++ b/matlab/minus_logged_prior_density.m
@@ -13,7 +13,7 @@ function [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparam
% f [double] value of minus the logged prior density.
% info [double] vector: second entry stores penalty, first entry the error code.
%
-% Copyright © 2009-2017 Dynare Team
+% Copyright © 2009-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -121,12 +121,11 @@ end
% 2. Check BK and steady state
%-----------------------------
-M_ = set_all_parameters(xparams,EstimatedParams,DynareModel);
-[dr,info,DynareModel,DynareResults] = resol(0,DynareModel,DynareOptions,DynareResults);
+[~,info] = resol(0,DynareModel,DynareOptions,DynareResults.dr,DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
- if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ...
+ if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
info(1) == 81 || info(1) == 84 || info(1) == 85
%meaningful second entry of output that can be used
@@ -142,6 +141,4 @@ if info(1)
end
end
-
-
fval = - priordens(xparams,pshape,p6,p7,p3,p4);
\ No newline at end of file
diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
index d863db5a38507bece3d6c444c65ae89c9a9c635e..1acc1829e2ccd41c6bdb55206d474224d3ea427f 100644
--- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
@@ -174,7 +174,7 @@ else
opts_regime.binding_indicator=zeros(smpl+2,M_.occbin.constraint_nbr);
end
occbin_options.opts_regime = opts_regime;
- [~, ~, ~, regimes_] = occbin.check_regimes([], [], [], opts_regime, M_, oo_, options_);
+ [~, ~, ~, regimes_] = occbin.check_regimes([], [], [], opts_regime, M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if length(occbin_.info)>4
if length(occbin_.info)==6 && options_.smoother_redux
TT=repmat(T,1,1,smpl+1);
@@ -506,7 +506,7 @@ while notsteady && t<smpl
opts_simul.init_regime = []; %regimes_(t);
opts_simul.waitbar=0;
options_.occbin.simul=opts_simul;
- [~, out, ss] = occbin.solver(M_,oo_,options_);
+ [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
for jnk=1:nk
if filter_covariance_flag
diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m
index 58dc3ee293933f095401065cb01124cd9a3cbea2..ecb861fba0c11054bed783676c5e8bd71c2bdb57 100644
--- a/matlab/non_linear_dsge_likelihood.m
+++ b/matlab/non_linear_dsge_likelihood.m
@@ -1,17 +1,17 @@
-function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults] = non_linear_dsge_likelihood(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults)
+function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,oo_] = non_linear_dsge_likelihood(xparam1,DynareDataset,DatasetInfo,options_,M_,EstimatedParameters,bayestopt_,BoundsInfo,oo_)
% Evaluates the posterior kernel of a dsge model using a non linear filter.
%
% INPUTS
% - xparam1 [double] n×1 vector, estimated parameters.
-% - DynareDataset [struct] Matlab's structure containing the dataset (initialized by dynare, aka dataset_).
-% - DatasetInfo [struct] Matlab's structure describing the dataset (initialized by dynare, aka dataset_info).
-% - DynareOptions [struct] Matlab's structure describing the options (initialized by dynare, aka options_).
-% - Model [struct] Matlab's structure describing the Model (initialized by dynare, aka M_).
-% - EstimatedParameters [struct] Matlab's structure describing the estimated_parameters (initialized by dynare, aka estim_params_).
-% - BayesInfo [struct] Matlab's structure describing the priors (initialized by dynare,aka bayesopt_).
-% - BoundsInfo [struct] Matlab's structure specifying the bounds on the paramater values (initialized by dynare,aka bayesopt_).
-% - DynareResults [struct] Matlab's structure gathering the results (initialized by dynare,aka oo_).
+% - DynareDataset [struct] Matlab's structure containing the dataset
+% - DatasetInfo [struct] Matlab's structure describing the dataset
+% - options_ [struct] Matlab's structure describing the options
+% - M_ [struct] Matlab's structure describing the M_
+% - EstimatedParameters [struct] Matlab's structure describing the estimated_parameters
+% - bayestopt_ [struct] Matlab's structure describing the priors
+% - BoundsInfo [struct] Matlab's structure specifying the bounds on the paramater values
+% - oo_ [struct] Matlab's structure gathering the results
%
% OUTPUTS
% - fval [double] scalar, value of the likelihood or posterior kernel.
@@ -21,12 +21,12 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,Model,DynareOptions,Bayes
% - Hess [double] Empty array.
% - ys [double] Empty array.
% - trend_coeff [double] Empty array.
-% - Model [struct] Updated Model structure described in INPUTS section.
-% - DynareOptions [struct] Updated DynareOptions structure described in INPUTS section.
-% - BayesInfo [struct] See INPUTS section.
-% - DynareResults [struct] Updated DynareResults structure described in INPUTS section.
+% - M_ [struct] Updated M_ structure described in INPUTS section.
+% - options_ [struct] Updated options_ structure described in INPUTS section.
+% - bayestopt_ [struct] See INPUTS section.
+% - oo_ [struct] Updated oo_ structure described in INPUTS section.
-% Copyright © 2010-2022 Dynare Team
+% Copyright © 2010-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -60,7 +60,7 @@ if ~isempty(xparam1)
end
% Issue an error if loglinear option is used.
-if DynareOptions.loglinear
+if options_.loglinear
error('non_linear_dsge_likelihood: It is not possible to use a non linear filter with the option loglinear!')
end
@@ -68,9 +68,9 @@ end
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
-Model = set_all_parameters(xparam1,EstimatedParameters,Model);
+M_ = set_all_parameters(xparam1,EstimatedParameters,M_);
-[fval,info,exit_flag,Q,H]=check_bounds_and_definiteness_estimation(xparam1, Model, EstimatedParameters, BoundsInfo);
+[fval,info,exit_flag,Q,H]=check_bounds_and_definiteness_estimation(xparam1, M_, EstimatedParameters, BoundsInfo);
if info(1)
return
end
@@ -79,8 +79,8 @@ end
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
-% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
-[dr, info, Model, DynareResults] = resol(0, Model, DynareOptions, DynareResults);
+% Linearize the model around the deterministic steadystate and extract the matrices of the state equation (T and R).
+[dr, info, M_.params] = resol(0, M_, options_, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 || ...
@@ -100,18 +100,18 @@ if info(1)
end
% Define a vector of indices for the observed variables. Is this really usefull?...
-BayesInfo.mf = BayesInfo.mf1;
+bayestopt_.mf = bayestopt_.mf1;
% Get needed informations for kalman filter routines.
-start = DynareOptions.presample+1;
+start = options_.presample+1;
Y = transpose(DynareDataset.data);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
-mf0 = BayesInfo.mf0;
-mf1 = BayesInfo.mf1;
+mf0 = bayestopt_.mf0;
+mf1 = bayestopt_.mf1;
restrict_variables_idx = dr.restrict_var_list;
state_variables_idx = restrict_variables_idx(mf0);
number_of_state_variables = length(mf0);
@@ -124,10 +124,10 @@ ReducedForm.H = H;
ReducedForm.mf0 = mf0;
ReducedForm.mf1 = mf1;
-if DynareOptions.order>3
+if options_.order>3
ReducedForm.use_k_order_solver = true;
ReducedForm.dr = dr;
- ReducedForm.udr = folded_to_unfolded_dr(dr, Model, DynareOptions);
+ ReducedForm.udr = folded_to_unfolded_dr(dr, M_, options_);
if pruning
error('Pruning is not available for orders > 3');
end
@@ -139,7 +139,7 @@ else
ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:);
- if DynareOptions.order==3
+ if options_.order==3
ReducedForm.ghxxx = dr.ghxxx(restrict_variables_idx,:);
ReducedForm.ghuuu = dr.ghuuu(restrict_variables_idx,:);
ReducedForm.ghxxu = dr.ghxxu(restrict_variables_idx,:);
@@ -150,22 +150,22 @@ else
end
% Set initial condition.
-switch DynareOptions.particle.initialization
+switch options_.particle.initialization
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
StateVectorMean = ReducedForm.constant(mf0);
[A,B] = kalman_transition_matrix(dr,dr.restrict_var_list,dr.restrict_columns);
- StateVectorVariance = lyapunov_symm(A, B*Q*B', DynareOptions.lyapunov_fixed_point_tol, ...
- DynareOptions.qz_criterium, DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
+ StateVectorVariance = lyapunov_symm(A, B*Q*B', options_.lyapunov_fixed_point_tol, ...
+ options_.qz_criterium, options_.lyapunov_complex_threshold, [], options_.debug);
StateVectorVariance = StateVectorVariance(mf0,mf0);
case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
StateVectorMean = ReducedForm.constant(mf0);
- old_DynareOptionsperiods = DynareOptions.periods;
- DynareOptions.periods = 5000;
- old_DynareOptionspruning = DynareOptions.pruning;
- DynareOptions.pruning = DynareOptions.particle.pruning;
- y_ = simult(DynareResults.steady_state, dr,Model,DynareOptions,DynareResults);
+ old_DynareOptionsperiods = options_.periods;
+ options_.periods = 5000;
+ old_DynareOptionspruning = options_.pruning;
+ options_.pruning = options_.particle.pruning;
+ y_ = simult(oo_.steady_state, dr,M_,options_);
y_ = y_(dr.order_var(state_variables_idx),2001:5000); %state_variables_idx is in dr-order while simult_ is in declaration order
- if any(any(isnan(y_))) || any(any(isinf(y_))) && ~ DynareOptions.pruning
+ if any(any(isnan(y_))) || any(any(isinf(y_))) && ~ options_.pruning
fval = Inf;
info(1) = 202;
info(4) = 0.1;
@@ -173,13 +173,13 @@ switch DynareOptions.particle.initialization
return;
end
StateVectorVariance = cov(y_');
- DynareOptions.periods = old_DynareOptionsperiods;
- DynareOptions.pruning = old_DynareOptionspruning;
+ options_.periods = old_DynareOptionsperiods;
+ options_.pruning = old_DynareOptionspruning;
clear('old_DynareOptionsperiods','y_');
case 3% Initial state vector covariance is a diagonal matrix (to be used
% if model has stochastic trends).
StateVectorMean = ReducedForm.constant(mf0);
- StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables);
+ StateVectorVariance = options_.particle.initial_state_prior_std*eye(number_of_state_variables);
otherwise
error('Unknown initialization option!')
end
@@ -197,9 +197,9 @@ end
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
-DynareOptions.warning_for_steadystate = 0;
+options_.warning_for_steadystate = 0;
[s1,s2] = get_dynare_random_generator_state();
-LIK = feval(DynareOptions.particle.algorithm, ReducedForm, Y, start, DynareOptions.particle, DynareOptions.threads, DynareOptions, Model);
+LIK = feval(options_.particle.algorithm, ReducedForm, Y, start, options_.particle, options_.threads, options_, M_);
set_dynare_random_generator_state(s1,s2);
if imag(LIK)
fval = Inf;
@@ -216,11 +216,11 @@ elseif isnan(LIK)
else
likelihood = LIK;
end
-DynareOptions.warning_for_steadystate = 1;
+options_.warning_for_steadystate = 1;
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
-lnprior = priordens(xparam1(:),BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
+lnprior = priordens(xparam1(:),bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
fval = (likelihood-lnprior);
if isnan(fval)
diff --git a/matlab/nonlinear-filters/solve_model_for_online_filter.m b/matlab/nonlinear-filters/solve_model_for_online_filter.m
index fb081b359d309dfcf3b9c356c6fd275728e7eded..a0b607bee3f6a1b037154bb26a73553ed5011cf4 100644
--- a/matlab/nonlinear-filters/solve_model_for_online_filter.m
+++ b/matlab/nonlinear-filters/solve_model_for_online_filter.m
@@ -125,8 +125,8 @@ Model.H = H;
%------------------------------------------------------------------------------
warning('off', 'MATLAB:nearlySingularMatrix')
-[~, ~, ~, info, Model, DynareResults] = ...
- dynare_resolve(Model, DynareOptions, DynareResults, 'restrict');
+[~, ~, ~, info, DynareResults.dr, Model.params] = ...
+ dynare_resolve(Model, DynareOptions, DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state, 'restrict');
warning('on', 'MATLAB:nearlySingularMatrix')
if info(1)~=0
@@ -197,7 +197,7 @@ if setinitialcondition
DynareOptions.periods = 5000;
old_DynareOptionspruning = DynareOptions.pruning;
DynareOptions.pruning = DynareOptions.particle.pruning;
- y_ = simult(dr.ys, dr, Model, DynareOptions, DynareResults);
+ y_ = simult(dr.ys, dr, Model, DynareOptions);
y_ = y_(dr.order_var(state_variables_idx),2001:DynareOptions.periods);
StateVectorVariance = cov(y_');
DynareOptions.periods = old_DynareOptionsperiods;
diff --git a/matlab/optimization/gmhmaxlik_core.m b/matlab/optimization/gmhmaxlik_core.m
index c75a946eb78e4a90d07abebd84f086798182ff36..9f404f59b69ce4e3cf19662e299a4db806d4ce13 100644
--- a/matlab/optimization/gmhmaxlik_core.m
+++ b/matlab/optimization/gmhmaxlik_core.m
@@ -73,9 +73,6 @@ function [PostMod,PostVar,Scale,PostMean] = gmhmaxlik_core(ObjFun,xparam1,mh_bou
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global bayestopt_ estim_params_ options_
-
-options_.lik_algo = 1;
npar = length(xparam1);
NumberOfIterations = options.number;
diff --git a/matlab/optimize_prior.m b/matlab/optimize_prior.m
index cfad0c0907354830a4e90e6339e4da8d899d0fb0..748495701c40fbfdbe5fa9ae44b90cbd6d5228b1 100644
--- a/matlab/optimize_prior.m
+++ b/matlab/optimize_prior.m
@@ -29,7 +29,7 @@ while look_for_admissible_initial_condition
xinit = xparam1+scale*randn(size(xparam1));
if all(xinit(:)>BayesInfo.p3) && all(xinit(:)<BayesInfo.p4)
ModelInfo = set_all_parameters(xinit,EstimationInfo,ModelInfo);
- [dr,INFO,ModelInfo,DynareResults] = resol(0,ModelInfo,DynareOptions,DynareResults);
+ [DynareResults.dr,INFO,ModelInfo.params] = resol(0,ModelInfo,DynareOptions,DynareResults.dr,DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
if ~INFO(1)
look_for_admissible_initial_condition = 0;
end
diff --git a/matlab/perfect-foresight-models/det_cond_forecast.m b/matlab/perfect-foresight-models/det_cond_forecast.m
index c6466bea636722ff4b6996d60fd35e9d642ce613..dd73cfb3efe7fe5e213ca7ded295e2421b577baf 100644
--- a/matlab/perfect-foresight-models/det_cond_forecast.m
+++ b/matlab/perfect-foresight-models/det_cond_forecast.m
@@ -47,7 +47,7 @@ if ~isfield(oo_,'dr') || ~isfield(oo_.dr,'ghx')
dr = struct();
oo_.dr=set_state_space(dr,M_,options_);
options_.order = 1;
- [dr,Info,M_,oo_] = resol(0,M_,options_,oo_);
+ [oo_.dr,Info,M_.params] = resol(0,M_,options_,oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
fprintf('done\n');
end
b_surprise = 0;
diff --git a/matlab/pm3.m b/matlab/pm3.m
index 3800cd51565f68f7cee9162332b06f7b84051adf..d08b068304bc008e0b55d1867adc2eccf0e77faa 100644
--- a/matlab/pm3.m
+++ b/matlab/pm3.m
@@ -1,4 +1,4 @@
-function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryName,var_type,dispString)
+function oo_=pm3(M_,options_,oo_,n1,n2,ifil,B,tit1,tit2,tit_tex,names1,names2,name3,DirectoryName,var_type,dispString)
% Computes, stores and plots the posterior moment statistics.
%
@@ -8,8 +8,7 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa
% ifil [scalar] number of moment files to load
% B [scalar] number of subdraws
% tit1 [string] Figure title
-% tit2 [string] not used
-% tit3 [string] Save name for figure
+% tit2 [string] Save name for figure
% tit_tex [cell array] TeX-Names for Variables
% names1 [cell array] Names of all variables in the moment matrix from
% which names2 is selected
@@ -20,6 +19,10 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa
% var_type [string] suffix of the filename from which to load moment
% matrix
% dispString [string] string to be displayes in the command window
+%
+% OUTPUTS
+% oo_ [structure] storing the results
+
% PARALLEL CONTEXT
% See also the comment in posterior_sampler.m funtion.
@@ -42,8 +45,6 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global options_ M_ oo_
-
nn = 3;
MaxNumberOfPlotsPerFigure = nn^2; % must be square
varlist = names2;
@@ -284,9 +285,7 @@ if strcmp(var_type,'_trend_coeff') || max(max(abs(Mean(:,:))))<=10^(-6) || all(a
fprintf(['%s: ' tit1 ', done!\n'],dispString);
return %not do plots
end
-%%
%% Finally I build the plots.
-%%
if ~options_.nograph && ~options_.no_graph.posterior
% Block of code executed in parallel, with the exception of file
@@ -297,8 +296,6 @@ if ~options_.nograph && ~options_.no_graph.posterior
%
% %%% The file .TeX! are not saved in parallel.
-
-
% Store the variable mandatory for local/remote parallel computing.
localVars=[];
@@ -313,8 +310,14 @@ if ~options_.nograph && ~options_.no_graph.posterior
end
localVars.MaxNumberOfPlotsPerFigure=MaxNumberOfPlotsPerFigure;
localVars.name3=name3;
- localVars.tit3=tit3;
+ localVars.tit2=tit2;
localVars.Mean=Mean;
+ localVars.TeX=options_.TeX;
+ localVars.nodisplay=options_.nodisplay;
+ localVars.graph_format=options_.graph_format;
+ localVars.dname=M_.dname;
+ localVars.fname=M_.fname;
+
% Like sequential execution!
nvar0=nvar;
@@ -332,16 +335,13 @@ if ~options_.nograph && ~options_.no_graph.posterior
if isRemoteOctave
fout = pm3_core(localVars,1,nvar,0);
else
- globalVars = struct('M_',M_, ...
- 'options_', options_, ...
- 'oo_', oo_);
+ globalVars = [];
[fout, nvar0, totCPU] = masterParallel(options_.parallel, 1, nvar, [],'pm3_core', localVars,globalVars, options_.parallel_info);
end
end
else
% For the time being in Octave enviroment the pm3.m is executed only in
% serial modality, to avoid problem with the plots.
-
fout = pm3_core(localVars,1,nvar,0);
end
@@ -365,8 +365,8 @@ if ~options_.nograph && ~options_.no_graph.posterior
if subplotnum == MaxNumberOfPlotsPerFigure || i == nvar
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
- fprintf(fidTeX,['\\includegraphics[width=%2.2f\\textwidth]{%s/Output/%s_' name3 '_%s}\n'],options_.figures.textwidth*min(subplotnum/nn,1),M_.dname,M_.fname, tit3{i});
- fprintf(fidTeX,'\\label{Fig:%s:%s}\n',name3,tit3{i});
+ fprintf(fidTeX,['\\includegraphics[width=%2.2f\\textwidth]{%s/Output/%s_' name3 '_%s}\n'],options_.figures.textwidth*min(subplotnum/nn,1),M_.dname,M_.fname, tit2{i});
+ fprintf(fidTeX,'\\label{Fig:%s:%s}\n',name3,tit2{i});
fprintf(fidTeX,'\\caption{%s}\n',tit1);
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
diff --git a/matlab/pm3_core.m b/matlab/pm3_core.m
index 343c708975b35b25242a3d90d407e0ece740f638..f820ed2337fc20a475fa46f143a6c428accd5109 100644
--- a/matlab/pm3_core.m
+++ b/matlab/pm3_core.m
@@ -1,13 +1,23 @@
function myoutput=pm3_core(myinputs,fpar,nvar,whoiam, ThisMatlab)
-
+% myoutput=pm3_core(myinputs,fpar,nvar,whoiam, ThisMatlab)
% PARALLEL CONTEXT
% Core functionality for pm3.m function, which can be parallelized.
% INPUTS
-% See the comment in posterior_sampler_core.m funtion.
-
+% o myimput [struc] The mandatory variables for local/remote
+% parallel computing obtained from prior_posterior_statistics.m
+% function.
+% o fpar and nvar [integer] first variable and number of variables
+% o whoiam [integer] In concurrent programming a modality to refer to the different threads running in parallel is needed.
+% The integer whoaim is the integer that
+% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
+% cluster.
+% o ThisMatlab [integer] Allows us to distinguish between the
+% 'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab,
+% ... Then it is the index number of this slave machine in the cluster.
+%
% OUTPUTS
-% o myoutput [struc]
+% o myoutput [struct] Contains file names
%
%
% SPECIAL REQUIREMENTS.
@@ -45,16 +55,18 @@ varlist=myinputs.varlist;
MaxNumberOfPlotsPerFigure=myinputs.MaxNumberOfPlotsPerFigure;
name3=myinputs.name3;
-tit3=myinputs.tit3;
+tit2=myinputs.tit2;
Mean=myinputs.Mean;
+options_.TeX=myinputs.TeX;
+options_.nodisplay=myinputs.nodisplay;
+options_.graph_format=myinputs.graph_format;
+fname=myinputs.fname;
+dname=myinputs.dname;
if whoiam
Parallel=myinputs.Parallel;
end
-
-global options_ M_ oo_
-
if options_.TeX
varlist_TeX=myinputs.varlist_TeX;
end
@@ -64,8 +76,6 @@ if whoiam
h = dyn_waitbar(prct0,'Parallel plots pm3 ...');
end
-
-
figunumber = 0;
subplotnum = 0;
hh_fig = dyn_figure(options_.nodisplay,'Name',[tit1 ' ' int2str(figunumber+1)]);
@@ -110,14 +120,14 @@ for i=fpar:nvar
if whoiam
if Parallel(ThisMatlab).Local==0
- DirectoryName = CheckPath('Output',M_.dname);
+ DirectoryName = CheckPath('Output',dname);
end
end
if subplotnum == MaxNumberOfPlotsPerFigure || i == nvar
- dyn_saveas(hh_fig,[M_.dname '/Output/' M_.fname '_' name3 '_' tit3{i}],options_.nodisplay,options_.graph_format);
+ dyn_saveas(hh_fig,[dname '/Output/' fname '_' name3 '_' tit2{i}],options_.nodisplay,options_.graph_format);
if RemoteFlag==1
- OutputFileName = [OutputFileName; {[M_.dname, filesep, 'Output',filesep], [M_.fname '_' name3 '_' deblank(tit3(i,:)) '.*']}];
+ OutputFileName = [OutputFileName; {[dname, filesep, 'Output',filesep], [fname '_' name3 '_' deblank(tit2(i,:)) '.*']}];
end
subplotnum = 0;
figunumber = figunumber+1;
@@ -127,12 +137,8 @@ for i=fpar:nvar
end
if whoiam
- % waitbarString = [ 'Variable ' int2str(i) '/' int2str(nvar) ' done.'];
- % fMessageStatus((i-fpar+1)/(nvar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
dyn_waitbar((i-fpar+1)/(nvar-fpar+1),h);
end
-
-
end
if whoiam
diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m
index 2c8afc8468afbc0a5b33b2359d791af265e722d1..3b523c0a7c392bda8dbfacea7d3cb63c13bd9ff9 100644
--- a/matlab/posterior_analysis.m
+++ b/matlab/posterior_analysis.m
@@ -1,4 +1,4 @@
-function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_)
+function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_,estim_params_)
% Copyright © 2008-2021 Dynare Team
%
@@ -29,7 +29,7 @@ switch info
if drsize*SampleSize>MaxMegaBytes
drsize=0;
end
- SampleAddress = selec_posterior_draws(SampleSize,drsize);
+ selec_posterior_draws(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state,estim_params_,SampleSize,drsize); %save draws to disk
oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
case {4,5}
oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m
index 4e3379a6443a26a1c9bf14240b75741ebec77d41..9b751c558bd0eda099f826962f2ae249c3d305c0 100644
--- a/matlab/posterior_sampler_core.m
+++ b/matlab/posterior_sampler_core.m
@@ -93,7 +93,6 @@ ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName];
save_tmp_file = sampler_options.save_tmp_file;
-options_.lik_algo = 1;
OpenOldFile = ones(nblck,1);
if strcmpi(ProposalFun,'rand_multivariate_normal')
sampler_options.n = npar;
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index 4bdddfd2d8182a43f4dab2330b7d26ca0e5edf9e..58f428c5a192a2e1ca26343b92f4f4a01123fa12 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -1,17 +1,19 @@
-function prior_posterior_statistics(type,dataset,dataset_info,dispString)
-% function prior_posterior_statistics(type,dataset,dataset_info,dispString)
+function oo_=prior_posterior_statistics(type,dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_,dispString)
+% oo_=prior_posterior_statistics(type,dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_,dispString))
% Computes Monte Carlo filter smoother and forecasts
%
% INPUTS
-% type: posterior
-% prior
-% gsa
-% dataset: data structure
-% dataset_info: dataset structure
-% dispString: string to display in the command window
-%
+% type [string] posterior, prior, or gsa
+% o dataset_ [structure] storing the dataset
+% o dataset_info [structure] Various information about the dataset
+% o M_ [structure] storing the model information
+% o oo_ [structure] storing the results
+% o options_ [structure] storing the options
+% o estim_params_ [structure] storing information about estimated parameters
+% o bayestopt_ [structure] storing information about priors
+% dispString: [string] display info in the command window
% OUTPUTS
-% none
+% oo_: [structure] storing the results
%
% SPECIAL REQUIREMENTS
% none
@@ -37,36 +39,23 @@ function prior_posterior_statistics(type,dataset,dataset_info,dispString)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-if nargin < 4
+if nargin < 9
dispString = 'prior_posterior_statistics';
end
-
-global options_ estim_params_ oo_ M_ bayestopt_
-
localVars=[];
-Y = transpose(dataset.data);
-gend = dataset.nobs;
-data_index = dataset_info.missing.aindex;
-missing_value = dataset_info.missing.state;
-mean_varobs = dataset_info.descriptive.mean;
-
+Y = transpose(dataset_.data);
+gend = dataset_.nobs;
-nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
-ncx = estim_params_.ncx;
-ncn = estim_params_.ncn;
-np = estim_params_.np ;
-npar = nvx+nvn+ncx+ncn+np;
+npar = estim_params_.nvx+nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np;
naK = length(options_.filter_step_ahead);
MaxNumberOfBytes=options_.MaxNumberOfBytes;
endo_nbr=M_.endo_nbr;
exo_nbr=M_.exo_nbr;
meas_err_nbr=length(M_.Correlation_matrix_ME);
-iendo = 1:endo_nbr;
horizon = options_.forecast;
-IdObs = bayestopt_.mfys;
if horizon
i_last_obs = gend+(1-M_.maximum_endo_lag:0);
end
@@ -130,13 +119,6 @@ varlist = options_.varlist;
if isempty(varlist)
varlist = sort(M_.endo_names(1:M_.orig_endo_nbr));
end
-nvar = length(varlist);
-SelecVariables = [];
-for i=1:nvar
- if ~isempty(strmatch(varlist{i}, M_.endo_names, 'exact'))
- SelecVariables = [SelecVariables; strmatch(varlist{i}, M_.endo_names, 'exact')];
- end
-end
n_variables_to_fill=13;
@@ -171,17 +153,17 @@ localVars.filter_covariance=filter_covariance;
localVars.smoothed_state_uncertainty=smoothed_state_uncertainty;
localVars.gend=gend;
localVars.Y=Y;
-localVars.data_index=data_index;
-localVars.missing_value=missing_value;
+localVars.data_index=dataset_info.missing.aindex;
+localVars.missing_value=dataset_info.missing.state;
localVars.varobs=options_.varobs;
-localVars.mean_varobs=mean_varobs;
+localVars.mean_varobs=dataset_info.descriptive.mean;
localVars.irun=irun;
localVars.endo_nbr=endo_nbr;
localVars.nvn=nvn;
localVars.naK=naK;
localVars.horizon=horizon;
-localVars.iendo=iendo;
-localVars.IdObs=IdObs;
+localVars.iendo=1:endo_nbr;
+localVars.IdObs=bayestopt_.mfys;
if horizon
localVars.i_last_obs=i_last_obs;
localVars.MAX_nforc1=MAX_nforc1;
@@ -212,6 +194,13 @@ localVars.MAX_momentsno = MAX_momentsno;
localVars.ifil=ifil;
localVars.DirectoryName = DirectoryName;
+localVars.M_=M_;
+localVars.oo_=oo_;
+localVars.options_=options_;
+localVars.estim_params_=estim_params_;
+localVars.bayestopt_=bayestopt_;
+
+
if strcmpi(type,'posterior')
record=load_last_mh_history_file(DirectoryName, M_.fname);
[nblck, npar] = size(record.LastParameters);
@@ -290,11 +279,7 @@ else
end
end
localVars.ifil = ifil;
- globalVars = struct('M_',M_, ...
- 'options_', options_, ...
- 'bayestopt_', bayestopt_, ...
- 'estim_params_', estim_params_, ...
- 'oo_', oo_);
+ globalVars = [];
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '.static.m']};
NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']};
@@ -330,28 +315,28 @@ if ~isnumeric(options_.parallel)
end
if options_.smoother
- pm3(endo_nbr,gend,ifil(1),B,'Smoothed variables',...
- '',varlist, M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(1),B,'Smoothed variables',...
+ varlist, M_.endo_names_tex,M_.endo_names,...
varlist,'SmoothedVariables',DirectoryName,'_smooth',dispString);
- pm3(exo_nbr,gend,ifil(2),B,'Smoothed shocks',...
- '',M_.exo_names,M_.exo_names_tex,M_.exo_names,...
+ oo_=pm3(M_,options_,oo_,exo_nbr,gend,ifil(2),B,'Smoothed shocks',...
+ M_.exo_names,M_.exo_names_tex,M_.exo_names,...
M_.exo_names,'SmoothedShocks',DirectoryName,'_inno',dispString);
- pm3(endo_nbr,1,ifil(9),B,'Trend_coefficients',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,1,ifil(9),B,'Trend_coefficients',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'TrendCoeff',DirectoryName,'_trend_coeff',dispString);
- pm3(endo_nbr,gend,ifil(10),B,'Smoothed constant',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(10),B,'Smoothed constant',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'Constant',DirectoryName,'_smoothed_constant',dispString);
- pm3(endo_nbr,gend,ifil(11),B,'Smoothed trend',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(11),B,'Smoothed trend',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'Trend',DirectoryName,'_smoothed_trend',dispString);
- pm3(endo_nbr,gend,ifil(1),B,'Updated Variables',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(1),B,'Updated Variables',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'UpdatedVariables',DirectoryName, ...
'_update',dispString);
if smoothed_state_uncertainty
- pm3(endo_nbr,endo_nbr,ifil(13),B,'State Uncertainty',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,endo_nbr,ifil(13),B,'State Uncertainty',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'StateUncertainty',DirectoryName,'_state_uncert',dispString);
end
@@ -360,24 +345,24 @@ if options_.smoother
meas_error_names{obs_iter,1}=['SE_EOBS_' M_.endo_names{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}];
texnames{obs_iter,1}=['\sigma^{ME}_' M_.endo_names_tex{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}];
end
- pm3(meas_err_nbr,gend,ifil(3),B,'Smoothed measurement errors',...
- '',meas_error_names,texnames,meas_error_names,...
- meas_error_names,'SmoothedMeasurementErrors',DirectoryName,'_error',dispString)
+ oo_=pm3(M_,options_,oo_,meas_err_nbr,gend,ifil(3),B,'Smoothed measurement errors',...
+ meas_error_names,texnames,meas_error_names,...
+ meas_error_names,'SmoothedMeasurementErrors',DirectoryName,'_error',dispString);
end
end
if options_.filtered_vars
- pm3(endo_nbr,gend,ifil(4),B,'One step ahead forecast (filtered variables)',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(4),B,'One step ahead forecast (filtered variables)',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead',dispString);
end
if options_.forecast
- pm3(endo_nbr,horizon,ifil(6),B,'Forecasted variables (mean)',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,horizon,ifil(6),B,'Forecasted variables (mean)',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'MeanForecast',DirectoryName,'_forc_mean',dispString);
- pm3(endo_nbr,horizon,ifil(7),B,'Forecasted variables (point)',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,horizon,ifil(7),B,'Forecasted variables (point)',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'PointForecast',DirectoryName,'_forc_point',dispString);
if ~isequal(M_.H,0) && ~isempty(intersect(options_.varobs,varlist))
texnames = cell(length(options_.varobs), 1);
@@ -387,15 +372,15 @@ if options_.forecast
texnames{obs_iter}=M_.endo_names_tex{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')};
end
varlist_forecast_ME=intersect(options_.varobs,varlist);
- pm3(meas_err_nbr,horizon,ifil(12),B,'Forecasted variables (point) with ME',...
- '',varlist_forecast_ME,texnames,obs_names,...
- varlist_forecast_ME,'PointForecastME',DirectoryName,'_forc_point_ME',dispString)
+ oo_=pm3(M_,options_,oo_,meas_err_nbr,horizon,ifil(12),B,'Forecasted variables (point) with ME',...
+ varlist_forecast_ME,texnames,obs_names,...
+ varlist_forecast_ME,'PointForecastME',DirectoryName,'_forc_point_ME',dispString);
end
end
if options_.filter_covariance
- pm3(endo_nbr,endo_nbr,ifil(8),B,'Filtered covariances',...
- '',varlist,M_.endo_names_tex,M_.endo_names,...
+ oo_=pm3(M_,options_,oo_,endo_nbr,endo_nbr,ifil(8),B,'Filtered covariances',...
+ varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'FilterCovariance',DirectoryName,'_filter_covar',dispString);
end
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index babb70c3b3e4374d9f7baef7360e3ae929c623b2..273cf0367b5490664cfa4cd4015ff7dd14434199 100644
--- a/matlab/prior_posterior_statistics_core.m
+++ b/matlab/prior_posterior_statistics_core.m
@@ -47,14 +47,17 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global options_ oo_ M_ bayestopt_ estim_params_
-
if nargin<4
whoiam=0;
end
% Reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
+M_=myinputs.M_;
+oo_=myinputs.oo_;
+options_=myinputs.options_;
+estim_params_=myinputs.estim_params_;
+bayestopt_=myinputs.bayestopt_;
type=myinputs.type;
run_smoother=myinputs.run_smoother;
@@ -91,7 +94,6 @@ if smoothed_state_uncertainty
MAX_n_smoothed_state_uncertainty=myinputs.MAX_n_smoothed_state_uncertainty;
end
-exo_nbr=myinputs.exo_nbr;
maxlag=myinputs.maxlag;
MAX_nsmoo=myinputs.MAX_nsmoo;
MAX_ninno=myinputs.MAX_ninno;
@@ -100,7 +102,6 @@ MAX_n_smoothed_trend=myinputs.MAX_n_smoothed_trend;
MAX_n_trend_coeff=myinputs.MAX_n_trend_coeff;
MAX_nerro = myinputs.MAX_nerro;
MAX_nruns=myinputs.MAX_nruns;
-MAX_momentsno = myinputs.MAX_momentsno;
ifil=myinputs.ifil;
if ~strcmpi(type,'prior')
@@ -214,10 +215,10 @@ for b=fpar:B
M_ = set_all_parameters(deep,estim_params_,M_);
if run_smoother
- [dr,info,M_,oo_] =compute_decision_rules(M_,opts_local,oo_);
+ [dr,info,M_.params] =compute_decision_rules(M_,opts_local,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if ismember(info(1),[3,4])
opts_local.qz_criterium = 1 + (opts_local.qz_criterium-1)*10; %loosen tolerance, useful for big models where BK conditions were tested on restricted state space
- [dr,info,M_,oo_] =compute_decision_rules(M_,opts_local,oo_);
+ [dr,info,M_.params] =compute_decision_rules(M_,opts_local,oo_);
end
if info(1)
message=get_error_message(info,opts_local);
@@ -228,14 +229,14 @@ for b=fpar:B
opts_local.occbin.simul.waitbar=0;
opts_local.occbin.smoother.waitbar = 0;
opts_local.occbin.smoother.linear_smoother=false; % speed-up
- [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,M_,oo_,bayestopt_] = ...
+ [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_] = ...
occbin.DSGE_smoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);
if oo_.occbin.smoother.error_flag(1)
message=get_error_message(oo_.occbin.smoother.error_flag,opts_local);
fprintf('\nprior_posterior_statistics: One of the draws failed with the error:\n%s\n',message)
end
else
- [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,M_,oo_,bayestopt_] = ...
+ [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_] = ...
DsgeSmoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);
end
@@ -311,7 +312,7 @@ for b=fpar:B
end
if horizon
yyyy = alphahat(iendo,i_last_obs);
- yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
+ yf = simulate_posterior_forecasts(yyyy,dr,horizon,false,M_.Sigma_e,1);
if options_.prefilter
% add mean
yf(:,IdObs) = yf(:,IdObs)+repmat(mean_varobs, ...
@@ -328,7 +329,7 @@ for b=fpar:B
else
yf = yf+repmat(SteadyState',horizon+maxlag,1);
end
- yf1 = forcst2(yyyy,horizon,dr,1);
+ yf1 = simulate_posterior_forecasts(yyyy,dr,horizon,false,M_.Sigma_e,1);
if options_.prefilter == 1
% add mean
yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
@@ -367,7 +368,7 @@ for b=fpar:B
stock_smoothed_uncert(dr.order_var,dr.order_var,:,irun(13)) = state_uncertainty;
end
else
- [T,R,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_);
+ [~,~,SteadyState,info] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
stock_param(irun(5),:) = deep;
stock_logpo(irun(5),1) = logpo;
@@ -540,3 +541,56 @@ if RemoteFlag==1
end
dyn_waitbar_close(h);
+
+
+function yf=simulate_posterior_forecasts(y0,dr,horizon,stochastic_indicator,Sigma_e,n)
+% function yf=forcst2(y0,horizon,dr,n)
+%
+% computes forecasts based on first order model solution, given shocks
+% drawn from the shock distribution, but not including measurement error
+% Inputs:
+% - y0 [endo_nbr by maximum_endo_lag] matrix of starting values
+% - dr [structure] structure with Dynare decision rules
+% - horizon [scalar] number of forecast periods
+% - stochastic_indicator [boolean] indicator whether to consider stochastic shocks
+% - Sigma_e [integer] covariance matrix of shocks
+% - n [scalar] number of repetitions
+%
+% Outputs:
+% - yf [horizon+ykmin_ by endo_nbr by n] array of forecasts
+
+if nargin< 4
+ stochastic_indicator=false;
+ n=1;
+end
+%select states
+k2 = dr.inv_order_var(dr.state_var);
+
+if stochastic_indicator
+ % eliminate shocks with 0 variance
+ i_exo_var = setdiff(1:length(Sigma_e),find(diag(Sigma_e) == 0));
+ nxs = length(i_exo_var);
+
+ chol_S = chol(Sigma_e(i_exo_var,i_exo_var));
+
+ if ~isempty(Sigma_e)
+ e = randn(nxs,n,horizon);
+ end
+
+ B1 = dr.ghu(:,i_exo_var)*chol_S';
+end
+endo_nbr=length(y0);
+
+yf = zeros(endo_nbr,1+horizon,n);
+yf(:,1,:,:) = repmat(y0,[1,1,n]);
+
+for iter=1:horizon
+ if stochastic_indicator
+ yf(:,iter+1,:) = dr.ghx*squeeze(yf(k2,iter,:))+B1*squeeze(e(:,:,iter));
+ else
+ yf(:,iter+1,:) = dr.ghx*squeeze(yf(k2,iter,:));
+ end
+end
+
+yf(dr.order_var,:,:) = yf;
+yf=permute(yf,[2 1 3]);
\ No newline at end of file
diff --git a/matlab/prior_sampler.m b/matlab/prior_sampler.m
index 3ef74fd45bda8292f4059a643c2ebaa8fcb4ce37..4c8b4cbcc8d5c30502280625a92c9e795e51b799 100644
--- a/matlab/prior_sampler.m
+++ b/matlab/prior_sampler.m
@@ -96,8 +96,7 @@ while iteration < NumberOfSimulations
loop_indx = loop_indx+1;
params = Prior.draw();
M_ = set_all_parameters(params, estim_params_, M_);
- [T, R, ~, INFO, M_, oo_] = dynare_resolve(M_, options_, oo_, 'restrict');
- dr = oo_.dr;
+ [T, R, ~, INFO, oo_.dr,M_.params] = dynare_resolve(M_, options_, oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state, 'restrict');
if ~INFO(1)
INFO=endogenous_prior_restrictions(T,R,M_,options_,oo_);
end
@@ -113,7 +112,7 @@ while iteration < NumberOfSimulations
switch INFO(1)
case 0
if drsave
- pdraws(file_line_number,3) = {dr};
+ pdraws(file_line_number,3) = {oo_.dr};
end
[sampled_prior_expectation,sampled_prior_covariance] = ...
recursive_prior_moments(sampled_prior_expectation,sampled_prior_covariance,params,iteration);
diff --git a/matlab/realtime_shock_decomposition.m b/matlab/realtime_shock_decomposition.m
index 728d83230e1b1c9504471a710da816b844239a63..33eaa97e3a52dde0d73a1e71cb0afba98ffb332d 100644
--- a/matlab/realtime_shock_decomposition.m
+++ b/matlab/realtime_shock_decomposition.m
@@ -116,7 +116,7 @@ nobs = options_.nobs;
if forecast_ && any(forecast_params)
M1=M_;
M1.params = forecast_params;
- [~,~,~,~,~,oo1] = dynare_resolve(M1,options_,oo_);
+ [~,~,~,~,~,dr1] = dynare_resolve(M1,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
gend0=0;
@@ -183,8 +183,8 @@ for j=presample+1:nobs
if forecast_
if any(forecast_params)
- Af = oo1.dr.ghx;
- Bf = oo1.dr.ghu;
+ Af = dr1.ghx;
+ Bf = dr1.ghu;
else
Af = A;
Bf = B;
diff --git a/matlab/resol.m b/matlab/resol.m
index a99d94a26ac570c54f6e79e193d49b4f56560938..85cb9b9d1f237f89acb7d9fff514a4fc5c884d5e 100644
--- a/matlab/resol.m
+++ b/matlab/resol.m
@@ -1,18 +1,20 @@
-function [dr, info, M_, oo_] = resol(check_flag, M_, options_, oo_)
-%[dr, info, M_, oo_] = resol(check_flag, M_, options_, oo_)
+function [dr, info, params] = resol(check_flag, M_, options_, dr_in, endo_steady_state, exo_steady_state, exo_det_steady_state)
+%[dr, info, params] = resol(check_flag, M_, options_, dr_in, endo_steady_state, exo_steady_state, exo_det_steady_state)
% Computes the perturbation based decision rules of the DSGE model (orders 1 to 3)
%
% INPUTS
% - check_flag [integer] scalar, equal to 0 if all the approximation is required, equal to 1 if only the eigenvalues are to be computed.
% - M_ [structure] Matlab's structure describing the model
% - options_ [structure] Matlab's structure describing the current options
-% - oo_ [structure] Matlab's structure containing the results
+% - dr_in [structure] model information structure
+% - endo_steady_state [vector] steady state value for endogenous variables
+% - exo_steady_state [vector] steady state value for exogenous variables
+% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - dr [structure] Reduced form model.
% - info [integer] scalar or vector, error code.
-% - M_ [structure] Matlab's structure describing the model
-% - oo_ [structure] Matlab's structure containing the results
+% - params [double] vector of potentially updated parameters
%
% REMARKS
% Possible values for the error codes are:
@@ -32,7 +34,7 @@ function [dr, info, M_, oo_] = resol(check_flag, M_, options_, oo_)
% info(1)=24 -> M_.params has been updated in the steadystate routine and has some NaNs.
% info(1)=30 -> Ergodic variance can't be computed.
-% Copyright © 2001-2023 Dynare Team
+% Copyright © 20012023 Dynare Team
%
% This file is part of Dynare.
%
@@ -49,35 +51,34 @@ function [dr, info, M_, oo_] = resol(check_flag, M_, options_, oo_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-if isfield(oo_,'dr')
- if isfield(oo_.dr,'kstate')
- dr.kstate = oo_.dr.kstate;
- end
- if isfield(oo_.dr,'inv_order_var')
- dr.inv_order_var = oo_.dr.inv_order_var;
- end
- if isfield(oo_.dr,'order_var')
- dr.order_var = oo_.dr.order_var;
- end
- if isfield(oo_.dr,'restrict_var_list')
- dr.restrict_var_list = oo_.dr.restrict_var_list;
- end
- if isfield(oo_.dr,'restrict_columns')
- dr.restrict_columns = oo_.dr.restrict_columns;
- end
- if isfield(oo_.dr,'obs_var')
- dr.obs_var = oo_.dr.obs_var;
- end
+%preserve only the following fields:
+if isfield(dr_in,'kstate')
+ dr.kstate = dr_in.kstate;
+end
+if isfield(dr_in,'inv_order_var')
+ dr.inv_order_var = dr_in.inv_order_var;
+end
+if isfield(dr_in,'order_var')
+ dr.order_var = dr_in.order_var;
+end
+if isfield(dr_in,'restrict_var_list')
+ dr.restrict_var_list = dr_in.restrict_var_list;
+end
+if isfield(dr_in,'restrict_columns')
+ dr.restrict_columns = dr_in.restrict_columns;
+end
+if isfield(dr_in,'obs_var')
+ dr.obs_var = dr_in.obs_var;
end
if M_.exo_nbr == 0
- oo_.exo_steady_state = [] ;
+ exo_steady_state = [] ;
end
-[dr.ys,M_.params,info] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,~options_.steadystate.nocheck);
+[dr.ys,M_.params,info] = evaluate_steady_state(endo_steady_state,[exo_steady_state; exo_det_steady_state],M_,options_,~options_.steadystate.nocheck);
+params=M_.params;
if info(1)
- oo_.dr = dr;
return
end
@@ -111,5 +112,4 @@ if options_.loglinear
end
end
-[dr, info] = stochastic_solvers(dr, check_flag, M_, options_, oo_);
-oo_.dr = dr;
+[dr, info] = stochastic_solvers(dr, check_flag, M_, options_, exo_steady_state, exo_det_steady_state);
\ No newline at end of file
diff --git a/matlab/rotated_slice_sampler.m b/matlab/rotated_slice_sampler.m
index 4362015d7d2b6498560cadb22455de001866e5b1..45f5c9a41676cf2bdb745d0742231946cd5863ce 100644
--- a/matlab/rotated_slice_sampler.m
+++ b/matlab/rotated_slice_sampler.m
@@ -24,7 +24,7 @@ function [theta, fxsim, neval] = rotated_slice_sampler(objective_function,theta,
% SPECIAL REQUIREMENTS
% none
-% Copyright © 2015-2017 Dynare Team
+% Copyright © 2015-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -51,10 +51,11 @@ end
if ~isempty(sampler_options.mode)
mm = sampler_options.mode;
n = length(mm);
+ distance=NaN(n,1);
for j=1:n
distance(j)=sqrt(sum((theta-mm(j).m).^2));
end
- [m, im] = min(distance);
+ [~, im] = min(distance);
r=im;
V1 = mm(r).m;
@@ -69,34 +70,6 @@ if ~isempty(sampler_options.mode)
end
resul=randperm(n-1,n-1);
V1 = V1(:,resul);
-
- %V1 = V1(:, randperm(n-1));
- % %d = chol(mm(r).invhess);
- % %V1 = transpose(feval(sampler_options.proposal_distribution, transpose(mm(r).m), d, npar));
- %
- % V1=eye(npar);
- % V1=V1(:,randperm(npar));
- % for j=1:2,
- % V1(:,j)=mm(r(j)).m-theta;
- % V1(:,j)=V1(:,j)/norm(V1(:,j));
- % end
- % % Gram-Schmidt
- % for j=2:npar,
- % for k=1:j-1,
- % V1(:,j)=V1(:,j)-V1(:,k)'*V1(:,j)*V1(:,k);
- % end
- % V1(:,j)=V1(:,j)/norm(V1(:,j));
- % end
- % for j=1:n,
- % distance(j)=sqrt(sum((theta-mm(j).m).^2));
- % end
- % [m, im] = min(distance);
- % if im==r,
- % fxsim=[];
- % return,
- % else
- % theta1=theta;
- % end
else
V1 = sampler_options.V1;
end
@@ -106,8 +79,6 @@ for it=1:npar
theta0 = theta;
neval(it) = 0;
xold = 0;
- % XLB = thetaprior(3);
- % XUB = thetaprior(4);
tb=sort([(thetaprior(:,1)-theta)./V1(:,it) (thetaprior(:,2)-theta)./V1(:,it)],2);
XLB=max(tb(:,1));
XUB=min(tb(:,2));
@@ -172,12 +143,4 @@ for it=1:npar
end
end
end
-
-% if ~isempty(sampler_options.mode),
-% dist1=sqrt(sum((theta-mm(r).m).^2));
-% if dist1>distance(r),
-% theta=theta1;
-% fxsim=[];
-% end
-% end
end
diff --git a/matlab/save_display_classical_smoother_results.m b/matlab/save_display_classical_smoother_results.m
new file mode 100644
index 0000000000000000000000000000000000000000..898c1927fcba6235dba7ad31f385b758bf6482ef
--- /dev/null
+++ b/matlab/save_display_classical_smoother_results.m
@@ -0,0 +1,252 @@
+function oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_)
+% oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_)
+% Inputs:
+% xparam1 [double] current values for the estimated parameters.
+% M_ [structure] storing the model information
+% oo_ [structure] storing the results
+% options_ [structure] storing the options
+% bayestopt_ [structure] storing information about priors
+% dataset_ [structure] storing the dataset
+% dataset_info [structure] information about the dataset (descriptive statistics and missing observations)
+% estim_params_ [structure] storing information about estimated parameters
+%
+% Outputs:
+% oo_ [structure] storing the results
+
+% Copyright © 2023 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+
+gend= dataset_.nobs;
+% Set the number of observed variables.
+n_varobs = length(options_.varobs);
+
+if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
+ [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_);
+ if ismember(info(1),[303,304,306])
+ fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n')
+ else
+ updated_variables = atT*nan;
+ measurement_error=[];
+ ys = oo_.dr.ys;
+ trend_coeff = zeros(length(options_.varobs_id),1);
+ bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf);
+ options_nk=options_.nk;
+ options_.nk=[]; %unset options_.nk and reset it later
+ [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff);
+ options_.nk=options_nk;
+ end
+else
+ if options_.occbin.smoother.status
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info);
+ if oo_.occbin.smoother.error_flag(1)==0
+ [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
+ else
+ fprintf('\nOccbin: smoother did not succeed. No results will be written to oo_.\n')
+ end
+ else
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_);
+ [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
+ end
+ [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
+end
+if ~options_.nograph
+ [nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);
+ if ~exist([M_.dname '/graphs'],'dir')
+ mkdir(M_.dname,'graphs');
+ end
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_SmoothedShocks.tex'],'w');
+ fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
+ fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+ fprintf(fidTeX,' \n');
+ end
+ for plt = 1:nbplt
+ fh = dyn_figure(options_.nodisplay,'Name','Smoothed shocks');
+ nstar0=min(nstar,M_.exo_nbr-(plt-1)*nstar);
+ if gend==1
+ marker_string{1,1}='-ro';
+ marker_string{2,1}='-ko';
+ else
+ marker_string{1,1}='-r';
+ marker_string{2,1}='-k';
+ end
+ for i=1:nstar0
+ k = (plt-1)*nstar+i;
+ subplot(nr,nc,i);
+ plot([1 gend],[0 0],marker_string{1,1},'linewidth',.5)
+ hold on
+ plot(1:gend,innov(k,:),marker_string{2,1},'linewidth',1)
+ hold off
+ name = M_.exo_names{k};
+ if ~isempty(options_.XTick)
+ set(gca,'XTick',options_.XTick)
+ set(gca,'XTickLabel',options_.XTickLabel)
+ end
+ if gend>1
+ xlim([1 gend])
+ end
+ if options_.TeX
+ title(['$' M_.exo_names_tex{k} '$'],'Interpreter','latex')
+ else
+ title(name,'Interpreter','none')
+ end
+ end
+ dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_SmoothedShocks' int2str(plt)],options_.nodisplay,options_.graph_format);
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fprintf(fidTeX,'\\begin{figure}[H]\n');
+ fprintf(fidTeX,'\\centering \n');
+ fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_SmoothedShocks%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
+ fprintf(fidTeX,'\\caption{Smoothed shocks.}');
+ fprintf(fidTeX,'\\label{Fig:SmoothedShocks:%s}\n',int2str(plt));
+ fprintf(fidTeX,'\\end{figure}\n');
+ fprintf(fidTeX,'\n');
+ end
+ end
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fprintf(fidTeX,'\n');
+ fprintf(fidTeX,'%% End of TeX file.\n');
+ fclose(fidTeX);
+ end
+end
+if estim_params_.nvn
+ number_of_plots_to_draw = 0;
+ index = [];
+ for obs_iter=1:n_varobs
+ if max(abs(measurement_error(obs_iter,:))) > options_.ME_plot_tol
+ number_of_plots_to_draw = number_of_plots_to_draw + 1;
+ index = cat(1,index,obs_iter);
+ end
+ end
+ if ~options_.nograph
+ [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw);
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_SmoothedObservationErrors.tex'],'w');
+ fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
+ fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+ fprintf(fidTeX,' \n');
+ end
+ for plt = 1:nbplt
+ fh = dyn_figure(options_.nodisplay,'Name','Smoothed observation errors');
+ nstar0=min(nstar,number_of_plots_to_draw-(plt-1)*nstar);
+ if gend==1
+ marker_string{1,1}='-ro';
+ marker_string{2,1}='-ko';
+ else
+ marker_string{1,1}='-r';
+ marker_string{2,1}='-k';
+ end
+ for i=1:nstar0
+ k = (plt-1)*nstar+i;
+ subplot(nr,nc,i);
+ plot([1 gend],[0 0],marker_string{1,1},'linewidth',.5)
+ hold on
+ plot(1:gend,measurement_error(index(k),:),marker_string{2,1},'linewidth',1)
+ hold off
+ name = options_.varobs{index(k)};
+ if gend>1
+ xlim([1 gend])
+ end
+ if ~isempty(options_.XTick)
+ set(gca,'XTick',options_.XTick)
+ set(gca,'XTickLabel',options_.XTickLabel)
+ end
+ if options_.TeX
+ idx = strmatch(options_.varobs{index(k)}, M_.endo_names, 'exact');
+ title(['$' M_.endo_names_tex{idx} '$'],'Interpreter','latex')
+ else
+ title(name,'Interpreter','none')
+ end
+ end
+ dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_SmoothedObservationErrors' int2str(plt)],options_.nodisplay,options_.graph_format);
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fprintf(fidTeX,'\\begin{figure}[H]\n');
+ fprintf(fidTeX,'\\centering \n');
+ fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_SmoothedObservationErrors%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
+ fprintf(fidTeX,'\\caption{Smoothed observation errors.}');
+ fprintf(fidTeX,'\\label{Fig:SmoothedObservationErrors:%s}\n',int2str(plt));
+ fprintf(fidTeX,'\\end{figure}\n');
+ fprintf(fidTeX,'\n');
+ end
+ end
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fprintf(fidTeX,'\n');
+ fprintf(fidTeX,'%% End of TeX file.\n');
+ fclose(fidTeX);
+ end
+ end
+end
+%%
+%% Historical and smoothed variabes
+%%
+if ~options_.nograph
+ [nbplt,nr,nc,lr,lc,nstar] = pltorg(n_varobs);
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_HistoricalAndSmoothedVariables.tex'],'w');
+ fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
+ fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+ fprintf(fidTeX,' \n');
+ end
+ for plt = 1:nbplt
+ fh = dyn_figure(options_.nodisplay,'Name','Historical and smoothed variables');
+ nstar0=min(nstar,n_varobs-(plt-1)*nstar);
+ if gend==1
+ marker_string{1,1}='-ro';
+ marker_string{2,1}='--ko';
+ else
+ marker_string{1,1}='-r';
+ marker_string{2,1}='--k';
+ end
+ for i=1:nstar0
+ k = (plt-1)*nstar+i;
+ subplot(nr,nc,i);
+ plot(1:gend,yf(k,:),marker_string{1,1},'linewidth',1)
+ hold on
+ plot(1:gend,dataset_info.rawdata(:,k),marker_string{2,1},'linewidth',1)
+ hold off
+ name = options_.varobs{k};
+ if ~isempty(options_.XTick)
+ set(gca,'XTick',options_.XTick)
+ set(gca,'XTickLabel',options_.XTickLabel)
+ end
+ if gend>1
+ xlim([1 gend])
+ end
+ if options_.TeX
+ idx = strmatch(options_.varobs{k}, M_.endo_names,'exact');
+ title(['$' M_.endo_names_tex{idx} '$'],'Interpreter','latex')
+ else
+ title(name,'Interpreter','none')
+ end
+ end
+ dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_HistoricalAndSmoothedVariables' int2str(plt)],options_.nodisplay,options_.graph_format);
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fprintf(fidTeX,'\\begin{figure}[H]\n');
+ fprintf(fidTeX,'\\centering \n');
+ fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_HistoricalAndSmoothedVariables%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
+ fprintf(fidTeX,'\\caption{Historical and smoothed variables.}');
+ fprintf(fidTeX,'\\label{Fig:HistoricalAndSmoothedVariables:%s}\n',int2str(plt));
+ fprintf(fidTeX,'\\end{figure}\n');
+ fprintf(fidTeX,'\n');
+ end
+ end
+ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fprintf(fidTeX,'\n');
+ fprintf(fidTeX,'%% End of TeX file.\n');
+ fclose(fidTeX);
+ end
+end
diff --git a/matlab/selec_posterior_draws.m b/matlab/selec_posterior_draws.m
index 0b5dadcc5ba1b38134b581167ccb211fb7b6dfe3..72b094ab7cf2d4732b1e653b6464e67dd67092b8 100644
--- a/matlab/selec_posterior_draws.m
+++ b/matlab/selec_posterior_draws.m
@@ -1,4 +1,4 @@
-function SampleAddress = selec_posterior_draws(SampleSize,drsize)
+function SampleAddress = selec_posterior_draws(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,estim_params_,SampleSize,drsize)
% Selects a sample of draws from the posterior distribution and if nargin>1
% saves the draws in _pdraws mat files (metropolis folder). If drsize>0
% the dr structure, associated to the parameters, is also saved in _pdraws.
@@ -6,8 +6,14 @@ function SampleAddress = selec_posterior_draws(SampleSize,drsize)
% _mh file cannot be opened twice.
%
% INPUTS
-% o SampleSize [integer] Size of the sample to build.
-% o drsize [double] structure dr is drsize megaoctets.
+% o M_ [structure] Matlab's structure describing the model
+% o options_ [structure] Matlab's structure describing the current options
+% o dr [structure] Reduced form model.
+% o endo_steady_state [vector] steady state value for endogenous variables
+% o exo_steady_state [vector] steady state value for exogenous variables
+% o exo_det_steady_state [vector] steady state value for exogenous deterministic variables
+% o SampleSize [integer] Size of the sample to build.
+% o drsize [double] structure dr is drsize megaoctets.
%
% OUTPUTS
% o SampleAddress [integer] A (SampleSize*4) array, each line specifies the
@@ -37,8 +43,6 @@ function SampleAddress = selec_posterior_draws(SampleSize,drsize)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global M_ options_ estim_params_ oo_
-
% Number of parameters:
npar = estim_params_.nvx;
npar = npar + estim_params_.nvn;
@@ -48,9 +52,9 @@ npar = npar + estim_params_.np;
% Select one task:
switch nargin
- case 1
+ case 8
info = 0;
- case 2
+ case 9
MAX_mega_bytes = 10;% Should be an option...
if drsize>0
info=2;
@@ -114,8 +118,8 @@ if info
end
pdraws(i,1) = {x2(SampleAddress(i,4),:)};
if info==2
- set_parameters(pdraws{i,1});
- [dr,~,M_,oo_] =compute_decision_rules(M_,options_,oo_);
+ M_ = set_parameters_locally(M_,pdraws{i,1});
+ [dr,~,M_.params] =compute_decision_rules(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
pdraws(i,2) = { dr };
end
old_mhfile = mhfile;
@@ -141,8 +145,8 @@ if info
end
pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
if info==2
- set_parameters(pdraws{linee,1});
- [dr,~,M_,oo_] = compute_decision_rules(M_,options_,oo_);
+ M_ = set_parameters_locally(M_,pdraws{i,1});
+ [dr,~,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
pdraws(linee,2) = { dr };
end
old_mhfile = mhfile;
diff --git a/matlab/simulated_moment_uncertainty.m b/matlab/simulated_moment_uncertainty.m
index 90e58de2d5d625515c26037cd2fdea4c9fdeebd5..63a0dcefe9ae5db378b73603ef4eb2140daf55d8 100644
--- a/matlab/simulated_moment_uncertainty.m
+++ b/matlab/simulated_moment_uncertainty.m
@@ -71,8 +71,7 @@ else
logged_steady_state_indicator=0;
end
-[dr,info,M_,oo_] = compute_decision_rules(M_,options_,oo_);
-oo_.dr=dr;
+[oo_.dr,info,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if info(1)
fprintf('\nsimulated_moment_uncertainty: model could not be solved')
print_info(info,0,options_);
@@ -95,7 +94,7 @@ end
for j=1:replic
- [ys, oo_] = simult(y0,oo_.dr,M_,options_,oo_);%do simulation
+ [ys, oo_.exo_simul] = simult(y0,oo_.dr,M_,options_);%do simulation
oo_=disp_moments(ys, options_.varobs,M_,options_,oo_); %get moments
dum=[oo_.mean; dyn_vech(oo_.var)];
sd = sqrt(diag(oo_.var));
diff --git a/matlab/simult.m b/matlab/simult.m
index 737c72ee628ce480bdf741dbfa547e3d0e2427ba..7dc2db62996f42e3c180f88f675e2e4f2ed1615f 100644
--- a/matlab/simult.m
+++ b/matlab/simult.m
@@ -1,8 +1,8 @@
-function [y_out,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareResults)
+function [y_out,exo_simul] =simult(y0, dr,M_,options_)
% Simulate a DSGE model (perturbation approach).
%@info:
-%! @deftypefn {Function File} {[@var{y_}, @var{DynareResults}] =} simult (@var{y0},@var{dr},@var{DynareModel},@var{DynareOptions},@var{DynareResults})
+%! @deftypefn {Function File} {[@var{y_}, @var{exo_simul}] =} simult (@var{y0},@var{dr},@var{M_},@var{options_})
%! @anchor{simult}
%! @sp 1
%! Simulate a DSGE model (perturbation approach).
@@ -14,12 +14,10 @@ function [y_out,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareRe
%! Vector of doubles, initial conditions.
%! @item dr
%! Matlab's structure describing decision and transition rules.
-%! @item DynareModel
+%! @item M_
%! Matlab's structure describing the model (initialized by dynare, see @ref{M_})
-%! @item DynareOptions
+%! @item options_
%! Matlab's structure describing the current options (initialized by dynare, see @ref{options_}).
-%! @item DynareResults
-%! Matlab's structure gathering the results (see @ref{oo_}).
%! @end table
%! @sp 2
%! @strong{Outputs}
@@ -27,8 +25,8 @@ function [y_out,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareRe
%! @table @ @var
%! @item y_out
%! Matrix of doubles, simulated time series for all the endogenous variables (one per row).
-%! @item DynareResults
-%! Matlab's structure gathering the results (see @ref{oo_}).
+%! @item exo_simul
+%! Matrix of doubles, random exogenous shocks
%! @end table
%! @sp 2
%! @strong{This function is called by:}
@@ -45,7 +43,7 @@ function [y_out,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareRe
%! @end deftypefn
%@eod:
-% Copyright © 2001-2019 Dynare Team
+% Copyright © 2001-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -62,30 +60,30 @@ function [y_out,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareRe
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-order = DynareOptions.order;
-replic = DynareOptions.simul_replic;
+order = options_.order;
+replic = options_.simul_replic;
if replic > 1
- if ~exist([DynareModel.dname '/Output'],'dir')
- mkdir(DynareModel.dname,'Output');
+ if ~exist([M_.dname '/Output'],'dir')
+ mkdir(M_.dname,'Output');
end
- fname = [DynareModel.dname filesep 'Output' filesep DynareModel.fname,'_simul'];
+ fname = [M_.dname filesep 'Output' filesep M_.fname,'_simul'];
fh = fopen(fname,'w+');
end
% eliminate shocks with 0 variance
-i_exo_var = setdiff([1:DynareModel.exo_nbr],find(diag(DynareModel.Sigma_e) == 0));
+i_exo_var = setdiff([1:M_.exo_nbr],find(diag(M_.Sigma_e) == 0));
nxs = length(i_exo_var);
-DynareResults.exo_simul = zeros(DynareOptions.periods,DynareModel.exo_nbr);
-chol_S = chol(DynareModel.Sigma_e(i_exo_var,i_exo_var));
+exo_simul = zeros(options_.periods,M_.exo_nbr);
+chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var));
for i=1:replic
- if ~isempty(DynareModel.Sigma_e)
+ if ~isempty(M_.Sigma_e)
% we fill the shocks row wise to have the same values
% independently of the length of the simulation
- DynareResults.exo_simul(:,i_exo_var) = randn(nxs,DynareOptions.periods)'*chol_S;
+ exo_simul(:,i_exo_var) = randn(nxs,options_.periods)'*chol_S;
end
- y_ = simult_(DynareModel,DynareOptions,y0,dr,DynareResults.exo_simul,order);
+ y_ = simult_(M_,options_,y0,dr,exo_simul,order);
% elimninating initial value
y_ = y_(:,2:end);
if replic > 1
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 0742ddc47ab24598b9fd2ad54d0ea033c5ff5e5e..6ea2f14d1b44eaf080396a8f5464c29e8ae0b9e0 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -59,8 +59,6 @@ end
test_for_deep_parameters_calibration(M_);
-dr = oo_.dr;
-
options_old = options_;
if options_.linear
options_.order = 1;
@@ -105,7 +103,7 @@ end
check_model(M_);
-oo_.dr=set_state_space(dr,M_,options_);
+oo_.dr=set_state_space(oo_.dr,M_,options_);
if PI_PCL_solver
[oo_.dr, info] = PCL_resol(oo_.steady_state,0);
@@ -113,14 +111,14 @@ elseif options_.discretionary_policy
if ~options_.order==1
error('discretionary_policy: only linear-quadratic problems can be solved');
end
- [~,info,M_,oo_] = discretionary_policy_1(options_.instruments,M_,options_,oo_);
+ [oo_.dr,info,M_.params] = discretionary_policy_1(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
else
if options_.logged_steady_state %if steady state was previously logged, undo this
oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state);
options_.logged_steady_state=0;
end
- [~,info,M_,oo_] = resol(0,M_,options_,oo_);
+ [oo_.dr,info,M_.params] = resol(0,M_,options_,oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
end
if options_.loglinear && isfield(oo_.dr,'ys') && options_.logged_steady_state==0 %log steady state for correct display of decision rule
@@ -209,8 +207,7 @@ if options_.periods > 0 && ~PI_PCL_solver
y0 = M_.endo_histval;
end
end
- [ys, oo_] = simult(y0,oo_.dr,M_,options_,oo_);
- oo_.endo_simul = ys;
+ [oo_.endo_simul, oo_.exo_simul] = simult(y0,oo_.dr,M_,options_);
if ~options_.minimal_workspace
dyn2vec(M_, oo_, options_);
end
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index 6b0c4a55e8842c27bb819f555d37941ea65ddb00..3234681f7d96ec5bf64914d3b94c662a1ff2ffae 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -1,5 +1,5 @@
-function [dr, info] = stochastic_solvers(dr, task, M_, options_, oo_)
-
+function [dr, info] = stochastic_solvers(dr, task, M_, options_, exo_steady_state, exo_det_steady_state)
+%[dr, info] = stochastic_solvers(dr, task, M_, options_, exo_steady_state, exo_det_steady_state)
% Computes the reduced form solution of a rational expectations model (first, second or third
% order approximation of the stochastic model around the deterministic steady state).
%
@@ -8,7 +8,8 @@ function [dr, info] = stochastic_solvers(dr, task, M_, options_, oo_)
% - task [integer] scalar, if task = 0 then decision rules are computed and if task = 1 then only eigenvales are computed.
% - M_ [struct] Definition of the model.
% - options_ [struct] Options.
-% - oo_ [struct] Results
+% - exo_steady_state [vector] steady state value for exogenous variables
+% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - dr [struct] Decision rules for stochastic simulations.
@@ -93,15 +94,11 @@ if options_.k_order_solver
end
klen = M_.maximum_lag + M_.maximum_lead + 1;
-exo_simul = [repmat(oo_.exo_steady_state',klen,1) repmat(oo_.exo_det_steady_state',klen,1)];
+exo_simul = [repmat(exo_steady_state',klen,1) repmat(exo_det_steady_state',klen,1)];
iyv = M_.lead_lag_incidence';
iyv = iyv(:);
iyr0 = find(iyv) ;
-if M_.exo_nbr == 0
- oo_.exo_steady_state = [] ;
-end
-
it_ = M_.maximum_lag + 1;
z = repmat(dr.ys,1,klen);
if local_order == 1
diff --git a/preprocessor b/preprocessor
index fe73d74f8a0f8cafda7062f6ec883bf831a30f50..097c4fdc77f865c40b1b4a6d625edaec39341819 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit fe73d74f8a0f8cafda7062f6ec883bf831a30f50
+Subproject commit 097c4fdc77f865c40b1b4a6d625edaec39341819
diff --git a/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod b/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod
index 068c46be991f4711101921fdb8edf0bebf534855..e1acfb16a3ea61e20778b2a39b285c44b5e92830 100644
--- a/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod
+++ b/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod
@@ -433,7 +433,7 @@ for jj = 1:2
xparam1_calib = [xparam1_calib; calib_params(indpmodel)];
M_ = set_all_parameters(xparam1_calib,estim_params_,M_);
end
- [~,info,M_,oo_] = resol(0,M_, options_, oo_);
+ [oo_.dr,info,M_.params] = resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
%For convenience we save the objects to compare into oo_.dr.
oo_.dr.Yss = oo_.dr.ys(oo_.dr.order_var);
oo_.dr.Sigma_e = M_.Sigma_e;
diff --git a/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod b/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod
index b868e7b15e8422d8ea06efc43889574f074490f7..da9951e128a8ee6e32acc3a337c4a2c76182ede6 100644
--- a/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod
+++ b/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod
@@ -115,7 +115,7 @@ identification(order=@{ORDER},nograph,no_identification_strength);
%make sure everything is computed at prior mean
xparam_prior = set_prior(estim_params_,M_,options_);
M_ = set_all_parameters(xparam_prior,estim_params_,M_);
-[~,info,M_,oo_] = resol(0,M_, options_, oo_);
+[oo_.dr,info,M_.params] = resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
indpmodel = estim_params_.param_vals(:,1);
indpstderr = estim_params_.var_exo(:,1);
diff --git a/tests/estimation/fs2000.mod b/tests/estimation/fs2000.mod
index eaaddf7de6dfd41901364ca6b37d6054486284e1..13f686415ba3b864c115b74591f068d8c4729dbd 100644
--- a/tests/estimation/fs2000.mod
+++ b/tests/estimation/fs2000.mod
@@ -132,4 +132,4 @@ oo_.MarginalDensity.LaplaceApproximation = Laplace; %reset correct Laplace
%test prior sampling
options_.prior_draws=100;
options_.no_graph.posterior=0;
-prior_posterior_statistics('prior',dataset_,dataset_info); %get smoothed and filtered objects and forecasts
+oo_=prior_posterior_statistics('prior',dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_); %get smoothed and filtered objects and forecasts
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod
index b0621a82bd0e0f8ef8e6bf5b2199a36a233ea08f..ac9acc7e8c24813296e88484edfc5d6a343b121f 100644
--- a/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod
@@ -129,7 +129,7 @@ y0 = [y0; oo_.SmoothedVariables.Mean.(M_.endo_names{endo_iter})(1)];
end;
%make sure decision rules were updated
-[oo_.dr,info,M_] = resol(0,M_,options_,oo_);
+[oo_.dr,info,M_.params] = resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
dr = oo_.dr;
iorder=1;
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod
index bc2f652f37984c3ba9c79337c39c155998621c07..4f92baa837fbaff56880d66c9f9886f1d8628b4d 100644
--- a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod
@@ -132,7 +132,7 @@ y0 = [y0; oo_.SmoothedVariables.(M_.endo_names{endo_iter})(1)];
end;
%make sure decision rules were updated
-[oo_.dr,info,M_] = resol(0,M_,options_,oo_);
+[oo_.dr,info,M_.params] = resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
dr = oo_.dr;
iorder=1;
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod
index ea87ccfe537001547bfdeb72b6382c8c5dac3fca..53c3881fa6e708f57d877186e4493ed1e43bcdba 100644
--- a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod
@@ -132,7 +132,7 @@ y0 = [y0; oo_.SmoothedVariables.(M_.endo_names{endo_iter})(1)];
end;
%make sure decision rules were updated
-[oo_.dr,info,M_] = resol(0,M_,options_,oo_);
+[oo_.dr,info,M_.params] = resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
dr = oo_.dr;
iorder=1;
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod
index 0277c9dcecaecb82c57c03b616c7a4de7d4720c3..f30d9d3fc9841e3fb801c5f00b0098eff161f63b 100644
--- a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod
@@ -146,7 +146,7 @@ y0 = [y0; oo_.SmoothedVariables.Mean.(M_.endo_names{endo_iter})(1)];
end;
%make sure decision rules were updated
-[oo_.dr,info,M_] = resol(0,M_,options_,oo_);
+[oo_.dr,info,M_.params] = resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
dr = oo_.dr;
iorder=1;
diff --git a/tests/kalman_filter_smoother/test_compute_Pinf_Pstar.mod b/tests/kalman_filter_smoother/test_compute_Pinf_Pstar.mod
index 24f2db61515f1880e7ec57f9fa9434316667ee01..cec2dea32da2dc72e6a9e2ae32bd1399ed120be3 100644
--- a/tests/kalman_filter_smoother/test_compute_Pinf_Pstar.mod
+++ b/tests/kalman_filter_smoother/test_compute_Pinf_Pstar.mod
@@ -95,7 +95,7 @@ calib_smoother(datafile=data_Pinf_Pstar,diffuse_filter);
mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf);
Q = M_.Sigma_e;
-[T,R,SteadyState,info,M_,oo_] = dynare_resolve(M_,options_,oo_);
+[T,R,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium);