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

Merge branch 'more_globals' into 'master' 

Eliminate some more globals

See merge request Dynare/dynare!2187
kalman-mex
Sébastien Villemot 2023-10-03 12:48:34 +00:00
parent da9789038d 7ff7561280
commit e4afa5eee0
58 changed files with 1028 additions and 990 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
examples/Gali_2015_prior_restrictions.m
View File

@ -1,12 +1,12 @@
function log_prior_val=Gali_2015_prior_restrictions(M_, oo_, options_, dataset_, dataset_info);
% function prior_val=Gali_2015_prior_restrictions(M_, oo_, options_, dataset_, dataset_info);
function log_prior_val=Gali_2015_prior_restrictions(M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_, dataset_, dataset_info);
% function prior_val=Gali_2015_prior_restrictions(M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_, dataset_, dataset_info);
% Example of a _prior_restrictions-file automatically called during
% estimation
% It imposes a prior of the slope of the New Keynesian Phillips Curve of
% 0.03. As the slope is a composite of other parameters with independent
% priors, a separate function is required to do this.
% Copyright © 2021 Dynare Team
% Copyright © 2021-2023 Dynare Team
%
% This file is part of Dynare.
%

2
matlab/+mom/Jtest.m
View File

@ -49,7 +49,7 @@ if options_mom_.mom.mom_nbr > length(xparam)
W_opt = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
oo_J = oo_;
oo_J.mom.Sw = chol(W_opt);
fval = feval(objective_function, xparam, Bounds, oo_J, estim_params_, M_, options_mom_, bayestopt_);
fval = feval(objective_function, xparam, Bounds, oo_J, estim_params_, M_, options_mom_);
else
fval = oo_.mom.Q;
end

11
matlab/+mom/get_data_moments.m
View File

@ -1,10 +1,11 @@
function [dataMoments, m_data] = get_data_moments(data, oo_, matched_moments_, options_mom_)
% [dataMoments, m_data] = get_data_moments(data, oo_, matched_moments_, options_mom_)
function [dataMoments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_, options_mom_)
% [dataMoments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_, options_mom_)
% This function computes the user-selected empirical moments from data
% =========================================================================
% INPUTS
% o data [T x varobs_nbr] data set
% o oo_: [structure] storage for results
% o obs_var: [integer] index of observables
% o inv_order_var: [integer] inverse decision rule order
% o matched_moments_: [structure] information about selected moments to match in estimation
% o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% -------------------------------------------------------------------------
@ -45,11 +46,11 @@ m_data = NaN(T,options_mom_.mom.mom_nbr);
% Product moment for each time period, i.e. each row t contains y_t1(l1)^p1*y_t2(l2)^p2*...
% note that here we already are able to treat leads and lags and any power product moments
for jm = 1:options_mom_.mom.mom_nbr
vars = oo_.dr.inv_order_var(matched_moments_{jm,1})';
vars = inv_order_var(matched_moments_{jm,1})';
leadlags = matched_moments_{jm,2}; % lags are negative numbers and leads are positive numbers
powers = matched_moments_{jm,3};
for jv = 1:length(vars)
jvar = (oo_.mom.obs_var == vars(jv));
jvar = (obs_var == vars(jv));
y = NaN(T,1); %Take care of T_eff instead of T for lags and NaN via mean with 'omitnan' option below
y( (1-min(leadlags(jv),0)) : (T-max(leadlags(jv),0)), 1) = data( (1+max(leadlags(jv),0)) : (T+min(leadlags(jv),0)), jvar).^powers(jv);
if jv==1

9
matlab/+mom/objective_function.m
View File

@ -1,5 +1,5 @@
function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_)
% [fval, info, exit_flag, df, junk1, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_)
function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_)
% [fval, info, exit_flag, df, junk1, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_)
% -------------------------------------------------------------------------
% This function evaluates the objective function for method of moments estimation
% =========================================================================
@ -10,7 +10,6 @@ function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(
% o estim_params_: [structure] describing the estimated_parameters
% o M_ [structure] describing the model
% o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% o bayestopt_: [structure] information about the prior
% -------------------------------------------------------------------------
% OUTPUTS
% o fval: [double] value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
@ -151,7 +150,7 @@ 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;
oo_.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_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, 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_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 1);
else
@ -247,7 +246,7 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
if options_mom_.prefilter
y_sim = bsxfun(@minus, y_sim, mean(y_sim,1));
end
oo_.mom.model_moments = mom.get_data_moments(y_sim, oo_, M_.matched_moments, options_mom_);
oo_.mom.model_moments = mom.get_data_moments(y_sim, oo_.mom.obs_var, oo_.dr.inv_order_var, M_.matched_moments, options_mom_);
end

6
matlab/+mom/run.m
View File

@ -435,7 +435,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
% Provide info on data moments handling
fprintf('Computing data moments. Note that NaN values in the moments (due to leads and lags or missing data) are replaced by the mean of the corresponding moment.\n');
% Get data moments for the method of moments
[oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_);
[oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, oo_.mom.obs_var, oo_.dr.inv_order_var, M_.matched_moments, options_mom_);
if ~isreal(dataset_.data)
error('method_of_moments: The data moments contain complex values!')
end
@ -509,7 +509,7 @@ try
oo_.mom.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix
end
tic_id = tic;
[fval, info, ~, ~, ~, oo_, M_] = feval(objective_function, xparam0, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_);
[fval, info, ~, ~, ~, oo_, M_] = feval(objective_function, xparam0, Bounds, oo_, estim_params_, M_, options_mom_);
elapsed_time = toc(tic_id);
if isnan(fval)
error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!')
@ -569,7 +569,7 @@ end
% -------------------------------------------------------------------------
if options_mom_.mode_check.status
mode_check(objective_function, xparam1, hessian_xparam1, options_mom_, M_, estim_params_, bayestopt_, Bounds, true,...
Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_);
Bounds, oo_, estim_params_, M_, options_mom_);
end

17
matlab/+occbin/IVF_core.m
View File

@ -1,6 +1,6 @@
function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% Computes thre
function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% Calls the solver to get the shocks explaining the data for the inversion filter
%
% Outputs:
% - filtered_errs [T by N_obs] filtered shocks
@ -11,7 +11,10 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o
%
% Inputs
% - M_ [structure] Matlab's structure describing the model (M_).
% - 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
% - options_ [structure] Matlab's structure describing the current options (options_).
% - err_index [double] index of shocks with strictly positive variance in M_.exo_names
% - filtered_errs_init [T by N_obs] initial values for the shocks
@ -39,7 +42,7 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,o
[sample_length, n_obs]= size(obs);
nerrs = size(err_index,1);
if nargin<8
if nargin<11
init_val = zeros(M_.endo_nbr,1);
end
@ -83,7 +86,7 @@ for this_period=1:sample_length
opts_simul.exo_pos=err_index(inan); %err_index is predefined mapping from observables to shocks
opts_simul.endo_init = init_val_old;
opts_simul.SHOCKS = filtered_errs_init(this_period,inan);
[err_vals_out, exitflag] = dynare_solve(@occbin.match_function, filtered_errs_init(this_period,inan)', options_.occbin.solver.maxit, options_.occbin.solver.solve_tolf, options_.occbin.solver.solve_tolx, options_, obs_list, current_obs, opts_simul, M_, oo_, options_);
[err_vals_out, exitflag] = dynare_solve(@occbin.match_function, filtered_errs_init(this_period,inan)', options_.occbin.solver.maxit, options_.occbin.solver.solve_tolf, options_.occbin.solver.solve_tolx, options_, obs_list, current_obs, opts_simul, M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_);
if exitflag
filtered_errs=NaN;
@ -97,7 +100,7 @@ for this_period=1:sample_length
opts_simul.SHOCKS = err_vals_out;
[ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_] = occbin.match_function(...
err_vals_out,obs_list,current_obs,opts_simul, M_,oo_,options_);
err_vals_out,obs_list,current_obs,opts_simul, M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_);
init_val = stateval(this_period,:); %update
if max(abs(err_vals_out))>1e8
error_code(1) = 306;

63
matlab/+occbin/IVF_posterior.m
View File

@ -1,18 +1,23 @@
function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults, atT, innov] = IVF_posterior(xparam1,...
dataset_,obs_info,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults)
% function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults, atT, innov] = IVF_posterior(xparam1,...
% dataset_,obs_info,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults)
function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
% dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% Computes Likelihood with inversion filter
%
% INPUTS
% - xparam1 [double] current values for the estimated parameters.
% - dataset_ [structure] dataset after transformations
% - DynareOptions [structure] Matlab's structure describing the current options (options_).
% - Model [structure] Matlab's structure describing the model (M_).
% - EstimatedParameters [structure] characterizing parameters to be estimated
% - BayesInfo [structure] describing the priors
% - dataset_info [structure] storing informations about the
% sample; not used but required for interface
% - options_ [structure] Matlab's structure describing the current options
% - M_ [structure] Matlab's structure describing the model
% - estim_params_ [structure] characterizing parameters to be estimated
% - bayestopt_ [structure] describing the priors
% - BoundsInfo [structure] containing prior bounds
% - DynareResults [structure] Matlab's structure containing the results (oo_).
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - fval [double] scalar, value of the likelihood or posterior kernel.
@ -22,14 +27,14 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti
% - Hess [double] Empty array.
% - SteadyState [double] Empty array.
% - trend [double] Empty array.
% - Model [struct] Updated Model structure described in INPUTS section.
% - DynareOptions [struct] Updated DynareOptions structure described in INPUTS section.
% - BayesInfo [struct] See INPUTS section.
% - DynareResults [struct] Updated DynareResults structure described in INPUTS section.
% - M_ [struct] Updated M_ structure described in INPUTS section.
% - options_ [struct] Updated options_ structure described in INPUTS section.
% - bayestopt_ [struct] See INPUTS section.
% - dr [structure] Reduced form model.
% - atT [double] (m*T) matrix, smoothed endogenous variables (a_{t|T}) (decision-rule order)
% - innov [double] (r*T) matrix, smoothed structural shocks (r>n is the umber of shocks).
% Copyright © 2021 Dynare Team
% Copyright © 2021-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -51,7 +56,7 @@ DLIK=[];
Hess=[];
trend_coeff = [];
obs = dataset_.data;
obs_list = DynareOptions.varobs(:);
obs_list = options_.varobs(:);
exit_flag = 1;
@ -65,18 +70,18 @@ end
%------------------------------------------------------------------------------
if ~isempty(xparam1)
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
[fval,info,exit_flag,Q,H]=check_bounds_and_definiteness_estimation(xparam1, Model, EstimatedParameters, BoundsInfo);
M_ = set_all_parameters(xparam1,estim_params_,M_);
[fval,info,exit_flag,Q,H]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, BoundsInfo);
if info(1)
return
end
end
err_index=DynareOptions.occbin.likelihood.IVF_shock_observable_mapping; % err_index= find(diag(Model.Sigma_e)~=0);
COVMAT1 = Model.Sigma_e(err_index,err_index);
err_index=options_.occbin.likelihood.IVF_shock_observable_mapping; % err_index= find(diag(M_.Sigma_e)~=0);
COVMAT1 = M_.Sigma_e(err_index,err_index);
% Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
[T,R,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict');
[T,R,SteadyState,info,dr, M_.params] = dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
@ -102,17 +107,17 @@ end
sample_length = size(obs,1);
filtered_errs_init = zeros(sample_length,sum(err_index));
[filtered_errs, resids, Emat, stateval, info] = occbin.IVF_core(Model,DynareResults,DynareOptions,err_index,filtered_errs_init,obs_list,obs);
[filtered_errs, resids, Emat, stateval, info] = occbin.IVF_core(M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,err_index,filtered_errs_init,obs_list,obs);
if info(1)
fval = Inf;
exit_flag = 0;
atT=NaN(size(stateval(:,DynareResults.dr.order_var)'));
innov=NaN(Model.exo_nbr,sample_length);
atT=NaN(size(stateval(:,dr.order_var)'));
innov=NaN(M_.exo_nbr,sample_length);
return
else
atT = stateval(:,DynareResults.dr.order_var)';
innov = zeros(Model.exo_nbr,sample_length);
innov(diag(Model.Sigma_e)~=0,:)=filtered_errs';
atT = stateval(:,dr.order_var)';
innov = zeros(M_.exo_nbr,sample_length);
innov(diag(M_.Sigma_e)~=0,:)=filtered_errs';
end
nobs=size(filtered_errs,1);
@ -121,7 +126,7 @@ nobs=size(filtered_errs,1);
% Calculate the selection matrix
%-------------------------------------
iobs=DynareOptions.varobs_id';
iobs=options_.varobs_id';
likei=NaN(nobs,1);
if ~any(any(isnan(obs)))
@ -145,7 +150,7 @@ else
end
end
like = 0.5*sum(likei(DynareOptions.presample+1:end));
like = 0.5*sum(likei(options_.presample+1:end));
if isinf(like)
fval = Inf;
@ -172,7 +177,7 @@ if maxresid>1e-3
end
if ~isempty(xparam1)
prior = -priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
prior = -priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
else
prior = 0;
end

26
matlab/+occbin/dynare_resolve.m
View File

@ -1,14 +1,17 @@
function [A,B,ys,info,dr,params,TT, RR, CC, A0, B0] ...
= dynare_resolve(M_,options_,oo_,regime_history, reduced_state_space, A, B)
% function [A,B,ys,info,M_,options_,oo_,TT, RR, CC, A0, B0] ...
% = dynare_resolve(M_,options_,oo_,regime_history, reduced_state_space, A, B)
= dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,regime_history, reduced_state_space, A, B)
% [A,B,ys,info,M_,options_,oo_,TT, RR, CC, A0, B0] ...
% = dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,regime_history, reduced_state_space, A, B)
% Computes the linear approximation and the matrices A and B of the
% transition equation. Mirrors dynare_resolve
%
% Inputs:
% - M_ [structure] Matlab's structure describing the model
% - options_ [structure] Matlab's structure containing the options
% - oo_ [structure] Matlab's structure containing the results
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - reduced_state_space [string]
% - A [double] State transition matrix
% - B [double] shock impact matrix
@ -26,7 +29,7 @@ function [A,B,ys,info,dr,params,TT, RR, CC, A0, B0] ...
% - A0 [double] State transition matrix (unrestricted state space)
% - B0 [double] shock impact matrix (unrestricted state space)
% Copyright © 2001-2021 Dynare Team
% Copyright © 2001-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -43,19 +46,18 @@ function [A,B,ys,info,dr,params,TT, RR, CC, A0, B0] ...
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if nargin<6
[A,B,ys,info,dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if nargin<9
[A,B,ys,info,dr,M_.params] = dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
else
ys = oo_.dr.ys;
dr = oo_.dr;
ys = dr.ys;
info = 0;
end
params=M_.params;
if ~info(1) && nargin>4 && ~isempty(regime_history)
if ~info(1) && nargin>7 && ~isempty(regime_history)
opts_regime.regime_history=regime_history;
opts_regime.binding_indicator=[];
[TT, RR, CC] = ...
occbin.check_regimes(A, B, [], opts_regime, M_, options_, dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
occbin.check_regimes(A, B, [], opts_regime, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
else
TT=A;
RR=B;
@ -64,7 +66,7 @@ end
A0=A;
B0=B;
if ~info(1) && nargin>4 && ischar(reduced_state_space) && ~isempty(reduced_state_space)
if ~info(1) && nargin>7 && ischar(reduced_state_space) && ~isempty(reduced_state_space)
iv = dr.restrict_var_list;
A=A(iv,iv);
B=B(iv,:);

51
matlab/+occbin/kalman_update_algo_1.m
View File

@ -1,5 +1,5 @@
function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,oo_,options_,occbin_options)
% function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,oo_,options_,occbin_options)
function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options)
% function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options)
% INPUTS
% - a [N by 1] t-1's state estimate
% - a1 [N by 2] state predictions made at t-1:t
@ -18,7 +18,10 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal
% - CC [N by 2] state space constant state transition matrix at t-1:t
% - regimes0 [structure] regime info at t-1:t
% - M_ [structure] Matlab's structure describing the model (M_).
% - oo_ [structure] Matlab's structure containing the results (oo_).
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - options_ [structure] Matlab's structure describing the current options (options_).
% - occbin_options_ [structure] Matlab's structure describing the Occbin options.
% - kalman_tol [double] tolerance for reciprocal condition number
@ -42,7 +45,7 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kal
% Philipp Pfeiffer, Marco Ratto (2021), Efficient and robust inference of models with occasionally binding
% constraints, Working Papers 2021-03, Joint Research Centre, European Commission
% Copyright © 2021-2022 Dynare Team
% Copyright © 2021-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -104,7 +107,7 @@ if ~options_.occbin.filter.use_relaxation
[a, a1, P, P1, v, alphahat, etahat, lik, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
else
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, base_regime,'reduced_state_space',T0,R0);
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, base_regime,'reduced_state_space',T0,R0);
regimes0(1)=base_regime;
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
@ -123,18 +126,18 @@ opts_simul.exo_pos=1:M_.exo_nbr;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
my_order_var = oo_.dr.order_var(oo_.dr.restrict_var_list);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
my_order_var = dr.order_var(dr.restrict_var_list);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
my_order_var = oo_.dr.order_var;
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
my_order_var = dr.order_var;
end
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
if out.error_flag
error_flag = out.error_flag;
@ -188,7 +191,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
end
regimes_(1).regimestart(end)=regimestart;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],'reduced_state_space', T0, R0);
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, [base_regime regimes_(1)],'reduced_state_space', T0, R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
@ -208,10 +211,10 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
if not(options_.occbin.filter.use_relaxation)
opts_simul.init_regime=regimes_(1);
@ -223,7 +226,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_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
@ -254,7 +257,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
regimes_(1).regimestart=[1 2];
end
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],'reduced_state_space',T0,R0);
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, [base_regime regimes_(1)],'reduced_state_space',T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
@ -273,7 +276,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_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
@ -307,7 +310,7 @@ if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~
regimes_(1).regimestart(end)=k;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],'reduced_state_space',T0,R0);
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, [base_regime regimes_(1)],'reduced_state_space',T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
@ -318,10 +321,10 @@ if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
% opts_simul.init_regime=regimes_(1);
if M_.occbin.constraint_nbr==1
@ -331,7 +334,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_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);

59
matlab/+occbin/kalman_update_algo_3.m
View File

@ -1,5 +1,5 @@
function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,oo_,options_,occbin_options,kalman_tol,nk)
% function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,oo_,options_,occbin_options,kalman_tol,nk)
function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options,kalman_tol,nk)
% function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options,kalman_tol,nk)
%
% INPUTS
% - a [N by 1] t-1's state estimate
@ -22,7 +22,10 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat,
% - regimes0 [structure] regime info at t-1:t
% - M_ [structure] Matlab's structure describing the model (M_).
% - options_ [structure] Matlab's structure describing the current options (options_).
% - oo_ [structure] Matlab's structure containing the results (oo_).
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - occbin_options_ [structure] Matlab's structure describing the Occbin options.
% - kalman_tol [double] tolerance for reciprocal condition number
% - nk [double] number of forecasting periods
@ -52,7 +55,7 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat,
% constraints, Working Papers 2021-03, Joint Research Centre, European Commission
% Copyright © 2021 Dynare Team
% Copyright © 2021-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -104,7 +107,7 @@ if ~options_.occbin.filter.use_relaxation
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
else
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, base_regime,myrestrict,T0,R0);
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, base_regime,myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
@ -120,16 +123,16 @@ opts_simul.exo_pos=1:M_.exo_nbr;
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
my_order_var = oo_.dr.order_var(oo_.dr.restrict_var_list);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
my_order_var = dr.order_var(dr.restrict_var_list);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
my_order_var = oo_.dr.order_var;
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
my_order_var = dr.order_var;
end
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
return;
@ -183,7 +186,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
% % % regimestart = regimes_(1).regimestart(end-1)+round(0.5*(newstart+oldstart))-1;
regimes_(1).regimestart(end)=regimestart;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],myrestrict,T0,R0);
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
@ -201,15 +204,15 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.SHOCKS(1,:) = etahat(:,2)';
% if opts_simul.restrict_state_space
% tmp=zeros(M_.endo_nbr,1);
% tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
% opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
% tmp(dr.restrict_var_list,1)=alphahat(:,1);
% opts_simul.endo_init = tmp(dr.inv_order_var,1);
% else
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
% end
if not(options_.occbin.filter.use_relaxation)
@ -222,7 +225,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_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
return;
@ -252,7 +255,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
regimes_(1).regimestart=[1 2];
end
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],myrestrict,T0,R0);
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
@ -271,7 +274,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_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
end
end
@ -301,7 +304,7 @@ if error_flag==0 && niter>options_.occbin.likelihood.max_number_of_iterations &&
regimes_(1).regimestart(end)=k;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,oo_, [base_regime regimes_(1)],myrestrict,T0,R0);
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
@ -309,10 +312,10 @@ if error_flag==0 && niter>options_.occbin.likelihood.max_number_of_iterations &&
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
@ -321,7 +324,7 @@ if error_flag==0 && 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_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if isequal(out.regime_history(1),regimes_(1))
error_flag=0;
break
@ -336,9 +339,9 @@ a = out.piecewise(1:nk+1,my_order_var)' - repmat(out.ys(my_order_var),1,nk+1);
T = ss.T(my_order_var,my_order_var,:);
R = ss.R(my_order_var,:,:);
C = ss.C(my_order_var,:);
TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,1);
RR = ss.R(oo_.dr.order_var,:,1);
CC = ss.C(oo_.dr.order_var,1);
TT = ss.T(dr.order_var,dr.order_var,1);
RR = ss.R(dr.order_var,:,1);
CC = ss.C(dr.order_var,1);
QQ = R(:,:,2)*QQQ(:,:,3)*transpose(R(:,:,2));
P(:,:,1) = P(:,:,2);
for j=1:nk

11
matlab/+occbin/match_function.m
View File

@ -1,7 +1,7 @@
function [resids, grad, state_out, E, M_, out] = match_function(err_0, obs_list,current_obs, opts_simul,...
M_, oo_, options_)
M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_)
% function [resids, grad, stateout, E, M_, out] = match_function(err_0, obs_list,current_obs, opts_simul,...
% M_, oo_, options_)
% M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_)
% Outputs:
% - resids [double] [n_exo by 1] vector of residuals
% - grad [double] [n by n_exo] gradient (response of observables to shocks)
@ -16,7 +16,10 @@ function [resids, grad, state_out, E, M_, out] = match_function(err_0, obs_list,
% - current_obs [double] [1 by n_obs] current value of observables
% - opts_simul [structure] Structure with simulation options
% - M_ [structure] Matlab's structure describing the model (M_).
% - 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
% - options_ [structure] Matlab's structure describing the current options (options_).
% Original authors: Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello, and Molin Zhong
@ -36,7 +39,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_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
nobs = size(obs_list,1);
resids = zeros(nobs,1);

2
matlab/+osr/objective.m
View File

@ -64,7 +64,7 @@ if ~options_.analytic_derivation
loss = full(weights(:)'*vx(:));
else
totparam_nbr=length(i_params);
oo_.dr.derivs = get_perturbation_params_derivs(M_, options_, [], oo_, i_params, [], [], 0); %analytic derivatives of perturbation matrices
oo_.dr.derivs = get_perturbation_params_derivs(M_, options_, [], oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, i_params, [], [], 0); %analytic derivatives of perturbation matrices
pruned_state_space = pruned_state_space_system(M_, options_, oo_.dr, i_var, 0, 0, 1);
vx = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y';

5
matlab/DsgeSmoother.m
View File

@ -1,4 +1,5 @@
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)
% [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
@ -115,7 +116,7 @@ else
[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,oo_.dr, M_.params,~,~,~, T0, R0] = ...
occbin.dynare_resolve(M_,options_,oo_,[],'restrict');
occbin.dynare_resolve(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state,[],'restrict');
varargin{length_varargin+1}=T0;
varargin{length_varargin+2}=R0;
end
@ -123,7 +124,7 @@ else
end
if options_.occbin.smoother.status
occbin_info.status = true;
occbin_info.info= [{options_,oo_,M_} varargin];
occbin_info.info= [{options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,M_} varargin];
else
occbin_info.status = false;
end

2
matlab/PosteriorIRF_core1.m
View File

@ -198,7 +198,7 @@ while fpar<B
end
if MAX_nirfs_dsgevar
IRUN = IRUN+1;
[~,~,~,~,~,~,~,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
[~,~,~,~,~,~,~,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight', M_.param_names));
DSGE_PRIOR_WEIGHT = floor(dataset_.nobs*(1+dsge_prior_weight));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.nobs*(dsge_prior_weight+1)));

2
matlab/SMC_samplers_initialization.m
View File

@ -85,7 +85,7 @@ fprintf(fidlog,'%% Session 1.\n');
fprintf(fidlog,' \n');
prior_draw(bayestopt_,options_.prior_trunc);
% Find initial values for the NumberOfParticles chains...
set_dynare_seed('default');
options_=set_dynare_seed_local_options(options_,'default');
fprintf(fidlog,[' Initial values of the parameters:\n']);
disp('Estimation::dsmh: Searching for initial values...');
ix2 = zeros(npar,NumberOfParticles);

2
matlab/backward/simul_backward_model_init.m
View File

@ -155,7 +155,7 @@ if nargin<6 || isempty(innovations)
covariance_matrix_upper_cholesky = chol(covariance_matrix);
% Set seed to its default state.
if DynareOptions.bnlms.set_dynare_seed_to_default
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
end
% Simulate structural innovations.
switch DynareOptions.bnlms.innovation_distribution

6
matlab/cli/prior.m
View File

@ -1,5 +1,5 @@
function varargout = prior(varargin)
% varargout = prior(varargin)
% Computes various prior statistics and display them in the command window.
%
% INPUTS
@ -11,7 +11,7 @@ function varargout = prior(varargin)
% SPECIAL REQUIREMENTS
% none
% Copyright © 2015-2019 Dynare Team
% Copyright © 2015-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -139,7 +139,7 @@ if ismember('moments', varargin) % Prior simulations (2nd order moments).
% Solve model
[T,R,~,info,oo__.dr, Model.params] = dynare_resolve(Model , options_ , oo__.dr, oo__.steady_state, oo__.exo_steady_state, oo__.exo_det_steady_state,'restrict');
if ~info(1)
info=endogenous_prior_restrictions(T,R,Model , options__ , oo__);
info=endogenous_prior_restrictions(T,R,Model , options__ , oo__.dr,oo__.steady_state,oo__.exo_steady_state,oo__.exo_det_steady_state);
end
if info
skipline()

93
matlab/dsge_conditional_likelihood_1.m
View File

@ -1,5 +1,36 @@
function [fval, info, exitflag, DLIK, Hess, SteadyState, trend_coeff, Model, DynareOptions, BayesInfo, DynareResults] = ...
dsge_conditional_likelihood_1(xparam1, DynareDataset, DatasetInfo, DynareOptions, Model, EstimatedParameters, BayesInfo, BoundsInfo, DynareResults, derivatives_info)
function [fval, info, exitflag, DLIK, Hess, SteadyState, trend_coeff, M_, options_, bayestopt_, dr] = ...
dsge_conditional_likelihood_1(xparam1, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, derivatives_info)
% [fval, info, exitflag, DLIK, Hess, SteadyState, trend_coeff, M_, options_, bayestopt_, dr] = ...
% dsge_conditional_likelihood_1(xparam1, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, derivatives_info)
%
% INPUTS
% - xparam1 [double] current values for the estimated parameters.
% - dataset_ [structure] dataset after transformations
% - dataset_info [structure] storing informations about the
% sample; not used but required for interface
% - options_ [structure] Matlab's structure describing the current options
% - M_ [structure] Matlab's structure describing the model
% - estim_params_ [structure] characterizing parameters to be estimated
% - bayestopt_ [structure] describing the priors
% - BoundsInfo [structure] containing prior bounds
% - 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
% - derivatives_info [structure] derivative info for identification
%
% OUTPUTS
% - fval [double] scalar, value of the likelihood or posterior kernel.
% - info [integer] 4×1 vector, informations resolution of the model and evaluation of the likelihood.
% - exit_flag [integer] scalar, equal to 1 (no issues when evaluating the likelihood) or 0 (not able to evaluate the likelihood).
% - DLIK [double] Vector with score of the likelihood
% - Hess [double] asymptotic hessian matrix.
% - SteadyState [double] steady state level for the endogenous variables
% - trend_coeff [double] Matrix of doubles, coefficients of the deterministic trend in the measurement equation.
% - M_ [struct] Updated M_ structure described in INPUTS section.
% - options_ [struct] Updated options_ structure described in INPUTS section.
% - bayestopt_ [struct] See INPUTS section.
% - dr [structure] Reduced form model.
% Copyright (C) 2017-2023 Dynare Team
%
@ -28,7 +59,7 @@ DLIK = [];
Hess = [];
% Exit with error if analytical_derivation option is used.
if DynareOptions.analytic_derivation
if options_.analytic_derivation
error('The analytic_derivation and conditional_likelihood are not compatible!')
end
@ -43,9 +74,9 @@ end
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
M_ = set_all_parameters(xparam1,estim_params_,M_);
[fval, info, exitflag, Q, H] = check_bounds_and_definiteness_estimation(xparam1, Model, EstimatedParameters, BoundsInfo);
[fval, info, exitflag, Q, H] = check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, BoundsInfo);
if info(1)
return
end
@ -53,7 +84,7 @@ end
iQ_upper_chol = chol(inv(Q));
% Return an error if the interface for measurement errors is used.
if ~isequal(H, zeros(size(H))) || EstimatedParameters.ncn || EstimatedParameters.ncx
if ~isequal(H, zeros(size(H))) || estim_params_.ncn || estim_params_.ncx
error('Option conditional_likelihood does not support declaration of measurement errors. You can specify the measurement errors in the model block directly by adding measurement equations.')
end
@ -62,8 +93,8 @@ end
%------------------------------------------------------------------------------
% 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');
[T, R, SteadyState, info,dr, M_.params] = ...
dynare_resolve(M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, 'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
@ -86,7 +117,7 @@ if info(1)
end
% check endogenous prior restrictions
info = endogenous_prior_restrictions(T, R, Model, DynareOptions, DynareResults);
info = endogenous_prior_restrictions(T, R, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
if info(1)
fval = Inf;
info(4)=info(2);
@ -95,40 +126,40 @@ 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;
% Define the constant vector of the measurement equation.
if ~DynareOptions.noconstant
if DynareOptions.loglinear
constant = log(SteadyState(BayesInfo.mfys));
if ~options_.noconstant
if options_.loglinear
constant = log(SteadyState(bayestopt_.mfys));
else
constant = SteadyState(BayesInfo.mfys);
constant = SteadyState(bayestopt_.mfys);
end
end
% Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend
[trend_addition, trend_coeff] = compute_trend_coefficients(Model, DynareOptions, DynareDataset.vobs, DynareDataset.nobs);
Y = bsxfun(@minus, transpose(DynareDataset.data), constant)-trend_addition;
if bayestopt_.with_trend
[trend_addition, trend_coeff] = compute_trend_coefficients(M_, options_, dataset_.vobs, dataset_.nobs);
Y = bsxfun(@minus, transpose(dataset_.data), constant)-trend_addition;
else
trend_coeff = zeros(DynareDataset.vobs, 1);
if ~DynareOptions.noconstant
Y = bsxfun(@minus, transpose(DynareDataset.data), constant);
trend_coeff = zeros(dataset_.vobs, 1);
if ~options_.noconstant
Y = bsxfun(@minus, transpose(dataset_.data), constant);
else
Y = transpose(DynareDataset.data);
Y = transpose(dataset_.data);
end
end
% Return an error if some observations are missing.
if DatasetInfo.missing.state
if dataset_info.missing.state
error('Option conditional_likelihood is not compatible with missing observations.')
end
% Get the selection matrix (vector of row indices for T and R)
Z = BayesInfo.mf;
Z = bayestopt_.mf;
% Get the number of observed variables.
pp = DynareDataset.vobs;
pp = dataset_.vobs;
% Get the number of variables in the state equations (state variables plus observed variables).
mm = size(T, 1);
@ -157,12 +188,12 @@ llik = zeros(size(Y, 2), 1);
Ytild = U\(L\Y);
Ttild = U\(L\T(Z,:));
for t = 1:DynareOptions.presample
for t = 1:options_.presample
epsilon = Ytild(:,t) - Ttild*S;
S = T*S + R*epsilon;
end
for t=(DynareOptions.presample+1):size(Y, 2)
for t=(options_.presample+1):size(Y, 2)
epsilon = Ytild(:,t) - Ttild*S;
upsilon = iQ_upper_chol*epsilon;
S = T*S + R*epsilon;
@ -177,17 +208,17 @@ likelihood = -sum(llik);
% 5. 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);
if DynareOptions.endogenous_prior==1
[lnpriormom] = endogenous_prior(Y, Pstar, BayesInfo, H);
if options_.endogenous_prior==1
[lnpriormom] = endogenous_prior(Y, Pstar, bayestopt_, H);
fval = (likelihood-lnprior-lnpriormom);
else
fval = (likelihood-lnprior);
end
if DynareOptions.prior_restrictions.status
tmp = feval(DynareOptions.prior_restrictions.routine, Model, DynareResults, DynareOptions, DynareDataset, DatasetInfo);
if options_.prior_restrictions.status
tmp = feval(options_.prior_restrictions.routine, M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_, dataset_, dataset_info);
fval = fval - tmp;
end

503
matlab/dsge_likelihood.m
View File

@ -1,119 +1,41 @@
function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults,derivatives_info)
% Evaluates the posterior kernel of a dsge model using the specified
function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr] = dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,derivatives_info)
% [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,oo_] = dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,oo_,derivatives_info)
% Evaluates the posterior kernel of a DSGE model using the specified
% kalman_algo; the resulting posterior includes the 2*pi constant of the
% likelihood function
%@info:
%! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults},@var{DLIK},@var{AHess}] =} dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults},@var{derivatives_flag})
%! @anchor{dsge_likelihood}
%! @sp 1
%! Evaluates the posterior kernel of a dsge model.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item xparam1
%! Vector of doubles, current values for the estimated parameters.
%! @item DynareDataset
%! Matlab's structure describing the dataset (initialized by dynare, see @ref{dataset_}).
%! @item DynareOptions
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%! @item Model
%! Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
%! @item EstimatedParamemeters
%! Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
%! @item BayesInfo
%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
%! @item DynareResults
%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%! @item derivates_flag
%! Integer scalar, flag for analytical derivatives of the likelihood.
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item fval
%! Double scalar, value of (minus) the likelihood.
%! @item info
%! Double vector, second entry stores penalty, first entry the error code.
%! @table @ @code
%! @item info==0
%! No error.
%! @item info==1
%! The model doesn't determine the current variables uniquely.
%! @item info==2
%! MJDGGES returned an error code.
%! @item info==3
%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
%! @item info==4
%! Blanchard & Kahn conditions are not satisfied: indeterminacy.
%! @item info==5
%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
%! @item info==6
%! The jacobian evaluated at the deterministic steady state is complex.
%! @item info==19
%! The steadystate routine has thrown an exception (inconsistent deep parameters).
%! @item info==20
%! Cannot find the steady state, info(4) contains the sum of square residuals (of the static equations).
%! @item info==21
%! The steady state is complex, info(4) contains the sum of square of imaginary parts of the steady state.
%! @item info==22
%! The steady has NaNs.
%! @item info==23
%! M_.params has been updated in the steadystate routine and has complex valued scalars.
%! @item info==24
%! M_.params has been updated in the steadystate routine and has some NaNs.
%! @item info==26
%! M_.params has been updated in the steadystate routine and has negative/0 values in loglinear model.
%! @item info==30
%! Ergodic variance can't be computed.
%! @item info==41
%! At least one parameter is violating a lower bound condition.
%! @item info==42
%! At least one parameter is violating an upper bound condition.
%! @item info==43
%! The covariance matrix of the structural innovations is not positive definite.
%! @item info==44
%! The covariance matrix of the measurement errors is not positive definite.
%! @item info==45
%! Likelihood is not a number (NaN).
%! @item info==46
%! Likelihood is a complex valued number.
%! @item info==47
%! Posterior kernel is not a number (logged prior density is NaN)
%! @item info==48
%! Posterior kernel is a complex valued number (logged prior density is complex).
%! @end table
%! @item exit_flag
%! Integer scalar, equal to zero if the routine return with a penalty (one otherwise).
%! @item DLIK
%! Vector of doubles, score of the likelihood.
%! @item AHess
%! Matrix of doubles, asymptotic hessian matrix.
%! @item SteadyState
%! Vector of doubles, steady state level for the endogenous variables.
%! @item trend_coeff
%! Matrix of doubles, coefficients of the deterministic trend in the measurement equation.
%! @item Model
%! Matlab's structure describing the model (initialized by dynare, see @ref{M_}).
%! @item DynareOptions
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%! @item BayesInfo
%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
%! @item DynareResults
%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{dynare_estimation_1}, @ref{mode_check}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{lyapunov_solver}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{get_perturbation_params_deriv}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
%! @end deftypefn
%@eod:
%
% INPUTS
% - xparam1 [double] current values for the estimated parameters.
% - dataset_ [structure] dataset after transformations
% - dataset_info [structure] storing informations about the
% sample; not used but required for interface
% - options_ [structure] Matlab's structure describing the current options
% - M_ [structure] Matlab's structure describing the model
% - estim_params_ [structure] characterizing parameters to be estimated
% - bayestopt_ [structure] describing the priors
% - BoundsInfo [structure] containing prior bounds
% - 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
% - derivatives_info [structure] derivative info for identification
%
% OUTPUTS
% - fval [double] scalar, value of the likelihood or posterior kernel.
% - info [integer] 4×1 vector, informations resolution of the model and evaluation of the likelihood.
% - exit_flag [integer] scalar, equal to 1 (no issues when evaluating the likelihood) or 0 (not able to evaluate the likelihood).
% - DLIK [double] Vector with score of the likelihood
% - Hess [double] asymptotic hessian matrix.
% - SteadyState [double] steady state level for the endogenous variables
% - trend_coeff [double] Matrix of doubles, coefficients of the deterministic trend in the measurement equation.
% - M_ [struct] Updated M_ structure described in INPUTS section.
% - options_ [struct] Updated options_ structure described in INPUTS section.
% - bayestopt_ [struct] See INPUTS section.
% - dr [structure] Reduced form model.
%
% This function is called by: dynare_estimation_1, mode_check
% This function calls: dynare_resolve, lyapunov_symm, lyapunov_solver, compute_Pinf_Pstar, kalman_filter_d, missing_observations_kalman_filter_d,
% univariate_kalman_filter_d, kalman_steady_state, get_perturbation_params_deriv, kalman_filter, missing_observations_kalman_filter, univariate_kalman_filter, priordens
% Copyright © 2004-2023 Dynare Team
%
@ -140,7 +62,7 @@ SteadyState = [];
trend_coeff = [];
exit_flag = 1;
info = zeros(4,1);
if DynareOptions.analytic_derivation
if options_.analytic_derivation
DLIK = NaN(1,length(xparam1));
else
DLIK = [];
@ -156,13 +78,13 @@ if ~isempty(xparam1)
end
% Set flag related to analytical derivatives.
analytic_derivation = DynareOptions.analytic_derivation;
analytic_derivation = options_.analytic_derivation;
if analytic_derivation
if DynareOptions.loglinear
if options_.loglinear
error('The analytic_derivation and loglinear options are not compatible')
end
if DynareOptions.endogenous_prior
if options_.endogenous_prior
error('The analytic_derivation and endogenous_prior options are not compatible')
end
end
@ -172,15 +94,15 @@ if nargout==1
end
if analytic_derivation
kron_flag=DynareOptions.analytic_derivation_mode;
kron_flag=options_.analytic_derivation_mode;
end
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
M_ = set_all_parameters(xparam1,estim_params_,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_, estim_params_, BoundsInfo);
if info(1)
return
end
@ -189,31 +111,31 @@ end
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
is_restrict_state_space = true;
if DynareOptions.occbin.likelihood.status
occbin_options = set_occbin_options(DynareOptions, Model);
if options_.occbin.likelihood.status
occbin_options = set_occbin_options(options_);
if occbin_options.opts_simul.restrict_state_space
[T,R,SteadyState,info,DynareResults.dr, Model.params,TTx,RRx,CCx, T0, R0] = ...
occbin.dynare_resolve(Model,DynareOptions,DynareResults,[],'restrict');
[T,R,SteadyState,info,dr, M_.params,TTx,RRx,CCx, T0, R0] = ...
occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,[],'restrict');
else
is_restrict_state_space = false;
oldoo.restrict_var_list = DynareResults.dr.restrict_var_list;
oldoo.restrict_columns = DynareResults.dr.restrict_columns;
DynareResults.dr.restrict_var_list = BayesInfo.smoother_var_list;
DynareResults.dr.restrict_columns = BayesInfo.smoother_restrict_columns;
oldoo.restrict_var_list = dr.restrict_var_list;
oldoo.restrict_columns = dr.restrict_columns;
dr.restrict_var_list = bayestopt_.smoother_var_list;
dr.restrict_columns = bayestopt_.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,DynareResults.dr, Model.params,TTx,RRx,CCx, T0, R0] = ...
occbin.dynare_resolve(Model,DynareOptions,DynareResults);
[T,R,SteadyState,info,M_,dr, M_.params,TTx,RRx,CCx, T0, R0] = ...
occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
DynareResults.dr.restrict_var_list = oldoo.restrict_var_list;
DynareResults.dr.restrict_columns = oldoo.restrict_columns;
dr.restrict_var_list = oldoo.restrict_var_list;
dr.restrict_columns = oldoo.restrict_columns;
end
occbin_.status = true;
occbin_.info= {DynareOptions, DynareResults, Model, occbin_options, TTx, RRx, CCx,T0,R0};
occbin_.info= {options_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, M_, 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,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict');
[T,R,SteadyState,info,dr, M_.params] = dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,'restrict');
occbin_.status = false;
end
@ -248,7 +170,7 @@ if info(1)
end
% check endogenous prior restrictions
info=endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
info=endogenous_prior_restrictions(T,R,M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if info(1)
fval = Inf;
info(4)=info(2);
@ -261,75 +183,75 @@ end
if is_restrict_state_space
%% Define a vector of indices for the observed variables. Is this really usefull?...
BayesInfo.mf = BayesInfo.mf1;
bayestopt_.mf = bayestopt_.mf1;
else
%get location of observed variables and requested smoothed variables in
%decision rules
BayesInfo.mf = BayesInfo.smoother_var_list(BayesInfo.smoother_mf);
bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf);
end
% Define the constant vector of the measurement equation.
if DynareOptions.noconstant
constant = zeros(DynareDataset.vobs,1);
if options_.noconstant
constant = zeros(dataset_.vobs,1);
else
if DynareOptions.loglinear
constant = log(SteadyState(BayesInfo.mfys));
if options_.loglinear
constant = log(SteadyState(bayestopt_.mfys));
else
constant = SteadyState(BayesInfo.mfys);
constant = SteadyState(bayestopt_.mfys);
end
end
% Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend
[trend_addition, trend_coeff]=compute_trend_coefficients(Model,DynareOptions,DynareDataset.vobs,DynareDataset.nobs);
trend = repmat(constant,1,DynareDataset.nobs)+trend_addition;
if bayestopt_.with_trend
[trend_addition, trend_coeff]=compute_trend_coefficients(M_,options_,dataset_.vobs,dataset_.nobs);
trend = repmat(constant,1,dataset_.nobs)+trend_addition;
else
trend_coeff = zeros(DynareDataset.vobs,1);
trend = repmat(constant,1,DynareDataset.nobs);
trend_coeff = zeros(dataset_.vobs,1);
trend = repmat(constant,1,dataset_.nobs);
end
% Get needed informations for kalman filter routines.
start = DynareOptions.presample+1;
Z = BayesInfo.mf; %selector for observed variables
no_missing_data_flag = ~DatasetInfo.missing.state;
start = options_.presample+1;
Z = bayestopt_.mf; %selector for observed variables
no_missing_data_flag = ~dataset_info.missing.state;
mm = length(T); %number of states
pp = DynareDataset.vobs; %number of observables
pp = dataset_.vobs; %number of observables
rr = length(Q); %number of shocks
kalman_tol = DynareOptions.kalman_tol;
diffuse_kalman_tol = DynareOptions.diffuse_kalman_tol;
riccati_tol = DynareOptions.riccati_tol;
Y = transpose(DynareDataset.data)-trend;
kalman_tol = options_.kalman_tol;
diffuse_kalman_tol = options_.diffuse_kalman_tol;
riccati_tol = options_.riccati_tol;
Y = transpose(dataset_.data)-trend;
smpl = size(Y,2);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
kalman_algo = DynareOptions.kalman_algo;
kalman_algo = options_.kalman_algo;
diffuse_periods = 0;
expanded_state_vector_for_univariate_filter=0;
singular_diffuse_filter = 0;
if DynareOptions.heteroskedastic_filter
Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,Model);
if options_.heteroskedastic_filter
Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,M_);
end
switch DynareOptions.lik_init
switch options_.lik_init
case 1% Standard initialization with the steady state of the state equation.
if kalman_algo~=2
% Use standard kalman filter except if the univariate filter is explicitely choosen.
kalman_algo = 1;
end
Pstar=lyapunov_solver(T,R,Q,DynareOptions);
Pstar=lyapunov_solver(T,R,Q,options_);
Pinf = [];
a = zeros(mm,1);
a=set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a=set_Kalman_starting_values(a,M_,dr,options_,bayestopt_);
a_0_given_tm1=T*a; %set state prediction for first Kalman step;
if DynareOptions.occbin.likelihood.status
Z =zeros(length(BayesInfo.mf),size(T,1));
for i = 1:length(BayesInfo.mf)
Z(i,BayesInfo.mf(i))=1;
if options_.occbin.likelihood.status
Z =zeros(length(bayestopt_.mf),size(T,1));
for i = 1:length(bayestopt_.mf)
Z(i,bayestopt_.mf(i))=1;
end
Zflag = 1;
else
@ -340,15 +262,15 @@ switch DynareOptions.lik_init
% Use standard kalman filter except if the univariate filter is explicitely choosen.
kalman_algo = 1;
end
Pstar = DynareOptions.Harvey_scale_factor*eye(mm);
Pstar = options_.Harvey_scale_factor*eye(mm);
Pinf = [];
a = zeros(mm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a = set_Kalman_starting_values(a,M_,dr,options_,bayestopt_);
a_0_given_tm1 = T*a; %set state prediction for first Kalman step;
if DynareOptions.occbin.likelihood.status
Z =zeros(length(BayesInfo.mf),size(T,1));
for i = 1:length(BayesInfo.mf)
Z(i,BayesInfo.mf(i))=1;
if options_.occbin.likelihood.status
Z =zeros(length(bayestopt_.mf),size(T,1));
for i = 1:length(bayestopt_.mf)
Z(i,bayestopt_.mf(i))=1;
end
Zflag = 1;
else
@ -362,13 +284,13 @@ switch DynareOptions.lik_init
error(['The model requires Diffuse filter, but you specified a different Kalman filter. You must set options_.kalman_algo ' ...
'to 0 (default), 3 or 4'])
end
[Pstar,Pinf] = compute_Pinf_Pstar(Z,T,R,Q,DynareOptions.qz_criterium);
Z =zeros(length(BayesInfo.mf),size(T,1));
for i = 1:length(BayesInfo.mf)
Z(i,BayesInfo.mf(i))=1;
[Pstar,Pinf] = compute_Pinf_Pstar(Z,T,R,Q,options_.qz_criterium);
Z =zeros(length(bayestopt_.mf),size(T,1));
for i = 1:length(bayestopt_.mf)
Z(i,bayestopt_.mf(i))=1;
end
Zflag = 1;
if DynareOptions.heteroskedastic_filter
if options_.heteroskedastic_filter
QQ=Qvec;
else
QQ=Q;
@ -377,19 +299,19 @@ switch DynareOptions.lik_init
if (kalman_algo==3)
% Multivariate Diffuse Kalman Filter
a = zeros(mm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a = set_Kalman_starting_values(a,M_,dr,options_,bayestopt_);
a_0_given_tm1 = T*a; %set state prediction for first Kalman step;
Pstar0 = Pstar; % store Pstar
if no_missing_data_flag
[dLIK,dlik,a_0_given_tm1,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, options_.presample, ...
T,R,QQ,H,Z,mm,pp,rr);
else
[dLIK,dlik,a_0_given_tm1,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ...
[dLIK,dlik,a_0_given_tm1,Pstar] = missing_observations_kalman_filter_d(dataset_info.missing.aindex,dataset_info.missing.number_of_observations,dataset_info.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, options_.presample, ...
T,R,QQ,H,Z,mm,pp,rr);
end
diffuse_periods = length(dlik);
@ -428,7 +350,7 @@ switch DynareOptions.lik_init
Pinf = blkdiag(Pinf,zeros(pp));
H1 = zeros(pp,1);
mmm = mm+pp;
if DynareOptions.heteroskedastic_filter
if options_.heteroskedastic_filter
clear QQ
for kv=1:size(Qvec,3)
QQ(:,:,kv) = blkdiag(Qvec(:,:,kv),H);
@ -441,14 +363,14 @@ switch DynareOptions.lik_init
end
a = zeros(mmm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a = set_Kalman_starting_values(a,M_,dr,options_,bayestopt_);
a_0_given_tm1 = T*a;
[dLIK,dlik,a_0_given_tm1,Pstar] = univariate_kalman_filter_d(DatasetInfo.missing.aindex,...
DatasetInfo.missing.number_of_observations,...
DatasetInfo.missing.no_more_missing_observations, ...
[dLIK,dlik,a_0_given_tm1,Pstar] = univariate_kalman_filter_d(dataset_info.missing.aindex,...
dataset_info.missing.number_of_observations,...
dataset_info.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, options_.presample, ...
T,R,QQ,H1,Z,mmm,pp,rr);
diffuse_periods = size(dlik,1);
end
@ -473,17 +395,17 @@ switch DynareOptions.lik_init
catch ME
disp(ME.message)
disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
DynareOptions.lik_init = 1;
Pstar=lyapunov_solver(T,R,Q,DynareOptions);
options_.lik_init = 1;
Pstar=lyapunov_solver(T,R,Q,options_);
end
Pinf = [];
a = zeros(mm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a = set_Kalman_starting_values(a,M_,dr,options_,bayestopt_);
a_0_given_tm1 = T*a;
if DynareOptions.occbin.likelihood.status
Z =zeros(length(BayesInfo.mf),size(T,1));
for i = 1:length(BayesInfo.mf)
Z(i,BayesInfo.mf(i))=1;
if options_.occbin.likelihood.status
Z =zeros(length(bayestopt_.mf),size(T,1));
for i = 1:length(bayestopt_.mf)
Z(i,bayestopt_.mf(i))=1;
end
Zflag = 1;
else
@ -498,22 +420,22 @@ switch DynareOptions.lik_init
indx_unstable = find(sum(abs(V),2)>1e-5);
stable = find(sum(abs(V),2)<1e-5);
nunit = length(eigenv) - nstable;
Pstar = DynareOptions.Harvey_scale_factor*eye(nunit);
Pstar = options_.Harvey_scale_factor*eye(nunit);
if kalman_algo ~= 2
kalman_algo = 1;
end
R_tmp = R(stable, :);
T_tmp = T(stable,stable);
Pstar_tmp=lyapunov_solver(T_tmp,R_tmp,Q,DynareOptions);
Pstar_tmp=lyapunov_solver(T_tmp,R_tmp,Q,options_);
Pstar(stable, stable) = Pstar_tmp;
Pinf = [];
a = zeros(mm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a = set_Kalman_starting_values(a,M_,dr,options_,bayestopt_);
a_0_given_tm1 = T*a;
if DynareOptions.occbin.likelihood.status
Z =zeros(length(BayesInfo.mf),size(T,1));
for i = 1:length(BayesInfo.mf)
Z(i,BayesInfo.mf(i))=1;
if options_.occbin.likelihood.status
Z =zeros(length(bayestopt_.mf),size(T,1));
for i = 1:length(bayestopt_.mf)
Z(i,bayestopt_.mf(i))=1;
end
Zflag = 1;
else
@ -524,10 +446,10 @@ switch DynareOptions.lik_init
end
if analytic_derivation
offset = EstimatedParameters.nvx;
offset = offset+EstimatedParameters.nvn;
offset = offset+EstimatedParameters.ncx;
offset = offset+EstimatedParameters.ncn;
offset = estim_params_.nvx;
offset = offset+estim_params_.nvn;
offset = offset+estim_params_.ncx;
offset = offset+estim_params_.ncn;
no_DLIK = 0;
full_Hess = analytic_derivation==2;
asy_Hess = analytic_derivation==-2;
@ -540,42 +462,42 @@ if analytic_derivation
end
DLIK = [];
AHess = [];
iv = DynareResults.dr.restrict_var_list;
if nargin<10 || isempty(derivatives_info)
[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);
iv = dr.restrict_var_list;
if nargin<13 || isempty(derivatives_info)
[A,B,nou,nou,dr, M_.params] = dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
if ~isempty(estim_params_.var_exo)
indexo=estim_params_.var_exo(:,1);
else
indexo=[];
end
if ~isempty(EstimatedParameters.param_vals)
indparam=EstimatedParameters.param_vals(:,1);
if ~isempty(estim_params_.param_vals)
indparam=estim_params_.param_vals(:,1);
else
indparam=[];
end
old_order = DynareOptions.order;
if DynareOptions.order > 1%not sure whether this check is necessary
DynareOptions.order = 1; fprintf('Reset order to 1 for analytical parameter derivatives.\n');
old_order = options_.order;
if options_.order > 1%not sure whether this check is necessary
options_.order = 1; fprintf('Reset order to 1 for analytical parameter derivatives.\n');
end
old_analytic_derivation_mode = DynareOptions.analytic_derivation_mode;
DynareOptions.analytic_derivation_mode = kron_flag;
old_analytic_derivation_mode = options_.analytic_derivation_mode;
options_.analytic_derivation_mode = kron_flag;
if full_Hess
DERIVS = get_perturbation_params_derivs(Model, DynareOptions, EstimatedParameters, DynareResults, indparam, indexo, [], true);
indD2T = reshape(1:Model.endo_nbr^2, Model.endo_nbr, Model.endo_nbr);
indD2Om = dyn_unvech(1:Model.endo_nbr*(Model.endo_nbr+1)/2);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indparam, indexo, [], true);
indD2T = reshape(1:M_.endo_nbr^2, M_.endo_nbr, M_.endo_nbr);
indD2Om = dyn_unvech(1:M_.endo_nbr*(M_.endo_nbr+1)/2);
D2T = DERIVS.d2KalmanA(indD2T(iv,iv),:);
D2Om = DERIVS.d2Om(dyn_vech(indD2Om(iv,iv)),:);
D2Yss = DERIVS.d2Yss(iv,:,:);
else
DERIVS = get_perturbation_params_derivs(Model, DynareOptions, EstimatedParameters, DynareResults, indparam, indexo, [], false);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indparam, indexo, [], false);
end
DT = zeros(Model.endo_nbr, Model.endo_nbr, size(DERIVS.dghx,3));
DT(:,Model.nstatic+(1:Model.nspred),:) = DERIVS.dghx;
DT = zeros(M_.endo_nbr, M_.endo_nbr, size(DERIVS.dghx,3));
DT(:,M_.nstatic+(1:M_.nspred),:) = DERIVS.dghx;
DT = DT(iv,iv,:);
DOm = DERIVS.dOm(iv,iv,:);
DYss = DERIVS.dYss(iv,:);
DynareOptions.order = old_order; %make sure order is reset (not sure if necessary)
DynareOptions.analytic_derivation_mode = old_analytic_derivation_mode;%make sure analytic_derivation_mode is reset (not sure if necessary)
options_.order = old_order; %make sure order is reset (not sure if necessary)
options_.analytic_derivation_mode = old_analytic_derivation_mode;%make sure analytic_derivation_mode is reset (not sure if necessary)
else
DT = derivatives_info.DT(iv,iv,:);
DOm = derivatives_info.DOm(iv,iv,:);
@ -606,18 +528,18 @@ if analytic_derivation
D2P=sparse(size(D2Om,1),size(D2Om,2)); %zeros([size(T),length(xparam1),length(xparam1)]);
jcount=0;
end
if DynareOptions.lik_init==1
for i=1:EstimatedParameters.nvx
k =EstimatedParameters.var_exo(i,1);
if options_.lik_init==1
for i=1:estim_params_.nvx
k =estim_params_.var_exo(i,1);
DQ(k,k,i) = 2*sqrt(Q(k,k));
dum = lyapunov_symm(T,DOm(:,:,i),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
dum = lyapunov_symm(T,DOm(:,:,i),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
% kk = find(abs(dum) < 1e-12);
% dum(kk) = 0;
DP(:,:,i)=dum;
if full_Hess
for j=1:i
jcount=jcount+1;
dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
% kk = (abs(dum) < 1e-12);
% dum(kk) = 0;
D2P(:,jcount)=dyn_vech(dum);
@ -626,18 +548,18 @@ if analytic_derivation
end
end
end
offset = EstimatedParameters.nvx;
for i=1:EstimatedParameters.nvn
k = EstimatedParameters.var_endo(i,1);
offset = estim_params_.nvx;
for i=1:estim_params_.nvn
k = estim_params_.var_endo(i,1);
DH(k,k,i+offset) = 2*sqrt(H(k,k));
if full_Hess
D2H(k,k,i+offset,i+offset) = 2;
end
end
offset = offset + EstimatedParameters.nvn;
if DynareOptions.lik_init==1
for j=1:EstimatedParameters.np
dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
offset = offset + estim_params_.nvn;
if options_.lik_init==1
for j=1:estim_params_.np
dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
% kk = find(abs(dum) < 1e-12);
% dum(kk) = 0;
DP(:,:,j+offset)=dum;
@ -651,7 +573,7 @@ if analytic_derivation
D2Tij = reshape(D2T(:,jcount),size(T));
D2Omij = dyn_unvech(D2Om(:,jcount));
tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij;
dum = lyapunov_symm(T,tmp,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],DynareOptions.debug);
dum = lyapunov_symm(T,tmp,options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
% dum(abs(dum)<1.e-12) = 0;
D2P(:,jcount) = dyn_vech(dum);
% D2P(:,:,j+offset,i) = dum;
@ -672,15 +594,15 @@ end
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
if DynareOptions.heteroskedastic_filter
if options_.heteroskedastic_filter
Q=Qvec;
end
singularity_has_been_detected = false;
% First test multivariate filter if specified; potentially abort and use univariate filter instead
if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
if no_missing_data_flag && ~DynareOptions.occbin.likelihood.status
if DynareOptions.fast_kalman_filter
if no_missing_data_flag && ~options_.occbin.likelihood.status
if options_.fast_kalman_filter
if diffuse_periods
%kalman_algo==3 requires no diffuse periods (stationary
%observables) as otherwise FE matrix will not be positive
@ -694,24 +616,24 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
[LIK,lik] = kalman_filter_fast(Y,diffuse_periods+1,size(Y,2), ...
a_0_given_tm1,Pstar, ...
kalman_tol, riccati_tol, ...
DynareOptions.presample, ...
options_.presample, ...
T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ...
analytic_deriv_info{:});
else
[LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ...
a_0_given_tm1,Pstar, ...
kalman_tol, riccati_tol, ...
DynareOptions.rescale_prediction_error_covariance, ...
DynareOptions.presample, ...
options_.rescale_prediction_error_covariance, ...
options_.presample, ...
T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ...
analytic_deriv_info{:});
end
else
[LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
[LIK,lik] = missing_observations_kalman_filter(dataset_info.missing.aindex,dataset_info.missing.number_of_observations,dataset_info.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
a_0_given_tm1, Pstar, ...
kalman_tol, DynareOptions.riccati_tol, ...
DynareOptions.rescale_prediction_error_covariance, ...
DynareOptions.presample, ...
kalman_tol, options_.riccati_tol, ...
options_.rescale_prediction_error_covariance, ...
options_.presample, ...
T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, occbin_);
if occbin_.status && isinf(LIK)
fval = Inf;
@ -728,7 +650,7 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
lik=lik1{1};
end
if isinf(LIK)
if DynareOptions.use_univariate_filters_if_singularity_is_detected
if options_.use_univariate_filters_if_singularity_is_detected
singularity_has_been_detected = true;
if kalman_algo == 1
kalman_algo = 2;
@ -743,7 +665,7 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
return
end
else
if DynareOptions.lik_init==3
if options_.lik_init==3
LIK = LIK + dLIK;
if analytic_derivation==0 && nargout>3
if ~singular_diffuse_filter
@ -785,7 +707,7 @@ if (kalman_algo==2) || (kalman_algo==4)
Z = [Z1, eye(pp)];
Zflag=1;
T = blkdiag(T,zeros(pp));
if DynareOptions.heteroskedastic_filter
if options_.heteroskedastic_filter
clear Q
for kv=1:size(Qvec,3)
Q(:,:,kv) = blkdiag(Qvec(:,:,kv),H);
@ -812,11 +734,11 @@ if (kalman_algo==2) || (kalman_algo==4)
if analytic_derivation
analytic_deriv_info{5}=DH;
end
[LIK, lik] = univariate_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
[LIK, lik] = univariate_kalman_filter(dataset_info.missing.aindex,dataset_info.missing.number_of_observations,dataset_info.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
a_0_given_tm1,Pstar, ...
DynareOptions.kalman_tol, ...
DynareOptions.riccati_tol, ...
DynareOptions.presample, ...
options_.kalman_tol, ...
options_.riccati_tol, ...
options_.presample, ...
T,Q,R,H1,Z,mmm,pp,rr,Zflag,diffuse_periods,analytic_deriv_info{:});
if analytic_derivation
LIK1=LIK;
@ -824,7 +746,7 @@ if (kalman_algo==2) || (kalman_algo==4)
lik1=lik;
lik=lik1{1};
end
if DynareOptions.lik_init==3
if options_.lik_init==3
LIK = LIK+dLIK;
if analytic_derivation==0 && nargout>3
lik = [dlik; lik];
@ -835,19 +757,12 @@ end
if analytic_derivation
if no_DLIK==0
DLIK = LIK1{2};
% [DLIK] = score(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
end
if full_Hess
Hess = -LIK1{3};
% [Hess, DLL] = get_Hessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,Z,kalman_tol,riccati_tol);
% Hess0 = getHessian(Y,T,DT,D2T, R*Q*transpose(R),DOm,D2Om,Z,DYss,D2Yss);
end
if asy_Hess
% if ~((kalman_algo==2) || (kalman_algo==4)),
% [Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
% else
Hess = LIK1{3};
% end
end
end
@ -882,10 +797,10 @@ likelihood = LIK;
% ------------------------------------------------------------------------------
if analytic_derivation
if full_Hess
[lnprior, dlnprior, d2lnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
[lnprior, dlnprior, d2lnprior] = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
Hess = Hess - d2lnprior;
else
[lnprior, dlnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
[lnprior, dlnprior] = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
end
if no_DLIK==0
DLIK = DLIK - dlnprior';
@ -896,22 +811,22 @@ if analytic_derivation
Hess = dlik'*dlik;
end
else
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);
end
if DynareOptions.endogenous_prior==1
if DynareOptions.lik_init==2 || DynareOptions.lik_init==3
if options_.endogenous_prior==1
if options_.lik_init==2 || options_.lik_init==3
error('Endogenous prior not supported with non-stationary models')
else
[lnpriormom] = endogenous_prior(Y,DatasetInfo,Pstar,BayesInfo,H);
[lnpriormom] = endogenous_prior(Y,dataset_info,Pstar,bayestopt_,H);
fval = (likelihood-lnprior-lnpriormom);
end
else
fval = (likelihood-lnprior);
end
if DynareOptions.prior_restrictions.status
tmp = feval(DynareOptions.prior_restrictions.routine, Model, DynareResults, DynareOptions, DynareDataset, DatasetInfo);
if options_.prior_restrictions.status
tmp = feval(options_.prior_restrictions.routine, M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_, dataset_, dataset_info);
fval = fval - tmp;
end
@ -931,9 +846,9 @@ if imag(fval)~=0
return
end
if ~DynareOptions.kalman.keep_kalman_algo_if_singularity_is_detected
% Update DynareOptions.kalman_algo.
DynareOptions.kalman_algo = kalman_algo;
if ~options_.kalman.keep_kalman_algo_if_singularity_is_detected
% Update options_.kalman_algo.
options_.kalman_algo = kalman_algo;
end
if analytic_derivation==0 && nargout>3
@ -941,14 +856,14 @@ if analytic_derivation==0 && nargout>3
DLIK=[-lnprior; lik(:)];
end
function a=set_Kalman_starting_values(a,M_,oo_,options_,bayestopt_)
% function a=set_Kalman_starting_values(a,M_,oo_,options_,bayestopt_)
function a=set_Kalman_starting_values(a,M_,dr,options_,bayestopt_)
% function a=set_Kalman_starting_values(a,M_,dr,options_,bayestopt_)
% Sets initial states guess for Kalman filter/smoother based on M_.filter_initial_state
%
% INPUTS
% o a [double] (p*1) vector of states
% o M_ [structure] decribing the model
% o oo_ [structure] storing the results
% o dr [structure] storing the decision rules
% o options_ [structure] describing the options
% o bayestopt_ [structure] describing the priors
%
@ -956,13 +871,13 @@ function a=set_Kalman_starting_values(a,M_,oo_,options_,bayestopt_)
% o a [double] (p*1) vector of set initial states
if isfield(M_,'filter_initial_state') && ~isempty(M_.filter_initial_state)
state_indices=oo_.dr.order_var(oo_.dr.restrict_var_list(bayestopt_.mf0));
state_indices=dr.order_var(dr.restrict_var_list(bayestopt_.mf0));
for ii=1:size(state_indices,1)
if ~isempty(M_.filter_initial_state{state_indices(ii),1})
if options_.loglinear && ~options_.logged_steady_state
a(bayestopt_.mf0(ii)) = log(eval(M_.filter_initial_state{state_indices(ii),2})) - log(oo_.dr.ys(state_indices(ii)));
a(bayestopt_.mf0(ii)) = log(eval(M_.filter_initial_state{state_indices(ii),2})) - log(dr.ys(state_indices(ii)));
elseif ~options_.loglinear && ~options_.logged_steady_state
a(bayestopt_.mf0(ii)) = eval(M_.filter_initial_state{state_indices(ii),2}) - oo_.dr.ys(state_indices(ii));
a(bayestopt_.mf0(ii)) = eval(M_.filter_initial_state{state_indices(ii),2}) - dr.ys(state_indices(ii));
else
error(['The steady state is logged. This should not happen. Please contact the developers'])
end
@ -970,21 +885,21 @@ if isfield(M_,'filter_initial_state') && ~isempty(M_.filter_initial_state)
end
end
function occbin_options = set_occbin_options(DynareOptions, Model)
function occbin_options = set_occbin_options(options_)
% this builds the opts_simul options field needed by occbin.solver
occbin_options.opts_simul = DynareOptions.occbin.simul;
occbin_options.opts_simul.curb_retrench = DynareOptions.occbin.likelihood.curb_retrench;
occbin_options.opts_simul.maxit = DynareOptions.occbin.likelihood.maxit;
occbin_options.opts_simul.periods = DynareOptions.occbin.likelihood.periods;
occbin_options.opts_simul.check_ahead_periods = DynareOptions.occbin.likelihood.check_ahead_periods;
occbin_options.opts_simul.max_check_ahead_periods = DynareOptions.occbin.likelihood.max_check_ahead_periods;
occbin_options.opts_simul.periodic_solution = DynareOptions.occbin.likelihood.periodic_solution;
occbin_options.opts_simul.restrict_state_space = DynareOptions.occbin.likelihood.restrict_state_space;
occbin_options.opts_simul = options_.occbin.simul;
occbin_options.opts_simul.curb_retrench = options_.occbin.likelihood.curb_retrench;
occbin_options.opts_simul.maxit = options_.occbin.likelihood.maxit;
occbin_options.opts_simul.periods = options_.occbin.likelihood.periods;
occbin_options.opts_simul.check_ahead_periods = options_.occbin.likelihood.check_ahead_periods;
occbin_options.opts_simul.max_check_ahead_periods = options_.occbin.likelihood.max_check_ahead_periods;
occbin_options.opts_simul.periodic_solution = options_.occbin.likelihood.periodic_solution;
occbin_options.opts_simul.restrict_state_space = options_.occbin.likelihood.restrict_state_space;
occbin_options.opts_simul.full_output = DynareOptions.occbin.likelihood.full_output;
occbin_options.opts_simul.piecewise_only = DynareOptions.occbin.likelihood.piecewise_only;
if ~isempty(DynareOptions.occbin.smoother.init_binding_indicator)
occbin_options.opts_simul.init_binding_indicator = DynareOptions.occbin.likelihood.init_binding_indicator;
occbin_options.opts_simul.init_regime_history=DynareOptions.occbin.likelihood.init_regime_history;
occbin_options.opts_simul.full_output = options_.occbin.likelihood.full_output;
occbin_options.opts_simul.piecewise_only = options_.occbin.likelihood.piecewise_only;
if ~isempty(options_.occbin.smoother.init_binding_indicator)
occbin_options.opts_simul.init_binding_indicator = options_.occbin.likelihood.init_binding_indicator;
occbin_options.opts_simul.init_regime_history=options_.occbin.likelihood.init_regime_history;
end

70
matlab/dsge_var_likelihood.m
View File

@ -1,17 +1,21 @@
function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_u_tilde,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DataSetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults)
function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_u_tilde,iXX,prior] = dsge_var_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_u_tilde,iXX,prior] = dsge_var_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% Evaluates the posterior kernel of the BVAR-DSGE model.
%
% INPUTS
% o xparam1 [double] Vector of model's parameters.
% o gend [integer] Number of observations (without conditionning observations for the lags).
% o DynareDataset [dseries] object storing the dataset
% o DataSetInfo [structure] storing informations about the sample.
% o DynareOptions [structure] describing the options
% o Model [structure] decribing the model
% o EstimatedParameters [structure] characterizing parameters to be estimated
% o BayesInfo [structure] describing the priors
% o dataset_ [dseries] object storing the dataset
% o dataset_info [structure] storing informations about the sample.
% o options_ [structure] describing the options
% o M_ [structure] decribing the model
% o estim_params_ [structure] characterizing parameters to be estimated
% o bayestopt_ [structure] describing the priors
% o BoundsInfo [structure] containing prior bounds
% o DynareResults [structure] storing the results
% 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
%
% OUTPUTS
% o fval [double] Value of the posterior kernel at xparam1.
@ -45,7 +49,7 @@ function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_
% SPECIAL REQUIREMENTS
% None.
% Copyright © 2006-2021 Dynare Team
% Copyright © 2006-2023Dynare Team
%
% This file is part of Dynare.
%
@ -82,42 +86,42 @@ if ~isempty(xparam1)
xparam1 = xparam1(:);
end
% Initialization of of the index for parameter dsge_prior_weight in Model.params.
dsge_prior_weight_idx = strmatch('dsge_prior_weight', Model.param_names);
% Initialization of of the index for parameter dsge_prior_weight in M_.params.
dsge_prior_weight_idx = strmatch('dsge_prior_weight', M_.param_names);
% Get the number of estimated (dsge) parameters.
nx = EstimatedParameters.nvx + EstimatedParameters.np;
nx = estim_params_.nvx + estim_params_.np;
% Get the number of observed variables in the VAR model.
NumberOfObservedVariables = DynareDataset.vobs;
NumberOfObservedVariables = dataset_.vobs;
% Get the number of observations.
NumberOfObservations = DynareDataset.nobs;
NumberOfObservations = dataset_.nobs;
% Get the number of lags in the VAR model.
NumberOfLags = DynareOptions.dsge_varlag;
NumberOfLags = options_.dsge_varlag;
% Get the number of parameters in the VAR model.
NumberOfParameters = NumberOfObservedVariables*NumberOfLags ;
if ~DynareOptions.noconstant
if ~options_.noconstant
NumberOfParameters = NumberOfParameters + 1;
end
% Get empirical second order moments for the observed variables.
mYY= DataSetInfo.mYY;
mYX= DataSetInfo.mYX;
mXX= DataSetInfo.mXX;
mYY= dataset_info.mYY;
mYX= dataset_info.mYX;
mXX= dataset_info.mXX;
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
M_ = set_all_parameters(xparam1,estim_params_,M_);
[fval,info,exit_flag,Q]=check_bounds_and_definiteness_estimation(xparam1, Model, EstimatedParameters, BoundsInfo);
[fval,info,exit_flag,Q]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, BoundsInfo);
if info(1)
return
end
% Get the weight of the dsge prior.
dsge_prior_weight = Model.params(dsge_prior_weight_idx);
dsge_prior_weight = M_.params(dsge_prior_weight_idx);
% Is the dsge prior proper?
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/NumberOfObservations
@ -125,7 +129,7 @@ if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/NumberOfObse
exit_flag = 0;
info(1) = 51;
info(2)=dsge_prior_weight;
info(3)=(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.nobs;
info(3)=(NumberOfParameters+NumberOfObservedVariables)/dataset_.nobs;
info(4)=abs(NumberOfObservations*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
return
end
@ -136,7 +140,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.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict');
[T,R,SteadyState,info] = dynare_resolve(M_,options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state,'restrict');
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
@ -157,11 +161,11 @@ if info(1)
end
% Define the mean/steady state vector.
if ~DynareOptions.noconstant
if DynareOptions.loglinear
constant = transpose(log(SteadyState(BayesInfo.mfys)));
if ~options_.noconstant
if options_.loglinear
constant = transpose(log(SteadyState(bayestopt_.mfys)));
else
constant = transpose(SteadyState(BayesInfo.mfys));
constant = transpose(SteadyState(bayestopt_.mfys));
end
else
constant = zeros(1,NumberOfObservedVariables);
@ -173,8 +177,8 @@ end
%------------------------------------------------------------------------------
% Compute the theoretical second order moments
tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
mf = BayesInfo.mf1;
tmp0 = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], options_.debug);
mf = bayestopt_.mf1;
% Get the non centered second order moments
TheoreticalAutoCovarianceOfTheObservedVariables = zeros(NumberOfObservedVariables,NumberOfObservedVariables,NumberOfLags+1);
@ -189,7 +193,7 @@ GYX = zeros(NumberOfObservedVariables,NumberOfParameters);
for i=1:NumberOfLags
GYX(:,(i-1)*NumberOfObservedVariables+1:i*NumberOfObservedVariables) = TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1);
end
if ~DynareOptions.noconstant
if ~options_.noconstant
GYX(:,end) = constant';
end
@ -202,7 +206,7 @@ for i = 1:NumberOfLags-1
GXX = GXX + kron(tmp2,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1)');
end
if ~DynareOptions.noconstant
if ~options_.noconstant
% Add one row and one column to GXX
GXX = [GXX , kron(ones(NumberOfLags,1),constant') ; [ kron(ones(1,NumberOfLags),constant) , 1] ];
end
@ -280,7 +284,7 @@ if imag(lik)~=0
end
% Add the (logged) prior density for the dsge-parameters.
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 = (lik-lnprior);
if isnan(fval)

180
matlab/dynare_estimation_1.m
View File

@ -93,7 +93,7 @@ if options_.order > 1
end
end
%% set objective function
%% set objective function
if ~options_.dsge_var
if options_.particle.status
objective_function = str2func('non_linear_dsge_likelihood');
@ -169,12 +169,12 @@ catch % if check fails, provide info on using calibration if present
rethrow(e);
end
%% Run smoother if no estimation or mode-finding are requested
%% 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
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_);
[~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if ismember(info(1),[303,304,306])
fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n')
else
@ -224,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, ~, 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_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
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
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;
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);
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);
end
options_.kalman_algo = kalman_algo0;
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
end
parameter_names = bayestopt_.name;
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_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
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_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
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_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state), nx, nx);
end
options_.kalman_algo = kalman_algo0;
end
end
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,7 +306,7 @@ if options_.mode_check.status && ~options_.mh_posterior_mode_estimation
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0;
mode_check(objective_function,xparam1,hh,options_,M_,estim_params_,bayestopt_,bounds,false,...
dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
options_.analytic_derivation = ana_deriv_old;
end
@ -344,12 +344,12 @@ end
if options_.particle.status && isfield(options_.particle,'posterior_sampler')
if strcmpi(options_.particle.posterior_sampler,'Herbst_Schorfheide')
Herbst_Schorfheide_sampler(objective_function,xparam1,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
Herbst_Schorfheide_sampler(objective_function,xparam1,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state)
elseif strcmpi(options_.particle.posterior_sampler,'DSMH')
DSMH_sampler(objective_function,xparam1,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
DSMH_sampler(objective_function,xparam1,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state)
end
end
if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
% display results table and store parameter estimates and standard errors in results
oo_ = display_estimation_results_table(xparam1, stdh, M_, options_, estim_params_, bayestopt_, oo_, prior_dist_names, 'Posterior', 'posterior');
@ -358,7 +358,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
estim_params_nbr = size(xparam1,1);
if ispd(invhess)
log_det_invhess = log(det(invhess./(stdh*stdh')))+2*sum(log(stdh));
likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
else
oo_.MarginalDensity.LaplaceApproximation = NaN;
@ -369,7 +369,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
end
if options_.dsge_var
[~,~,~,~,~,~,~,oo_.dsge_var.posterior_mode.PHI_tilde,oo_.dsge_var.posterior_mode.SIGMA_u_tilde,oo_.dsge_var.posterior_mode.iXX,oo_.dsge_var.posterior_mode.prior] =...
feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
elseif ~any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
@ -391,7 +391,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if options_.mh_tune_jscale.status
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings')
options_.mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_, M_, objective_function, xparam1, bounds,...
dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_);
dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
bayestopt_.jscale(:) = options_.mh_jscale;
fprintf('mh_jscale has been set equal to %s\n', num2str(options_.mh_jscale));
else
@ -449,9 +449,9 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
elseif options_.load_mh_file && options_.load_results_after_load_mh
%% load fields from previous MCMC run stored in results-file
field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density
'posterior_mean','posterior_hpdinf','posterior_hpdsup','posterior_median','posterior_variance','posterior_std','posterior_deciles','posterior_density',...% fields set by GetPosteriorParametersStatistics
'prior_density',...%fields set by PlotPosteriorDistributions
};
'posterior_mean','posterior_hpdinf','posterior_hpdsup','posterior_median','posterior_variance','posterior_std','posterior_deciles','posterior_density',...% fields set by GetPosteriorParametersStatistics
'prior_density',...%fields set by PlotPosteriorDistributions
};
for field_iter=1:size(field_names,2)
if isfield(oo_load_mh.oo_,field_names{1,field_iter})
oo_.(field_names{1,field_iter})=oo_load_mh.oo_.(field_names{1,field_iter});
@ -479,23 +479,23 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
error('%s: I cannot compute the posterior moments for the endogenous variables!',dispString)
end
if options_.load_mh_file && options_.mh_replic==0 %user wants to recompute results
[MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
if ~info
generic_post_data_file_name={'Posterior2ndOrderMoments','decomposition','PosteriorVarianceDecomposition','correlation','PosteriorCorrelations','conditional decomposition','PosteriorConditionalVarianceDecomposition'};
for ii=1:length(generic_post_data_file_name)
delete_stale_file([MetropolisFolder filesep M_.fname '_' generic_post_data_file_name{1,ii} '*']);
end
% restore compatibility for loading pre-4.6.2 runs where estim_params_ was not saved; see 6e06acc7 and !1944
NumberOfDrawsFiles = length(dir([M_.dname '/metropolis/' M_.fname '_posterior_draws*' ]));
if NumberOfDrawsFiles>0
temp=load([M_.dname '/metropolis/' M_.fname '_posterior_draws1']);
if ~isfield(temp,'estim_params_')
for file_iter=1:NumberOfDrawsFiles
save([M_.dname '/metropolis/' M_.fname '_posterior_draws' num2str(file_iter)],'estim_params_','-append')
end
end
end
end
[MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
if ~info
generic_post_data_file_name={'Posterior2ndOrderMoments','decomposition','PosteriorVarianceDecomposition','correlation','PosteriorCorrelations','conditional decomposition','PosteriorConditionalVarianceDecomposition'};
for ii=1:length(generic_post_data_file_name)
delete_stale_file([MetropolisFolder filesep M_.fname '_' generic_post_data_file_name{1,ii} '*']);
end
% restore compatibility for loading pre-4.6.2 runs where estim_params_ was not saved; see 6e06acc7 and !1944
NumberOfDrawsFiles = length(dir([M_.dname '/metropolis/' M_.fname '_posterior_draws*' ]));
if NumberOfDrawsFiles>0
temp=load([M_.dname '/metropolis/' M_.fname '_posterior_draws1']);
if ~isfield(temp,'estim_params_')
for file_iter=1:NumberOfDrawsFiles
save([M_.dname '/metropolis/' M_.fname '_posterior_draws' num2str(file_iter)],'estim_params_','-append')
end
end
end
end
end
oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
end
@ -527,7 +527,7 @@ 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
|| ~options_.smoother ) && ~options_.partial_information % to be fixed
%% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
end

35
matlab/endogenous_prior_restrictions.m
View File

@ -1,13 +1,16 @@
function [info, info_irf, info_moment, data_irf, data_moment] = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults)
function [info, info_irf, info_moment, data_irf, data_moment] = endogenous_prior_restrictions(T,R,Model,DynareOptions,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
%[info, info_irf, info_moment, data_irf, data_moment] = endogenous_prior_restrictions(T,R,Model,DynareOptions,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% Check for prior (sign) restrictions on irf's and theoretical moments
%
% INPUTS
% T [double] n*n state space matrix
% R [double] n*k matrix of shocks
% Model [structure]
% DynareOptions [structure]
% DynareResults [structure]
% T [double] n*n state space matrix
% R [double] n*k matrix of shocks
% Model [structure]
% DynareOptions [structure]
% 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
% info [double] check if prior restrictions are matched by the
% model and related info
@ -15,7 +18,7 @@ function [info, info_irf, info_moment, data_irf, data_moment] = endogenous_prior
% info_moment [double] array of test checks for all individual moment restrictions
%
% Copyright © 2013-2018 Dynare Team
% Copyright © 2013-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -45,14 +48,13 @@ if ~isempty(endo_prior_restrictions.irf)
data_irf=cell(size(endo_prior_restrictions.irf,1),1);
if DynareOptions.order>1
error('The algorithm for prior (sign) restrictions on irf''s is currently restricted to first-order decision rules')
return
end
varlist = Model.endo_names(DynareResults.dr.order_var);
varlist = Model.endo_names(dr.order_var);
if isempty(T)
[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);
[T,R,SteadyState,infox,dr, Model.params] = dynare_resolve(Model,DynareOptions,dr, endo_steady_state, exo_steady_state, 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);
varlist = varlist(dr.restrict_var_list);
end
% check if endo_prior_restrictions.irf{:,1} variables are in varlist
varlistok=1;
@ -62,8 +64,8 @@ if ~isempty(endo_prior_restrictions.irf)
end
end
if ~varlistok
varlist = Model.endo_names(DynareResults.dr.order_var);
[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);
varlist = Model.endo_names(dr.order_var);
[T,R,SteadyState,infox,dr, Model.params] = dynare_resolve(Model,DynareOptions,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
end
end
NT=1;
@ -105,7 +107,6 @@ end
if ~isempty(endo_prior_restrictions.moment)
if DynareOptions.order>1
error('The algorithm for prior (sign) restrictions on moments is currently restricted to first-order decision rules')
return
end
data_moment=cell(size(endo_prior_restrictions.moment,1),1);
var_list_ = endo_prior_restrictions.moment{1,1};
@ -137,7 +138,7 @@ if ~isempty(endo_prior_restrictions.moment)
end
end
DynareOptions.ar = max(abs(NTmin),NTmax);
[gamma_y,stationary_vars] = th_autocovariances(DynareResults.dr, ivar, Model, DynareOptions,1);
[gamma_y,stationary_vars] = th_autocovariances(dr, ivar, Model, DynareOptions,1);
for t=NTmin:NTmax
RR = gamma_y{abs(t)+1};
if t==0

2
matlab/ep/extended_path_initialization.m
View File

@ -97,7 +97,7 @@ end
% Set seed.
if ep.set_dynare_seed_to_default
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
end
% hybrid correction

7
matlab/evaluate_likelihood.m
View File

@ -1,4 +1,5 @@
function [llik,parameters] = evaluate_likelihood(parameters,M_,estim_params_,oo_,options_,bayestopt_)
% [llik,parameters] = evaluate_likelihood(parameters,M_,estim_params_,oo_,options_,bayestopt_)
% Evaluate the logged likelihood at parameters.
%
% INPUTS
@ -22,7 +23,7 @@ function [llik,parameters] = evaluate_likelihood(parameters,M_,estim_params_,oo_
% [2] This function use persistent variables for the dataset and the description of the missing observations. Consequently, if this function
% is called more than once (by changing the value of parameters) the sample *must not* change.
% Copyright © 2009-2017 Dynare Team
% Copyright © 2009-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -87,9 +88,9 @@ end
if options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter
llik = -occbin.IVF_posterior(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
llik = -occbin.IVF_posterior(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
ldens = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_);
llik = llik - ldens;

5
matlab/evaluate_smoother.m
View File

@ -1,4 +1,5 @@
function [oo_,M_,options_,bayestopt_,Smoothed_variables_declaration_order_deviation_form,initial_date]=evaluate_smoother(parameters,var_list,M_,oo_,options_,bayestopt_,estim_params_)
% [oo_,M_,options_,bayestopt_,Smoothed_variables_declaration_order_deviation_form,initial_date]=evaluate_smoother(parameters,var_list,M_,oo_,options_,bayestopt_,estim_params_)
% Evaluate the smoother at parameters.
%
% INPUTS
@ -38,7 +39,7 @@ function [oo_,M_,options_,bayestopt_,Smoothed_variables_declaration_order_deviat
% [1] This function use persistent variables for the dataset and the description of the missing observations. Consequently, if this function
% is called more than once (by changing the value of parameters) the sample *must not* change.
% Copyright © 2010-2020 Dynare Team
% Copyright © 2010-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -105,7 +106,7 @@ end
if options_.occbin.smoother.status
if options_.occbin.smoother.inversion_filter
[~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(parameters,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_);
[~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov] = occbin.IVF_posterior(parameters,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if ismember(info(1),[303,304,306])
oo_.occbin.smoother.error_flag=1;
else

8
matlab/fjaco.m
View File

@ -10,7 +10,7 @@ function fjac = fjaco(f,x,varargin)
% OUTPUT
% fjac : finite difference Jacobian
%
% Copyright © 2010-2017,2019-2020 Dynare Team
% Copyright © 2010-2017,2019-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -31,7 +31,7 @@ ff=feval(f,x,varargin{:});
tol = eps.^(1/3); %some default value
if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective')
tol= varargin{5}.dynatol.x;
tol= varargin{4}.dynatol.x;
end
h = tol.*max(abs(x),1);
xh1=x+h; xh0=x-h;
@ -57,9 +57,9 @@ feval(f,x,varargin{:});
%Auxiliary functions
function disp_info_error_identification_perturbation(info,j)
% there are errors in the solution algorithm
probl_par = get_the_name(j,varargin{5}.TeX,varargin{3},varargin{2},varargin{5});
probl_par = get_the_name(j,varargin{4}.TeX,varargin{3},varargin{2},varargin{4});
skipline()
message = get_error_message(info,varargin{5});
message = get_error_message(info,varargin{4});
fprintf('Parameter error in numerical two-sided difference method:\n')
fprintf('Cannot solve the model for %s (info = %d, %s)\n', probl_par, info(1), message);
fprintf('Possible solutions:\n')

123
matlab/get_identification_jacobians.m
View File

@ -1,12 +1,11 @@
function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs)
% function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs)
function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% previously getJJ.m in Dynare 4.5
% Sets up the Jacobians needed for identification analysis
% =========================================================================
% INPUTS
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% oo: [structure] storing the reduced-form solution results
% options: [structure] storing the options
% options_ident: [structure] storing the options for identification
% indpmodel: [modparam_nbr by 1] index of estimated parameters in M_.params;
@ -22,6 +21,10 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM
% in the estimated_params block; if estimated_params block
% is not available, then no corr parameters are selected
% indvobs: [obs_nbr by 1] index of observed (VAROBS) variables
% 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
%
@ -150,94 +153,94 @@ obs_nbr = length(indvobs);
d2flag = 0; % do not compute second parameter derivatives
% Get Jacobians (wrt selected params) of steady state, dynamic model derivatives and perturbation solution matrices for all endogenous variables
oo.dr.derivs = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag);
dr.derivs = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag);
[I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files
yy0 = oo.dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order
Yss = oo.dr.ys(oo.dr.order_var); % steady state in DR order
yy0 = dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order
Yss = dr.ys(dr.order_var); % steady state in DR order
if order == 1
[~, g1 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1);
[~, g1 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1);
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
DYNAMIC = [Yss;
vec(g1)]; %add steady state and put rows of g1 in DR order
dDYNAMIC = [oo.dr.derivs.dYss;
reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state
dDYNAMIC = [dr.derivs.dYss;
reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state
REDUCEDFORM = [Yss;
vec(oo.dr.ghx);
dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu))]; %in DR order
vec(dr.ghx);
dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu))]; %in DR order
dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2, totparam_nbr);
for j=1:totparam_nbr
dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j));
dyn_vech(oo.dr.derivs.dOm(:,:,j))];
dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j));
dyn_vech(dr.derivs.dOm(:,:,j))];
end
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
elseif order == 2
[~, g1, g2 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1);
[~, g1, g2 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1);
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
DYNAMIC = [Yss;
vec(g1);
vec(g2)]; %add steady state and put rows of g1 and g2 in DR order
dDYNAMIC = [oo.dr.derivs.dYss;
reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order
reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg2 in DR order
dDYNAMIC = [dr.derivs.dYss;
reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)); %reshape dg1 in DR order
reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3))]; %reshape dg2 in DR order
REDUCEDFORM = [Yss;
vec(oo.dr.ghx);
dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu));
vec(oo.dr.ghxx);
vec(oo.dr.ghxu);
vec(oo.dr.ghuu);
vec(oo.dr.ghs2)]; %in DR order
vec(dr.ghx);
dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu));
vec(dr.ghxx);
vec(dr.ghxu);
vec(dr.ghuu);
vec(dr.ghs2)]; %in DR order
dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2+endo_nbr*nspred^2+endo_nbr*nspred*exo_nbr+endo_nbr*exo_nbr^2+endo_nbr, totparam_nbr);
for j=1:totparam_nbr
dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j));
dyn_vech(oo.dr.derivs.dOm(:,:,j));
vec(oo.dr.derivs.dghxx(:,:,j));
vec(oo.dr.derivs.dghxu(:,:,j));
vec(oo.dr.derivs.dghuu(:,:,j));
vec(oo.dr.derivs.dghs2(:,j))];
dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j));
dyn_vech(dr.derivs.dOm(:,:,j));
vec(dr.derivs.dghxx(:,:,j));
vec(dr.derivs.dghxu(:,:,j));
vec(dr.derivs.dghuu(:,:,j));
vec(dr.derivs.dghs2(:,j))];
end
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
elseif order == 3
[~, g1, g2, g3 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1);
[~, g1, g2, g3 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1);
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
DYNAMIC = [Yss;
vec(g1);
vec(g2);
vec(g3)]; %add steady state and put rows of g1 and g2 in DR order
dDYNAMIC = [oo.dr.derivs.dYss;
reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order
reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3));
reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg3 in DR order
dDYNAMIC = [dr.derivs.dYss;
reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)); %reshape dg1 in DR order
reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3));
reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3))]; %reshape dg3 in DR order
REDUCEDFORM = [Yss;
vec(oo.dr.ghx);
dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu));
vec(oo.dr.ghxx); vec(oo.dr.ghxu); vec(oo.dr.ghuu); vec(oo.dr.ghs2);
vec(oo.dr.ghxxx); vec(oo.dr.ghxxu); vec(oo.dr.ghxuu); vec(oo.dr.ghuuu); vec(oo.dr.ghxss); vec(oo.dr.ghuss)]; %in DR order
vec(dr.ghx);
dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu));
vec(dr.ghxx); vec(dr.ghxu); vec(dr.ghuu); vec(dr.ghs2);
vec(dr.ghxxx); vec(dr.ghxxu); vec(dr.ghxuu); vec(dr.ghuuu); vec(dr.ghxss); vec(dr.ghuss)]; %in DR order
dREDUCEDFORM = zeros(size(REDUCEDFORM,1)-endo_nbr, totparam_nbr);
for j=1:totparam_nbr
dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j));
dyn_vech(oo.dr.derivs.dOm(:,:,j));
vec(oo.dr.derivs.dghxx(:,:,j)); vec(oo.dr.derivs.dghxu(:,:,j)); vec(oo.dr.derivs.dghuu(:,:,j)); vec(oo.dr.derivs.dghs2(:,j))
vec(oo.dr.derivs.dghxxx(:,:,j)); vec(oo.dr.derivs.dghxxu(:,:,j)); vec(oo.dr.derivs.dghxuu(:,:,j)); vec(oo.dr.derivs.dghuuu(:,:,j)); vec(oo.dr.derivs.dghxss(:,:,j)); vec(oo.dr.derivs.dghuss(:,:,j))];
dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j));
dyn_vech(dr.derivs.dOm(:,:,j));
vec(dr.derivs.dghxx(:,:,j)); vec(dr.derivs.dghxu(:,:,j)); vec(dr.derivs.dghuu(:,:,j)); vec(dr.derivs.dghs2(:,j))
vec(dr.derivs.dghxxx(:,:,j)); vec(dr.derivs.dghxxu(:,:,j)); vec(dr.derivs.dghxuu(:,:,j)); vec(dr.derivs.dghuuu(:,:,j)); vec(dr.derivs.dghxss(:,:,j)); vec(dr.derivs.dghuss(:,:,j))];
end
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
end
% Get (pruned) state space representation:
pruned = pruned_state_space_system(M, options, oo.dr, indvobs, nlags, useautocorr, 1);
pruned = pruned_state_space_system(M, options, dr, indvobs, nlags, useautocorr, 1);
MEAN = pruned.E_y;
dMEAN = pruned.dE_y;
%storage for Jacobians used in dsge_likelihood.m for analytical Gradient and Hession of likelihood (only at order=1)
derivatives_info = struct();
if order == 1
dT = zeros(endo_nbr,endo_nbr,totparam_nbr);
dT(:,(nstatic+1):(nstatic+nspred),:) = oo.dr.derivs.dghx;
dT(:,(nstatic+1):(nstatic+nspred),:) = dr.derivs.dghx;
derivatives_info.DT = dT;
derivatives_info.DOm = oo.dr.derivs.dOm;
derivatives_info.DYss = oo.dr.derivs.dYss;
derivatives_info.DOm = dr.derivs.dOm;
derivatives_info.DYss = dr.derivs.dYss;
end
%% Compute dMOMENTS
@ -255,7 +258,7 @@ if ~no_identification_moments
if kronflag == -1
%numerical derivative of autocovariogram
dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1]
dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M, options, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1]
dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables
else
dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr);
@ -312,7 +315,7 @@ if ~no_identification_spectrum
IA = eye(size(pruned.A,1));
if kronflag == -1
%numerical derivative of spectral density
dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=2]
dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M, options, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2]
kk = 0;
for ig = 1:length(freqs)
kk = kk+1;
@ -389,14 +392,14 @@ if ~no_identification_minimal
dMINIMAL = [];
else
% Derive and check minimal state vector of first-order
SYS.A = oo.dr.ghx(pruned.indx,:);
SYS.dA = oo.dr.derivs.dghx(pruned.indx,:,:);
SYS.B = oo.dr.ghu(pruned.indx,:);
SYS.dB = oo.dr.derivs.dghu(pruned.indx,:,:);
SYS.C = oo.dr.ghx(pruned.indy,:);
SYS.dC = oo.dr.derivs.dghx(pruned.indy,:,:);
SYS.D = oo.dr.ghu(pruned.indy,:);
SYS.dD = oo.dr.derivs.dghu(pruned.indy,:,:);
SYS.A = dr.ghx(pruned.indx,:);
SYS.dA = dr.derivs.dghx(pruned.indx,:,:);
SYS.B = dr.ghu(pruned.indx,:);
SYS.dB = dr.derivs.dghu(pruned.indx,:,:);
SYS.C = dr.ghx(pruned.indy,:);
SYS.dC = dr.derivs.dghx(pruned.indy,:,:);
SYS.D = dr.ghu(pruned.indy,:);
SYS.dD = dr.derivs.dghu(pruned.indy,:,:);
[CheckCO,minnx,SYS] = get_minimal_state_representation(SYS,1);
if CheckCO == 0
@ -415,7 +418,7 @@ if ~no_identification_minimal
dminB = reshape(dminB,size(dminB,1)*size(dminB,2),size(dminB,3));
dminC = reshape(dminC,size(dminC,1)*size(dminC,2),size(dminC,3));
dminD = reshape(dminD,size(dminD,1)*size(dminD,2),size(dminD,3));
dvechSig = reshape(oo.dr.derivs.dSigma_e,exo_nbr*exo_nbr,totparam_nbr);
dvechSig = reshape(dr.derivs.dSigma_e,exo_nbr*exo_nbr,totparam_nbr);
indvechSig= find(tril(ones(exo_nbr,exo_nbr)));
dvechSig = dvechSig(indvechSig,:);
Inx = eye(minnx);

76
matlab/get_perturbation_params_derivs.m
View File

@ -1,5 +1,5 @@
function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag)
% DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag)
function DERIVS = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag)
% DERIVS = get_perturbation_params_derivs(M, options, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag)
% previously getH.m in dynare 4.5
% -------------------------------------------------------------------------
% Computes derivatives (with respect to parameters) of
@ -16,17 +16,24 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i
%
% =========================================================================
% INPUTS
% M: [structure] storing the model information
% options: [structure] storing the options
% estim_params: [structure] storing the estimation information
% oo: [structure] storing the solution results
% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M.params;
% corresponds to model parameters (no stderr and no corr) in estimated_params block
% indpstderr: [stderrparam_nbr by 1] index of selected (estimated) standard errors,
% i.e. for all exogenous variables where 'stderr' is given in the estimated_params block
% indpcorr: [corrparam_nbr by 2] matrix of selected (estimated) correlations,
% i.e. for all exogenous variables where 'corr' is given in the estimated_params block
% d2flag: [boolean] flag to compute second-order parameter derivatives of steady state and first-order Kalman transition matrices
% M: [structure] storing the model information
% options: [structure] storing the options
% estim_params: [structure] storing the estimation information
% dr [structure] Reduced form model.
% endo_steady_state [vector] steady state value for endogenous variables
% exo_steady_state [vector] steady state value for exogenous variables
% exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M.params;
% corresponds to model parameters (no stderr and no corr)
% in estimated_params block
% indpstderr: [stderrparam_nbr by 1] index of selected (estimated) standard errors,
% i.e. for all exogenous variables where 'stderr' is given
% in the estimated_params block
% indpcorr: [corrparam_nbr by 2] matrix of selected (estimated) correlations,
% i.e. for all exogenous variables where 'corr' is given in
% the estimated_params block
% d2flag: [boolean] flag to compute second-order parameter derivatives of steady state
% and first-order Kalman transition matrices
% -------------------------------------------------------------------------
% OUTPUTS
% DERIVS: Structure with the following fields:
@ -154,26 +161,25 @@ end
qz_criterium = options.qz_criterium;
threads_BC = options.threads.kronecker.sparse_hessian_times_B_kronecker_C;
% Get fields from oo
exo_steady_state = oo.exo_steady_state;
ghx = oo.dr.ghx;
ghu = oo.dr.ghu;
% Get fields from dr
ghx = dr.ghx;
ghu = dr.ghu;
if order > 1
ghxx = oo.dr.ghxx;
ghxu = oo.dr.ghxu;
ghuu = oo.dr.ghuu;
ghs2 = oo.dr.ghs2;
ghxx = dr.ghxx;
ghxu = dr.ghxu;
ghuu = dr.ghuu;
ghs2 = dr.ghs2;
end
if order > 2
ghxxx = oo.dr.ghxxx;
ghxxu = oo.dr.ghxxu;
ghxuu = oo.dr.ghxuu;
ghuuu = oo.dr.ghuuu;
ghxss = oo.dr.ghxss;
ghuss = oo.dr.ghuss;
ghxxx = dr.ghxxx;
ghxxu = dr.ghxxu;
ghxuu = dr.ghxuu;
ghuuu = dr.ghuuu;
ghxss = dr.ghxss;
ghuss = dr.ghuss;
end
order_var = oo.dr.order_var;
ys = oo.dr.ys;
order_var = dr.order_var;
ys = dr.ys;
% Some checks
if exo_det_nbr > 0
@ -289,7 +295,7 @@ if analytic_derivation_mode == -1
% - perturbation solution matrices: dghx, dghu, dghxx, dghxu, dghuu, dghs2, dghxxx, dghxxu, dghxuu, dghuuu, dghxss, dghuss
%Parameter Jacobian of covariance matrix and solution matrices (wrt selected stderr, corr and model paramters)
dSig_gh = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params, M, oo, options);
dSig_gh = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Sigma_e = (1:exo_nbr^2);
ind_ghx = ind_Sigma_e(end) + (1:endo_nbr*nspred);
ind_ghu = ind_ghx(end) + (1:endo_nbr*exo_nbr);
@ -342,7 +348,7 @@ if analytic_derivation_mode == -1
end
%Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only)
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options);
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Yss = 1:endo_nbr;
if options.discretionary_policy || options.ramsey_policy
ind_g1 = ind_Yss(end) + (1:M.eq_nbr*yy0ex0_nbr);
@ -368,7 +374,7 @@ if analytic_derivation_mode == -1
% Hessian (wrt paramters) of steady state and first-order solution matrices ghx and Om
% note that hessian_sparse.m (contrary to hessian.m) does not take symmetry into account, but focuses already on unique values
options.order = 1; %make sure only first order
d2Yss_KalmanA_Om = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params, M, oo, options);
d2Yss_KalmanA_Om = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
options.order = order; %make sure to set back
ind_KalmanA = ind_Yss(end) + (1:endo_nbr^2);
DERIVS.d2KalmanA = d2Yss_KalmanA_Om(ind_KalmanA, indp2tottot2); %only unique elements
@ -399,7 +405,7 @@ if analytic_derivation_mode == -2
% computation of d2Yss and d2g1
% note that hessian_sparse does not take symmetry into account, i.e. compare hessian_sparse.m to hessian.m, but focuses already on unique values, which are duplicated below
options.order = 1; %d2flag requires only first order
d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M, oo, options); % d2flag requires only first-order
d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % d2flag requires only first-order
options.order = order; %make sure to set back the order
d2Yss = reshape(full(d2Yss_g1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation
for j=1:endo_nbr
@ -425,7 +431,7 @@ if analytic_derivation_mode == -2
end
%Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only)
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options);
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Yss = 1:endo_nbr;
ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr);
dYss = dYss_g(ind_Yss,:); %in tensor notation, wrt selected model parameters only
@ -592,7 +598,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1)
[II, JJ] = ind2sub([yy0ex0_nbr yy0ex0_nbr], find(g2(j,:))); %g2 is [endo_nbr by yy0ex0_nbr^2]
for i = 1:yy0ex0_nbr
is = find(II==i);
is = is(find(JJ(is)<=yy0_nbr)); %focus only on oo.dr.ys(I) derivatives as exogenous variables are 0 in a stochastic context
is = is(find(JJ(is)<=yy0_nbr)); %focus only on dr.ys(I) derivatives as exogenous variables are 0 in a stochastic context
if ~isempty(is)
tmp_g2 = full(g2(j,find(g2(j,:))));
dg1(j,i,:) = tmp_g2(is)*dyy0(JJ(is),:); %put into tensor notation

27
matlab/get_perturbation_params_derivs_numerical_objective.m
View File

@ -1,5 +1,5 @@
function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options)
%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options)
function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, options, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, options, dr, steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Objective function to compute numerically the Jacobians used for get_perturbation_params_derivs
% =========================================================================
@ -9,8 +9,11 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% outputflag: [string] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% oo: [structure] storing the solution results
% options: [structure] storing the options
% dr [structure] Reduced form model.
% endo_steady_state [vector] steady state value for endogenous variables
% exo_steady_state [vector] steady state value for exogenous variables
% exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% -------------------------------------------------------------------------
%
% OUTPUT
@ -51,7 +54,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);
[oo.dr,info,M.params] = compute_decision_rules(M,options,oo.dr, oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state);
[dr,info,M.params] = compute_decision_rules(M,options,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
Sigma_e = M.Sigma_e;
if info(1) > 0
@ -59,16 +62,16 @@ if info(1) > 0
out = [];
return
else
ys = oo.dr.ys; %steady state of model variables in declaration order
ghx = oo.dr.ghx; ghu = oo.dr.ghu;
ys = dr.ys; %steady state of model variables in declaration order
ghx = dr.ghx; ghu = dr.ghu;
if options.order > 1
ghxx = oo.dr.ghxx; ghxu = oo.dr.ghxu; ghuu = oo.dr.ghuu; ghs2 = oo.dr.ghs2;
ghxx = dr.ghxx; ghxu = dr.ghxu; ghuu = dr.ghuu; ghs2 = dr.ghs2;
end
if options.order > 2
ghxxx = oo.dr.ghxxx; ghxxu = oo.dr.ghxxu; ghxuu = oo.dr.ghxuu; ghxss = oo.dr.ghxss; ghuuu = oo.dr.ghuuu; ghuss = oo.dr.ghuss;
ghxxx = dr.ghxxx; ghxxu = dr.ghxxu; ghxuu = dr.ghxuu; ghxss = dr.ghxss; ghuuu = dr.ghuuu; ghuss = dr.ghuss;
end
end
Yss = ys(oo.dr.order_var); %steady state of model variables in DR order
Yss = ys(dr.order_var); %steady state of model variables in DR order
%% out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]
if strcmp(outputflag,'perturbation_solution')
@ -85,13 +88,13 @@ end
if strcmp(outputflag,'dynamic_model')
[I,~] = find(M.lead_lag_incidence'); %I is used to evaluate dynamic model files
if options.order == 1
[~, g1] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
[~, g1] = feval([M.fname,'.dynamic'], ys(I), exo_steady_state', M.params, ys, 1);
out = [Yss; g1(:)];
elseif options.order == 2
[~, g1, g2] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
[~, g1, g2] = feval([M.fname,'.dynamic'], ys(I), exo_steady_state', M.params, ys, 1);
out = [Yss; g1(:); g2(:)];
elseif options.order == 3
[~, g1, g2, g3] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
[~, g1, g2, g3] = feval([M.fname,'.dynamic'], ys(I), exo_steady_state', M.params, ys, 1);
g3 = unfold_g3(g3, length(ys(I))+M.exo_nbr);
out = [Yss; g1(:); g2(:); g3(:)];
end

2
matlab/global_initialization.m
View File

@ -142,7 +142,7 @@ priordens([],[],[],[],[],[],1);
dyn_first_order_solver();
% Set dynare random generator and seed.
set_dynare_seed('default');
options_=set_dynare_seed_local_options(options_,'default');
% Load user configuration file.
if isfield(options_, 'global_init_file')

5
matlab/gsa/map_calibration.m
View File

@ -1,11 +1,12 @@
function map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResults, EstimatedParameters, BayesInfo)
% map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResults, EstimatedParameters, BayesInfo)
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright © 2014-2016 European Commission
% Copyright © 2014-2018 Dynare Team
% Copyright © 2014-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -107,7 +108,7 @@ if init
[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);
[info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults.dr,DynareResults.steady_state,DynareResults.exo_steady_state,DynareResults.exo_det_steady_state);
if ~isempty(info_irf)
for ij=1:nbr_irf_restrictions
mat_irf{ij}(j,:)=data_irf{ij}(:,2)';

8
matlab/gsa/stab_map_.m
View File

@ -1,7 +1,5 @@
function x0 = stab_map_(OutputDirectoryName,opt_gsa)
%
% function x0 = stab_map_(OutputDirectoryName)
%
% x0 = stab_map_(OutputDirectoryName,opt_gsa)
% Mapping of stability regions in the prior ranges applying
% Monte Carlo filtering techniques.
%
@ -35,7 +33,7 @@ function x0 = stab_map_(OutputDirectoryName,opt_gsa)
% marco.ratto@ec.europa.eu
% Copyright © 2012-2016 European Commission
% Copyright © 2012-2018 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -296,7 +294,7 @@ if fload==0
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_);
info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
infox(j,1)=info(1);
if info(1)
inorestriction(j)=j;

8
matlab/identification_analysis.m
View File

@ -73,7 +73,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% * stoch_simul
% * vec
% =========================================================================
% Copyright © 2008-2021 Dynare Team
% Copyright © 2008-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -138,7 +138,7 @@ error_indicator.identification_spectrum=0;
if info(1) == 0 %no errors in solution
% Compute parameter Jacobians for identification analysis
[MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs);
[MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if isempty(dMINIMAL)
% Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options
error_indicator.identification_minimal = 1;
@ -202,7 +202,7 @@ if info(1) == 0 %no errors in solution
options_ident_local.no_identification_spectrum = 1; %do not recompute dSPECTRUM
options_ident_local.ar = nlags; %store new lag number
options_.ar = nlags; %store new lag number
[~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs);
[~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %new index with non-zero rows
end
@ -296,7 +296,7 @@ if info(1) == 0 %no errors in solution
derivatives_info.no_DLIK = 1;
bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
%note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1.
[fval, info, cost_flag, DLIK, AHess, ys, trend_coeff, M_, options_, bayestopt_, oo_] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_, derivatives_info); %non-used output variables need to be set for octave for some reason
[fval, info, cost_flag, DLIK, AHess, ys, trend_coeff, M_, options_, bayestopt_, oo_.dr] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state. derivatives_info); %non-used output variables need to be set for octave for some reason
%note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block
options_.analytic_derivation = analytic_derivation; %reset option
AHess = -AHess; %take negative of hessian

32
matlab/identification_numerical_objective.m
View File

@ -1,21 +1,23 @@
function out = identification_numerical_objective(params, outputflag, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar, useautocorr, nlags, grid_nbr)
%function out = identification_numerical_objective(params, outputflag, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar, useautocorr, nlags, grid_nbr)
function out = identification_numerical_objective(params, outputflag, estim_params, M, options, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state)
%function out = identification_numerical_objective(params, outputflag, %estim_params, M, options, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Objective function to compute numerically the Jacobians used for identification analysis
% Previously this function was called thet2tau.m
% =========================================================================
% INPUTS
% params: [vector] parameter values at which to evaluate objective function
% stderr parameters come first, corr parameters second, model parameters third
% outputflag: [integer] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% oo: [structure] storing the reduced form solution results
% options: [structure] storing the options
% indpmodel: [vector] Index of model parameters
% indpstderr: [vector] Index of stderr parameters
% indpcorr: [matrix] Index of corr parameters
% indvar: [vector] Index of selected or observed variables
% params: [vector] parameter values at which to evaluate objective function
% stderr parameters come first, corr parameters second, model parameters third
% outputflag: [integer] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% options: [structure] storing the options
% indpmodel: [vector] Index of model parameters
% indpstderr: [vector] Index of stderr parameters
% indvar: [vector] Index of selected or observed variables
% 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
% out: dependent on outputflag
@ -76,9 +78,9 @@ else
end
%% compute Kalman transition matrices and steady state with updated parameters
[oo.dr,info,M.params] = compute_decision_rules(M,options,oo.dr, oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state);
[dr,info,M.params] = compute_decision_rules(M,options,dr, steady_state, exo_steady_state, exo_det_steady_state);
options = rmfield(options,'options_ident');
pruned = pruned_state_space_system(M, options, oo.dr, indvar, nlags, useautocorr, 0);
pruned = pruned_state_space_system(M, options, dr, indvar, nlags, useautocorr, 0);
%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
if outputflag == 1

174
matlab/initial_estimation_checks.m
View File

@ -1,21 +1,21 @@
function DynareResults = initial_estimation_checks(objective_function,xparam1,DynareDataset,DatasetInfo,Model,EstimatedParameters,DynareOptions,BayesInfo,BoundsInfo,DynareResults)
% function DynareResults = initial_estimation_checks(objective_function,xparam1,DynareDataset,DatasetInfo,Model,EstimatedParameters,DynareOptions,BayesInfo,BoundsInfo,DynareResults)
function oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,BoundsInfo,oo_)
% function oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,BoundsInfo,oo_)
% Checks data (complex values, ML evaluation, initial values, BK conditions,..)
%
% INPUTS
% objective_function [function handle] of the objective function
% xparam1: [vector] of parameters to be estimated
% DynareDataset: [dseries] object storing the dataset
% DataSetInfo: [structure] storing informations about the sample.
% Model: [structure] decribing the model
% EstimatedParameters [structure] characterizing parameters to be estimated
% DynareOptions [structure] describing the options
% BayesInfo [structure] describing the priors
% xparam1 [vector] of parameters to be estimated
% dataset_ [dseries] object storing the dataset
% dataset_info [structure] storing informations about the sample.
% M_ [structure] decribing the model
% estim_params_ [structure] characterizing parameters to be estimated
% options_ [structure] describing the options
% bayestopt_ [structure] describing the priors
% BoundsInfo [structure] containing prior bounds
% DynareResults [structure] storing the results
% oo_ [structure] storing the results
%
% OUTPUTS
% DynareResults structure of temporary results
% oo_ [structure] storing the results
%
% SPECIAL REQUIREMENTS
% none
@ -39,81 +39,81 @@ function DynareResults = initial_estimation_checks(objective_function,xparam1,Dy
%get maximum number of simultaneously observed variables for stochastic
%singularity check
maximum_number_non_missing_observations=max(sum(~isnan(DynareDataset.data(2:end,:)),2));
init_number_non_missing_observations=sum(~isnan(DynareDataset.data(1,:)),2);
maximum_number_non_missing_observations=max(sum(~isnan(dataset_.data(2:end,:)),2));
init_number_non_missing_observations=sum(~isnan(dataset_.data(1,:)),2);
if DynareOptions.heteroskedastic_filter
if DynareOptions.order>1
if options_.heteroskedastic_filter
if options_.order>1
error('initial_estimation_checks:: heteroskedastic shocks are only supported with the Kalman filter/smoother')
end
observations_by_period=sum(~isnan(DynareDataset.data),2);
base_shocks = find(diag(Model.Sigma_e));
zero_shocks = ~ismember(1:Model.exo_nbr,base_shocks);
observations_by_period=sum(~isnan(dataset_.data),2);
base_shocks = find(diag(M_.Sigma_e));
zero_shocks = ~ismember(1:M_.exo_nbr,base_shocks);
non_zero_shocks_by_period=repmat(length(base_shocks),size(observations_by_period));
% check periods for which base shocks are scaled to zero
non_zero_shocks_by_period = non_zero_shocks_by_period-sum(Model.heteroskedastic_shocks.Qscale(base_shocks,1:DynareOptions.nobs)==0,1)';
non_zero_shocks_by_period = non_zero_shocks_by_period-sum(M_.heteroskedastic_shocks.Qscale(base_shocks,1:options_.nobs)==0,1)';
% check periods for which base shocks are set to zero
non_zero_shocks_by_period = non_zero_shocks_by_period-sum(Model.heteroskedastic_shocks.Qvalue(base_shocks,1:DynareOptions.nobs)==0,1)';
non_zero_shocks_by_period = non_zero_shocks_by_period-sum(M_.heteroskedastic_shocks.Qvalue(base_shocks,1:options_.nobs)==0,1)';
% check periods for which other shocks are set to a positive number
non_zero_shocks_by_period = non_zero_shocks_by_period+sum(Model.heteroskedastic_shocks.Qvalue(zero_shocks,1:DynareOptions.nobs)>0,1)';
non_zero_shocks_by_period = non_zero_shocks_by_period+sum(M_.heteroskedastic_shocks.Qvalue(zero_shocks,1:options_.nobs)>0,1)';
end
if DynareOptions.order>1
if any(any(isnan(DynareDataset.data)))
if options_.order>1
if any(any(isnan(dataset_.data)))
error('initial_estimation_checks:: particle filtering does not support missing observations')
end
if DynareOptions.prefilter==1
if options_.prefilter==1
error('initial_estimation_checks:: particle filtering does not support the prefilter option')
end
if BayesInfo.with_trend
if bayestopt_.with_trend
error('initial_estimation_checks:: particle filtering does not support trends')
end
if DynareOptions.order>3 && DynareOptions.particle.pruning==1
if options_.order>3 && options_.particle.pruning==1
error('initial_estimation_checks:: the particle filter with order>3 does not support pruning')
end
if DynareOptions.particle.pruning~=DynareOptions.pruning
if options_.particle.pruning~=options_.pruning
warning('initial_estimation_checks:: the pruning settings differ between the particle filter and the one used for IRFs/simulations. Make sure this is intended.')
end
end
if DynareOptions.occbin.likelihood.status || DynareOptions.occbin.smoother.status
if DynareOptions.prefilter
if options_.occbin.likelihood.status || options_.occbin.smoother.status
if options_.prefilter
error('initial_estimation_checks:: Occbin is incompatible with the prefilter option due to the sample mean generally not corresponding to the steady state with an occasionally binding constraint.')
end
if ~DynareOptions.occbin.likelihood.inversion_filter && (DynareOptions.kalman_algo==2 || DynareOptions.kalman_algo==4)
if ~options_.occbin.likelihood.inversion_filter && (options_.kalman_algo==2 || options_.kalman_algo==4)
error('initial_estimation_checks:: Occbin is incompatible with the selected univariate Kalman filter.')
end
if DynareOptions.fast_kalman_filter
if options_.fast_kalman_filter
error('initial_estimation_checks:: Occbin is incompatible with the fast Kalman filter.')
end
end
if (DynareOptions.occbin.likelihood.status && DynareOptions.occbin.likelihood.inversion_filter) || (DynareOptions.occbin.smoother.status && DynareOptions.occbin.smoother.inversion_filter)
err_index= find(diag(Model.Sigma_e)~=0);
if length(err_index)~=length(DynareOptions.varobs)
if (options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter) || (options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter)
err_index= find(diag(M_.Sigma_e)~=0);
if length(err_index)~=length(options_.varobs)
fprintf('initial_estimation_checks:: The IVF requires exactly as many shocks as observables.')
end
var_index=find(any(isnan(DynareDataset.data)));
var_index=find(any(isnan(dataset_.data)));
if ~isempty(var_index)
fprintf('initial_estimation_checks:: The IVF requires exactly as many shocks as observables.\n')
fprintf('initial_estimation_checks:: The data series %s contains NaN, I am therefore dropping shock %s for these time points.\n',...
DynareOptions.varobs{var_index},Model.exo_names{DynareOptions.occbin.likelihood.IVF_shock_observable_mapping(var_index)})
options_.varobs{var_index},M_.exo_names{options_.occbin.likelihood.IVF_shock_observable_mapping(var_index)})
end
end
if DynareOptions.order>1 || (DynareOptions.order==1 && ~ischar(DynareOptions.mode_compute) && DynareOptions.mode_compute==11)
if DynareOptions.order==1 && DynareOptions.mode_compute==11
if options_.order>1 || (options_.order==1 && ~ischar(options_.mode_compute) && options_.mode_compute==11)
if options_.order==1 && options_.mode_compute==11
disp_string='mode_compute=11';
else
disp_string='particle filtering';
end
if Model.H==0
if M_.H==0
error('initial_estimation_checks:: %s requires measurement error on the observables',disp_string)
else
if sum(diag(Model.H)>0)<length(DynareOptions.varobs)
if sum(diag(M_.H)>0)<length(options_.varobs)
error('initial_estimation_checks:: %s requires as many measurement errors as observed variables',disp_string)
else
[~,flag]=chol(Model.H);
[~,flag]=chol(M_.H);
if flag
error('initial_estimation_checks:: the measurement error matrix must be positive definite')
end
@ -121,31 +121,31 @@ if DynareOptions.order>1 || (DynareOptions.order==1 && ~ischar(DynareOptions.mod
end
end
non_zero_ME=length(EstimatedParameters.H_entries_to_check_for_positive_definiteness);
non_zero_ME=length(estim_params_.H_entries_to_check_for_positive_definiteness);
print_init_check_warning=false;
if maximum_number_non_missing_observations>Model.exo_nbr+non_zero_ME
if maximum_number_non_missing_observations>M_.exo_nbr+non_zero_ME
error(['initial_estimation_checks:: Estimation can''t take place because there are less declared shocks than observed variables!'])
end
if init_number_non_missing_observations>Model.exo_nbr+non_zero_ME
if DynareOptions.no_init_estimation_check_first_obs
if init_number_non_missing_observations>M_.exo_nbr+non_zero_ME
if options_.no_init_estimation_check_first_obs
print_init_check_warning=true;
else
error(['initial_estimation_checks:: Estimation can''t take place because there are less declared shocks than observed variables in first period!'])
end
end
if DynareOptions.heteroskedastic_filter
if options_.heteroskedastic_filter
if any(observations_by_period>(non_zero_shocks_by_period+non_zero_ME))
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance: Check heteroskedastic block and shocks calibration!'])
end
else
if maximum_number_non_missing_observations>length(find(diag(Model.Sigma_e)))+non_zero_ME
if maximum_number_non_missing_observations>length(find(diag(M_.Sigma_e)))+non_zero_ME
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!'])
end
end
if init_number_non_missing_observations>length(find(diag(Model.Sigma_e)))+non_zero_ME
if DynareOptions.no_init_estimation_check_first_obs
if init_number_non_missing_observations>length(find(diag(M_.Sigma_e)))+non_zero_ME
if options_.no_init_estimation_check_first_obs
print_init_check_warning=true;
else
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance in first period!'])
@ -157,54 +157,54 @@ if print_init_check_warning
fprintf('ESTIMATION_CHECKS: it may lead to a crash or provide undesired/wrong results later on!\n');
end
if (any(BayesInfo.pshape >0 ) && DynareOptions.mh_replic) && DynareOptions.mh_nblck<1
if (any(bayestopt_.pshape >0 ) && options_.mh_replic) && options_.mh_nblck<1
error(['initial_estimation_checks:: Bayesian estimation cannot be conducted with mh_nblocks=0.'])
end
% check and display warnings if steady-state solves static model (except if diffuse_filter == 1) and if steady-state changes estimated parameters
[DynareResults.steady_state] = check_steady_state_changes_parameters(Model,EstimatedParameters,DynareResults,DynareOptions, [DynareOptions.diffuse_filter==0 DynareOptions.diffuse_filter==0] );
[oo_.steady_state] = check_steady_state_changes_parameters(M_,estim_params_,oo_,options_, [options_.diffuse_filter==0 options_.diffuse_filter==0] );
% check and display warning if negative values of stderr or corr params are outside unit circle for Bayesian estimation
if any(BayesInfo.pshape)
check_prior_stderr_corr(EstimatedParameters,BayesInfo);
if any(bayestopt_.pshape)
check_prior_stderr_corr(estim_params_,bayestopt_);
end
% display warning if some parameters are still NaN
test_for_deep_parameters_calibration(Model);
test_for_deep_parameters_calibration(M_);
[lnprior,~,~,info]= priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
[lnprior,~,~,info]= priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
if any(info)
fprintf('The prior density evaluated at the initial values is Inf for the following parameters: %s\n',BayesInfo.name{info,1})
fprintf('The prior density evaluated at the initial values is Inf for the following parameters: %s\n',bayestopt_.name{info,1})
error('The initial value of the prior is -Inf')
end
if isfield(Model,'filter_initial_state') && ~isempty(Model.filter_initial_state)
state_indices=DynareResults.dr.order_var(DynareResults.dr.restrict_var_list(BayesInfo.mf0));
if isfield(M_,'filter_initial_state') && ~isempty(M_.filter_initial_state)
state_indices=oo_.dr.order_var(oo_.dr.restrict_var_list(bayestopt_.mf0));
for ii=1:size(state_indices,1)
if ~isempty(Model.filter_initial_state{state_indices(ii),1})
if ~isempty(M_.filter_initial_state{state_indices(ii),1})
try
evaluate_expression(Model.filter_initial_state{state_indices(ii),2},Model,DynareResults)
evaluate_expression(M_.filter_initial_state{state_indices(ii),2},M_,oo_)
catch
fprintf('Unable to evaluate the expression\n %s \nfor the filter_initial_state of variable %s\n',Model.filter_initial_state{state_indices(ii),2},Model.endo_names(state_indices(ii),:))
fprintf('Unable to evaluate the expression\n %s \nfor the filter_initial_state of variable %s\n',M_.filter_initial_state{state_indices(ii),2},M_.endo_names(state_indices(ii),:))
end
end
end
end
if ~isreal(DynareDataset.data)
if ~isreal(dataset_.data)
error('initial_estimation_checks: the data contains complex values.')
end
% Evaluate the likelihood.
ana_deriv = DynareOptions.analytic_derivation;
DynareOptions.analytic_derivation=0;
if ~isequal(DynareOptions.mode_compute,11) || ...
(isequal(DynareOptions.mode_compute,11) && isequal(DynareOptions.order,1))
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation=0;
if ~isequal(options_.mode_compute,11) || ...
(isequal(options_.mode_compute,11) && isequal(options_.order,1))
%shut off potentially automatic switch to diffuse filter for the
%purpose of checking stochastic singularity
use_univariate_filters_if_singularity_is_detected_old=DynareOptions.use_univariate_filters_if_singularity_is_detected;
DynareOptions.use_univariate_filters_if_singularity_is_detected=0;
[fval,info] = feval(objective_function,xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults);
use_univariate_filters_if_singularity_is_detected_old=options_.use_univariate_filters_if_singularity_is_detected;
options_.use_univariate_filters_if_singularity_is_detected=0;
[fval,info] = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if info(1)==50
fprintf('\ninitial_estimation_checks:: The forecast error variance in the multivariate Kalman filter became singular.\n')
fprintf('initial_estimation_checks:: This is often a sign of stochastic singularity, but can also sometimes happen by chance\n')
@ -213,30 +213,24 @@ if ~isequal(DynareOptions.mode_compute,11) || ...
error('initial_estimation_checks:: The forecast error variance in the multivariate Kalman filter became singular.')
end
if info(1)==201 || info(1)==202
message=get_error_message(info,DynareOptions);
message=get_error_message(info,options_);
error('initial_estimation_checks:: %s.',message)
end
%reset options
DynareOptions.use_univariate_filters_if_singularity_is_detected=use_univariate_filters_if_singularity_is_detected_old;
options_.use_univariate_filters_if_singularity_is_detected=use_univariate_filters_if_singularity_is_detected_old;
else
info=0;
fval = 0;
end
if DynareOptions.debug
DynareResults.likelihood_at_initial_parameters=fval;
if options_.debug
oo_.likelihood_at_initial_parameters=fval;
end
DynareOptions.analytic_derivation=ana_deriv;
options_.analytic_derivation=ana_deriv;
if DynareOptions.mode_compute==13
if options_.mode_compute==13
error('Options mode_compute=13 is only compatible with quadratic objective functions')
end
% if DynareOptions.mode_compute==5
% if ~strcmp(func2str(objective_function),'dsge_likelihood')
% error('Options mode_compute=5 is not compatible with non linear filters or Dsge-VAR models!')
% end
% end
if isnan(fval)
error('The initial value of the likelihood is NaN')
elseif imag(fval)
@ -244,34 +238,34 @@ elseif imag(fval)
end
if info(1) > 0
if DynareOptions.order>1
[eigenvalues_] = check(Model,DynareOptions, DynareResults);
if any(abs(1-abs(eigenvalues_))<abs(DynareOptions.qz_criterium-1))
if options_.order>1
[eigenvalues_] = check(M_,options_, oo_);
if any(abs(1-abs(eigenvalues_))<abs(options_.qz_criterium-1))
error('Your model has at least one unit root and you are using a nonlinear filter. Please set nonlinear_filter_initialization=3.')
end
else
disp('Error in computing likelihood for initial parameter values')
print_info(info, DynareOptions.noprint, DynareOptions)
print_info(info, options_.noprint, options_)
end
end
if DynareOptions.prefilter==1
if (~DynareOptions.loglinear && any(abs(DynareResults.steady_state(BayesInfo.mfys))>1e-9)) || (DynareOptions.loglinear && any(abs(log(DynareResults.steady_state(BayesInfo.mfys)))>1e-9))
if options_.prefilter==1
if (~options_.loglinear && any(abs(oo_.steady_state(bayestopt_.mfys))>1e-9)) || (options_.loglinear && any(abs(log(oo_.steady_state(bayestopt_.mfys)))>1e-9))
disp(['You are trying to estimate a model with a non zero steady state for the observed endogenous'])
disp(['variables using demeaned data!'])
error('You should change something in your mod file...')
end
end
if ~isequal(DynareOptions.mode_compute,11)
if ~isequal(options_.mode_compute,11)
disp(['Initial value of the log posterior (or likelihood): ' num2str(-fval)]);
end
if DynareOptions.mh_tune_jscale.status && (DynareOptions.mh_tune_jscale.maxiter<DynareOptions.mh_tune_jscale.stepsize)
if options_.mh_tune_jscale.status && (options_.mh_tune_jscale.maxiter<options_.mh_tune_jscale.stepsize)
warning('You specified mh_tune_jscale, but the maximum number of iterations is smaller than the step size. No update will take place.')
end
if ~isempty(DynareOptions.conditional_variance_decomposition) && ~DynareOptions.moments_varendo
if ~isempty(options_.conditional_variance_decomposition) && ~options_.moments_varendo
disp('The conditional_variance_decomposition-option will be ignored. You need to set moments_varendo');
end

30
matlab/kalman/likelihood/missing_observations_kalman_filter.m
View File

@ -1,4 +1,5 @@
function [LIK, lik, a, P] = missing_observations_kalman_filter(data_index,number_of_observations,no_more_missing_observations,Y,start,last,a,P,kalman_tol,riccati_tol,rescale_prediction_error_covariance,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,occbin_)
% [LIK, lik, a, P] = missing_observations_kalman_filter(data_index,number_of_observations,no_more_missing_observations,Y,start,last,a,P,kalman_tol,riccati_tol,rescale_prediction_error_covariance,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,occbin_)
% Computes the likelihood of a state space model in the case with missing observations.
%
% INPUTS
@ -32,7 +33,7 @@ function [LIK, lik, a, P] = missing_observations_kalman_filter(data_index,numbe
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright © 2004-2021 Dynare Team
% Copyright © 2004-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -100,9 +101,12 @@ if occbin_.status
vv = zeros(pp,last);
options_=occbin_.info{1};
oo_=occbin_.info{2};
M_=occbin_.info{3};
occbin_options=occbin_.info{4};
dr=occbin_.info{2};
endo_steady_state=occbin_.info{3};
exo_steady_state=occbin_.info{4};
exo_det_steady_state=occbin_.info{5};
M_=occbin_.info{6};
occbin_options=occbin_.info{7};
opts_regime.regime_history = occbin_options.opts_simul.init_regime;
opts_regime.binding_indicator = occbin_options.opts_simul.init_binding_indicator;
if t>1
@ -113,13 +117,13 @@ 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_, 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};
CC=occbin_.info{7};
T0=occbin_.info{8};
R0=occbin_.info{9};
[~, ~, ~, regimes_] = occbin.check_regimes([], [], [], opts_regime, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
if length(occbin_.info)>7
TT=occbin_.info{8};
RR=occbin_.info{9};
CC=occbin_.info{10};
T0=occbin_.info{11};
R0=occbin_.info{12};
TT = cat(3,TT,T);
RR = cat(3,RR,R);
CC = cat(2,CC,zeros(mm,1));
@ -254,12 +258,12 @@ while notsteady && t<=last
RR01 = cat(3,R,RR(:,:,1));
CC01 = zeros(size(CC,1),2);
CC01(:,2) = CC(:,1);
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a00, a10, P00, P10, data_index0, Z, v0, Y0, H, Qt, T0, R0, TT01, RR01, CC01, regimes_(t:t+1), M_, oo_, options_, occbin_options);
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a00, a10, P00, P10, data_index0, Z, v0, Y0, H, Qt, T0, R0, TT01, RR01, CC01, regimes_(t:t+1), M_, dr, endo_steady_state,exo_steady_state,exo_det_steady_state, options_, occbin_options);
else
if isqvec
Qt = Qvec(:,:,t-1:t+1);
end
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a0(:,t-1),a1(:,t-1:t),P0(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,vv(:,t-1:t),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,oo_,options_,occbin_options);
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a0(:,t-1),a1(:,t-1:t),P0(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,vv(:,t-1:t),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
end
if info
if options_.debug

47
matlab/missing_DiffuseKalmanSmootherH3_Z.m
View File

@ -163,9 +163,12 @@ else
isoccbin = 1;
Qt = repmat(Q,[1 1 3]);
options_=occbin_.info{1};
oo_=occbin_.info{2};
M_=occbin_.info{3};
occbin_options=occbin_.info{4};
dr=occbin_.info{2};
endo_steady_state=occbin_.info{3};
exo_steady_state=occbin_.info{4};
exo_det_steady_state=occbin_.info{5};
M_=occbin_.info{6};
occbin_options=occbin_.info{7};
opts_regime = occbin_options.opts_regime;
% first_period_occbin_update = inf;
if isfield(opts_regime,'regime_history') && ~isempty(opts_regime.regime_history)
@ -174,29 +177,29 @@ 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_, 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
[~, ~, ~, regimes_] = occbin.check_regimes([], [], [], opts_regime, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
if length(occbin_.info)>7
if length(occbin_.info)==9 && options_.smoother_redux
TT=repmat(T,1,1,smpl+1);
RR=repmat(R,1,1,smpl+1);
CC=repmat(zeros(mm,1),1,smpl+1);
T0=occbin_.info{5};
R0=occbin_.info{6};
T0=occbin_.info{8};
R0=occbin_.info{9};
else
TT=occbin_.info{5};
RR=occbin_.info{6};
CC=occbin_.info{7};
TT=occbin_.info{8};
RR=occbin_.info{9};
CC=occbin_.info{10};
% TT = cat(3,TT,T);
% RR = cat(3,RR,R);
% CC = cat(2,CC,zeros(mm,1));
if options_.smoother_redux
my_order_var = oo_.dr.restrict_var_list;
my_order_var = dr.restrict_var_list;
CC = CC(my_order_var,:);
RR = RR(my_order_var,:,:);
TT = TT(my_order_var,my_order_var,:);
T0=occbin_.info{8};
R0=occbin_.info{9};
T0=occbin_.info{11};
R0=occbin_.info{12};
end
if size(TT,3)<(smpl+1)
TT=repmat(T,1,1,smpl+1);
@ -368,12 +371,12 @@ while notsteady && t<smpl
RR01 = cat(3,R,RR(:,:,1));
CC01 = zeros(size(CC,1),2);
CC01(:,2) = CC(:,1);
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a10,P0,P10,data_index0,Z,v0,Fi0,Ki0,Y0,H,Qt,T0,R0,TT01,RR01,CC01,regimes_(t:t+1),M_,oo_,options_,occbin_options,kalman_tol,nk);
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a10,P0,P10,data_index0,Z,v0,Fi0,Ki0,Y0,H,Qt,T0,R0,TT01,RR01,CC01,regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
else
if isqvec
Qt = Qvec(:,:,t-1:t+1);
end
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a(:,t-1),a1(:,t-1:t),P(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,v(:,t-1:t),Fi(:,t-1),Ki(:,:,t-1),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,oo_,options_,occbin_options,kalman_tol,nk);
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a(:,t-1),a1(:,t-1:t),P(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,v(:,t-1:t),Fi(:,t-1),Ki(:,:,t-1),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if ~error_flag
regimes_(t:t+2)=tmp;
@ -498,15 +501,15 @@ while notsteady && t<smpl
opts_simul.SHOCKS = zeros(nk,M_.exo_nbr);
if smoother_redux
tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list)=a(:,t);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var);
tmp(dr.restrict_var_list)=a(:,t);
opts_simul.endo_init = tmp(dr.inv_order_var);
else
opts_simul.endo_init = a(oo_.dr.inv_order_var,t);
opts_simul.endo_init = a(dr.inv_order_var,t);
end
opts_simul.init_regime = []; %regimes_(t);
opts_simul.waitbar=0;
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~, out, ss] = occbin.solver(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
end
for jnk=1:nk
if filter_covariance_flag
@ -518,9 +521,9 @@ while notsteady && t<smpl
if jnk>1
if isoccbin && (t>=first_period_occbin_update || isinf(first_period_occbin_update))
if smoother_redux
aK(jnk,:,t+jnk) = out.piecewise(jnk,oo_.dr.order_var(oo_.dr.restrict_var_list)) - out.ys(oo_.dr.order_var(oo_.dr.restrict_var_list))';
aK(jnk,:,t+jnk) = out.piecewise(jnk,dr.order_var(dr.restrict_var_list)) - out.ys(dr.order_var(dr.restrict_var_list))';
else
aK(jnk,oo_.dr.inv_order_var,t+jnk) = out.piecewise(jnk,:) - out.ys';
aK(jnk,dr.inv_order_var,t+jnk) = out.piecewise(jnk,:) - out.ys';
end
else
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));

23
matlab/non_linear_dsge_likelihood.m
View File

@ -1,17 +1,20 @@
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_)
function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,dr] = non_linear_dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,EstimatedParameters,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,dr] = non_linear_dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,EstimatedParameters,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% 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
% - DatasetInfo [struct] Matlab's structure describing the dataset
% - dataset_ [struct] Matlab's structure containing the dataset
% - dataset_info [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
% - dr [structure] Reduced form model.
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - fval [double] scalar, value of the likelihood or posterior kernel.
@ -23,8 +26,8 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,oo
% - trend_coeff [double] Empty array.
% - 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.
% - bayestopt_ [struct] See INPUTS section.
% - dr [struct] decision rule structure described in INPUTS section.
% Copyright © 2010-2023 Dynare Team
%
@ -80,7 +83,7 @@ end
%------------------------------------------------------------------------------
% 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);
[dr, info, M_.params] = resol(0, M_, options_, dr , endo_steady_state, exo_steady_state, exo_det_steady_state);
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 || ...
@ -104,7 +107,7 @@ bayestopt_.mf = bayestopt_.mf1;
% Get needed informations for kalman filter routines.
start = options_.presample+1;
Y = transpose(DynareDataset.data);
Y = transpose(dataset_.data);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
@ -163,7 +166,7 @@ switch options_.particle.initialization
options_.periods = 5000;
old_DynareOptionspruning = options_.pruning;
options_.pruning = options_.particle.pruning;
y_ = simult(oo_.steady_state, dr,M_,options_);
y_ = simult(endo_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_))) && ~ options_.pruning
fval = Inf;

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

@ -83,7 +83,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
%end
% Set seed for randn().
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables;

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

@ -77,7 +77,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn().
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables+log(det(H));

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

@ -78,7 +78,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot, 2);
Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn().
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
% Initialization of the likelihood.
lik = NaN(T, 1);

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

@ -74,7 +74,7 @@ else
end
if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
end
% Get covariance matrices

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

@ -79,7 +79,7 @@ else
end
if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
end
% Get covariance matrices

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

@ -105,7 +105,7 @@ else
end
if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
end
% Get covariance matrices

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

@ -39,7 +39,7 @@ function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxili
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Set seed for randn().
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
pruning = DynareOptions.particle.pruning;
second_resample = DynareOptions.particle.resampling.status.systematic;
variance_update = true;

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

@ -101,7 +101,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn().
set_dynare_seed('default');
DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
% Initialization of the weights across particles.
weights = ones(1,number_of_particles)/number_of_particles ;

2
matlab/parallel/InitializeComputationalEnvironment.m
View File

@ -64,7 +64,7 @@ isHybridMatlabOctave = isHybridMatlabOctave && ~isoctave;
options_.parallel_info.isHybridMatlabOctave = isHybridMatlabOctave;
if isHybridMatlabOctave
% Reset dynare random generator and seed.
set_dynare_seed('default');
options_=set_dynare_seed_local_options(options_,'default');
end

5
matlab/posterior_sampler.m
View File

@ -104,7 +104,10 @@ localVars = struct('TargetFun', TargetFun, ...
'M_',M_, ...
'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ...
'oo_', oo_,...
'dr', oo_.dr,...
'endo_steady_state', oo_.steady_state,...
'exo_steady_state', oo_.exo_steady_state,...
'exo_det_steady_state', oo_.exo_det_steady_state,...
'varargin',[]);
if strcmp(sampler_options.posterior_sampling_method,'tailored_random_block_metropolis_hastings')

14
matlab/posterior_sampler_core.m
View File

@ -81,7 +81,11 @@ bayestopt_ = myinputs.bayestopt_;
estim_params_ = myinputs.estim_params_;
options_ = myinputs.options_;
M_ = myinputs.M_;
oo_ = myinputs.oo_;
dr = myinputs.dr;
endo_steady_state = myinputs.endo_steady_state;
exo_steady_state=myinputs.exo_steady_state;
exo_det_steady_state=myinputs.exo_det_steady_state;
% Necessary only for remote computing!
if whoiam
% initialize persistent variables in priordens()
@ -126,15 +130,15 @@ for curr_block = fblck:nblck
%
% Set the random number generator type (the seed is useless but needed by the function)
if ~isoctave
set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
else
set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
end
% Set the state of the RNG
set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal);
catch
% If the state set by master is incompatible with the slave, we only reseed
set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
end
mh_recover_flag=0;
if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
@ -191,7 +195,7 @@ for curr_block = fblck:nblck
sampler_options.curr_block = curr_block;
while draw_iter <= nruns(curr_block)
[par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun, last_draw(curr_block,:), last_posterior(curr_block), sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
[par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun, last_draw(curr_block,:), last_posterior(curr_block), sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
x2(draw_index_current_file,:) = par;
last_draw(curr_block,:) = par;

8
matlab/posterior_sampler_initialization.m
View File

@ -176,7 +176,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
end
% Find initial values for the NumberOfBlocks chains...
if NumberOfBlocks > 1 || options_.mh_initialize_from_previous_mcmc.status% Case 1: multiple chains
set_dynare_seed('default');
options_=set_dynare_seed_local_options(options_,'default');
fprintf(fidlog,[' Initial values of the parameters:\n']);
fprintf('%s: Searching for initial values...\n', dispString);
if ~options_.mh_initialize_from_previous_mcmc.status
@ -202,7 +202,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
end
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
ix2(j,new_estimated_parameters) = candidate(new_estimated_parameters);
ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if ~isfinite(ilogpo2(j)) % if returned log-density is
% Inf or Nan (penalized value)
validate = 0;
@ -256,7 +256,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
candidate = transpose(xparam1(:));%
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
ix2 = candidate;
ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
fprintf('%s: Initialization at the posterior mode.\n\n',dispString);
fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']);
for i=1:length(ix2(1,:))
@ -295,7 +295,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
else
for j=1:NumberOfBlocks
% we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
set_dynare_seed(options_.DynareRandomStreams.seed+j);
options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+j);
% record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
[record.InitialSeeds(j).Unifor,record.InitialSeeds(j).Normal] = get_dynare_random_generator_state();
end

2
matlab/prior_sampler.m
View File

@ -98,7 +98,7 @@ while iteration < NumberOfSimulations
M_ = set_all_parameters(params, estim_params_, M_);
[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_);
INFO=endogenous_prior_restrictions(T,R,M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
file_line_number = file_line_number + 1;
iteration = iteration + 1;

2
matlab/save_display_classical_smoother_results.m
View File

@ -36,7 +36,7 @@ gend= dataset_.nobs;
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_);
[~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if ismember(info(1),[303,304,306])
fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n')
else

98
matlab/set_dynare_seed.m
View File

@ -1,8 +1,9 @@
function set_dynare_seed(a,b)
% Set seeds depending on matlab (octave) version. This routine is called in dynare_config and can be called by the
% user in the mod file.
%
% Copyright © 2010-2020 Dynare Team
function set_dynare_seed(varargin)
% set_dynare_seed(varargin)
% Set seeds depending on Matlab (Octave) version. This routine is
% a wrapper for set_dynare_seed_local_options
% Copyright © 2010-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -18,92 +19,11 @@ function set_dynare_seed(a,b)
%
% 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_
if ~nargin
if nargin<1
error('set_dynare_seed:: I need at least one input argument!')
end
matlab_random_streams = ~(isoctave || options_.parallel_info.isHybridMatlabOctave);
if matlab_random_streams% Use new matlab interface.
if nargin==1
if ischar(a) && strcmpi(a,'default')
options_.DynareRandomStreams.algo = 'mt19937ar';
options_.DynareRandomStreams.seed = 0;
s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
reset(RandStream.setGlobalStream(s));
return
end
if ischar(a) && strcmpi(a,'reset')
s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
reset(RandStream.setGlobalStream(s));
return
end
if ~ischar(a) || (ischar(a) && strcmpi(a, 'clock'))
options_.DynareRandomStreams.algo = 'mt19937ar';
if ischar(a)
options_.DynareRandomStreams.seed = rem(floor(now*24*60*60), 2^32);
else
options_.DynareRandomStreams.seed = a;
end
s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
reset(RandStream.setGlobalStream(s));
return
end
error('set_dynare_seed:: something is wrong in the calling sequence!')
elseif nargin==2
if ~ischar(a) || ~( strcmpi(a,'mcg16807') || ...
strcmpi(a,'mlfg6331_64') || ...
strcmpi(a,'mrg32k3a') || ...
strcmpi(a,'mt19937ar') || ...
strcmpi(a,'shr3cong') || ...
strcmpi(a,'swb2712') )
disp('set_dynare_seed:: First argument must be string designing the uniform random number algorithm!')
RandStream.list
skipline()
disp('set_dynare_seed:: Change the first input accordingly...')
skipline()
error(' ')
end
if ~isint(b)
error('set_dynare_seed:: The second input argument must be an integer!')
end
options_.DynareRandomStreams.algo = a;
options_.DynareRandomStreams.seed = b;
s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
reset(RandStream.setGlobalStream(s));
end
else% Use old matlab interface.
if nargin==1
if ischar(a) && strcmpi(a,'default')
if isoctave
options_.DynareRandomStreams.algo = 'state';
else
options_.DynareRandomStreams.algo = 'twister';
end
options_.DynareRandomStreams.seed = 0;
rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
randn('state',options_.DynareRandomStreams.seed);
return
end
if ischar(a) && strcmpi(a,'reset')
rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
randn('state',options_.DynareRandomStreams.seed);
return
end
if (~ischar(a) && isint(a)) || (ischar(a) && strcmpi(a,'clock'))
if ischar(a)
options_.DynareRandomStreams.seed = floor(now*24*60*60);
else
options_.DynareRandomStreams.seed = a;
end
rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
randn('state',options_.DynareRandomStreams.seed);
return
end
error('set_dynare_seed:: Something is wrong in the calling sequence!')
else
error('set_dynare_seed:: Cannot use more than one input argument with your version of Matlab/Octave!')
end
end
options_=set_dynare_seed_local_options(options_,varargin{:});

121
matlab/set_dynare_seed_local_options.m Normal file
View File

@ -0,0 +1,121 @@
function options_=set_dynare_seed_local_options(options_,a,b)
% options_=set_dynare_seed_local_options(options_,a,b)
% Set seeds depending on Matlab (octave) version
% Inputs:
% o options_ options structure
% o a first input argument,
% for single argument input, either
% [number] seed
% [string] 'default'
% 'reset'
% 'clock'
% for two inputs:
% [string] indicating feasible RNG algorithm (default: 'mt19937ar')
% o b second input argument
% [number] seed for algorithm
% specified with first input
% Copyright © 2010-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/>.
if nargin<2
error('set_dynare_seed:: I need at least two input arguments!')
end
matlab_random_streams = ~(isoctave || options_.parallel_info.isHybridMatlabOctave);
if matlab_random_streams% Use new matlab interface.
if nargin==2
if ischar(a) && strcmpi(a,'default')
options_.DynareRandomStreams.algo = 'mt19937ar';
options_.DynareRandomStreams.seed = 0;
s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
reset(RandStream.setGlobalStream(s));
return
end
if ischar(a) && strcmpi(a,'reset')
s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
reset(RandStream.setGlobalStream(s));
return
end
if ~ischar(a) || (ischar(a) && strcmpi(a, 'clock'))
options_.DynareRandomStreams.algo = 'mt19937ar';
if ischar(a)
options_.DynareRandomStreams.seed = rem(floor(now*24*60*60), 2^32);
else
options_.DynareRandomStreams.seed = a;
end
s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
reset(RandStream.setGlobalStream(s));
return
end
error('set_dynare_seed:: something is wrong in the calling sequence!')
elseif nargin==3
if ~ischar(a) || ~( strcmpi(a,'mcg16807') || ...
strcmpi(a,'mlfg6331_64') || ...
strcmpi(a,'mrg32k3a') || ...
strcmpi(a,'mt19937ar') || ...
strcmpi(a,'shr3cong') || ...
strcmpi(a,'swb2712') )
disp('set_dynare_seed:: First argument must be string designing the uniform random number algorithm!')
RandStream.list
skipline()
disp('set_dynare_seed:: Change the first input accordingly...')
skipline()
error(' ')
end
if ~isint(b)
error('set_dynare_seed:: The second input argument must be an integer!')
end
options_.DynareRandomStreams.algo = a;
options_.DynareRandomStreams.seed = b;
s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
reset(RandStream.setGlobalStream(s));
end
else% Use old matlab interface.
if nargin==2
if ischar(a) && strcmpi(a,'default')
if isoctave
options_.DynareRandomStreams.algo = 'state';
else
options_.DynareRandomStreams.algo = 'twister';
end
options_.DynareRandomStreams.seed = 0;
rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
randn('state',options_.DynareRandomStreams.seed);
return
end
if ischar(a) && strcmpi(a,'reset')
rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
randn('state',options_.DynareRandomStreams.seed);
return
end
if (~ischar(a) && isint(a)) || (ischar(a) && strcmpi(a,'clock'))
if ischar(a)
options_.DynareRandomStreams.seed = floor(now*24*60*60);
else
options_.DynareRandomStreams.seed = a;
end
rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
randn('state',options_.DynareRandomStreams.seed);
return
end
error('set_dynare_seed:: Something is wrong in the calling sequence!')
else
error('set_dynare_seed:: Cannot use more than one input argument with your version of Matlab/Octave!')
end
end

2
matlab/simul_static_model.m
View File

@ -53,7 +53,7 @@ if nargin<2 || isempty(innovations)
covariance_matrix_upper_cholesky = chol(covariance_matrix);
% Set seed to its default state.
if options_.bnlms.set_dynare_seed_to_default
set_dynare_seed('default');
options_=set_dynare_seed_local_options(options_,'default');
end
% Simulate structural innovations.
switch options_.bnlms.innovation_distribution

2
tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod
View File

@ -460,7 +460,7 @@ for jj = 1:2
for id_kronflag = 1:length(KRONFLAG)
fprintf('***** %s: d2flag=%d and kronflag=%d *****\n',strparamset, d2flag,KRONFLAG(id_kronflag))
options_.analytic_derivation_mode = KRONFLAG(id_kronflag);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, d2flag);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag);
for id_var = 1:size(lst_dvars,2)
dx = norm( vec(nSYM.(sprintf('%s',lst_dvars{id_var}))) - vec(DERIVS.(sprintf('%s',lst_dvars{id_var}))), Inf);
fprintf('Max absolute deviation for %s: %e\n', lst_dvars{id_var}, dx);

2
tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod
View File

@ -170,7 +170,7 @@ KRONFLAGS = [-1 -2 0 1];
for k = 1:length(KRONFLAGS)
fprintf('KRONFLAG=%d\n',KRONFLAGS(k));
options_.analytic_derivation_mode = KRONFLAGS(k);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, d2flag);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag);
oo_.dr.dg_0 = permute(1/2*DERIVS.dghs2,[1 3 2]);
oo_.dr.dg_1 = cat(2,DERIVS.dghx,DERIVS.dghu) + 3/6*cat(2,DERIVS.dghxss,DERIVS.dghuss);
oo_.dr.dg_2 = 1/2*cat(2,DERIVS.dghxx,DERIVS.dghxu,DERIVS.dghuu);

6
tests/estimation/system_prior_restriction/Gali_2015_prior_restrictions.m
View File

@ -1,12 +1,12 @@
function log_prior_val=Gali_2015_prior_restrictions(M_, oo_, options_, dataset_, dataset_info);
% function prior_val=Gali_2015_prior_restrictions(M_, oo_, options_, dataset_, dataset_info);
function log_prior_val=Gali_2015_prior_restrictions(M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_, dataset_, dataset_info);
% function prior_val=Gali_2015_prior_restrictions(M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_, dataset_, dataset_info);
% Example of a _prior_restrictions-file automatically called during
% estimation
% It imposes a prior of the slope of the New Keynesian Phillips Curve of
% 0.03. As the slope is a composite of other parameters with independent
% priors, a separate function is required to do this.
% Copyright © 2021 Dynare Team
% Copyright © 2021-2023 Dynare Team
%
% This file is part of Dynare.
%