From 3a115d4fcc1e88a37b22c9ae5ea257004c68dfe4 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 27 Sep 2023 08:06:37 +0200 Subject: [PATCH] Remove full oo_ input from likelihood functions --- examples/Gali_2015_prior_restrictions.m | 6 +- matlab/+occbin/IVF_core.m | 17 +- matlab/+occbin/IVF_posterior.m | 25 ++- matlab/+occbin/kalman_update_algo_1.m | 51 ++--- matlab/+occbin/kalman_update_algo_3.m | 59 +++--- matlab/+occbin/match_function.m | 11 +- matlab/DsgeSmoother.m | 3 +- matlab/PosteriorIRF_core1.m | 2 +- matlab/cli/prior.m | 6 +- matlab/dsge_conditional_likelihood_1.m | 14 +- matlab/dsge_likelihood.m | 70 +++---- matlab/dsge_var_likelihood.m | 12 +- matlab/dynare_estimation_1.m | 180 +++++++++--------- matlab/endogenous_prior_restrictions.m | 35 ++-- matlab/evaluate_likelihood.m | 7 +- matlab/evaluate_smoother.m | 5 +- matlab/gsa/map_calibration.m | 5 +- matlab/gsa/stab_map_.m | 8 +- matlab/identification_analysis.m | 4 +- matlab/initial_estimation_checks.m | 2 +- .../missing_observations_kalman_filter.m | 30 +-- matlab/missing_DiffuseKalmanSmootherH3_Z.m | 47 ++--- matlab/non_linear_dsge_likelihood.m | 13 +- matlab/posterior_sampler.m | 5 +- matlab/posterior_sampler_core.m | 8 +- matlab/posterior_sampler_initialization.m | 4 +- matlab/prior_sampler.m | 2 +- .../save_display_classical_smoother_results.m | 2 +- .../Gali_2015_prior_restrictions.m | 6 +- 29 files changed, 341 insertions(+), 298 deletions(-) diff --git a/examples/Gali_2015_prior_restrictions.m b/examples/Gali_2015_prior_restrictions.m index b40c60627..84a6e4fd2 100644 --- a/examples/Gali_2015_prior_restrictions.m +++ b/examples/Gali_2015_prior_restrictions.m @@ -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. % diff --git a/matlab/+occbin/IVF_core.m b/matlab/+occbin/IVF_core.m index 82cacf2a1..30666cd32 100644 --- a/matlab/+occbin/IVF_core.m +++ b/matlab/+occbin/IVF_core.m @@ -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; diff --git a/matlab/+occbin/IVF_posterior.m b/matlab/+occbin/IVF_posterior.m index 726f08256..07dabc48a 100644 --- a/matlab/+occbin/IVF_posterior.m +++ b/matlab/+occbin/IVF_posterior.m @@ -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,Model,DynareOptions,BayesInfo,dr, atT, innov] = IVF_posterior(xparam1,... + dataset_,dataset_info,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state) +% [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOptions,BayesInfo,dr, atT, innov] = IVF_posterior(xparam1,... +% dataset_,dataset_info,DynareOptions,Model,EstimatedParameters,BayesInfo,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 +% - dataset_info [structure] storing informations about the +% sample; not used but required for interface % - 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 % - 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. @@ -25,7 +30,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti % - 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. +% - 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). @@ -76,7 +81,7 @@ err_index=DynareOptions.occbin.likelihood.IVF_shock_observable_mapping; % err_in COVMAT1 = Model.Sigma_e(err_index,err_index); % Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R). -[T,R,SteadyState,info,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, Model.params] = dynare_resolve(Model,DynareOptions,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,15 +107,15 @@ 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(Model,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,DynareOptions,err_index,filtered_errs_init,obs_list,obs); if info(1) fval = Inf; exit_flag = 0; - atT=NaN(size(stateval(:,DynareResults.dr.order_var)')); + atT=NaN(size(stateval(:,dr.order_var)')); innov=NaN(Model.exo_nbr,sample_length); return else - atT = stateval(:,DynareResults.dr.order_var)'; + atT = stateval(:,dr.order_var)'; innov = zeros(Model.exo_nbr,sample_length); innov(diag(Model.Sigma_e)~=0,:)=filtered_errs'; end diff --git a/matlab/+occbin/kalman_update_algo_1.m b/matlab/+occbin/kalman_update_algo_1.m index b0318e7f4..e8d4a487d 100644 --- a/matlab/+occbin/kalman_update_algo_1.m +++ b/matlab/+occbin/kalman_update_algo_1.m @@ -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_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, 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_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, [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_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, [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_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, [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); diff --git a/matlab/+occbin/kalman_update_algo_3.m b/matlab/+occbin/kalman_update_algo_3.m index 3fdfa7bb7..b07a6d227 100644 --- a/matlab/+occbin/kalman_update_algo_3.m +++ b/matlab/+occbin/kalman_update_algo_3.m @@ -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_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state, 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_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state, [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_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state, [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_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state, [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 diff --git a/matlab/+occbin/match_function.m b/matlab/+occbin/match_function.m index c1059470b..f91cc0547 100644 --- a/matlab/+occbin/match_function.m +++ b/matlab/+occbin/match_function.m @@ -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); diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index df5f9f931..470d573e7 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -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 @@ -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 diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m index 91d5b24ce..096c12c78 100644 --- a/matlab/PosteriorIRF_core1.m +++ b/matlab/PosteriorIRF_core1.m @@ -198,7 +198,7 @@ while fpar3 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 % @@ -957,13 +957,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 @@ -971,7 +971,7 @@ if isfield(M_,'filter_initial_state') && ~isempty(M_.filter_initial_state) end end -function occbin_options = set_occbin_options(options_, M_) +function occbin_options = set_occbin_options(options_) % this builds the opts_simul options field needed by occbin.solver occbin_options.opts_simul = options_.occbin.simul; diff --git a/matlab/dsge_var_likelihood.m b/matlab/dsge_var_likelihood.m index 910ef5ad3..bd41bd9fa 100644 --- a/matlab/dsge_var_likelihood.m +++ b/matlab/dsge_var_likelihood.m @@ -1,4 +1,5 @@ -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,DynareDataset,DataSetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,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,DynareDataset,DataSetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state) % Evaluates the posterior kernel of the BVAR-DSGE model. % % INPUTS @@ -11,7 +12,10 @@ function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_ % o EstimatedParameters [structure] characterizing parameters to be estimated % o BayesInfo [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. % @@ -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(Model,DynareOptions, 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) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 1e216266b..368ddb7cf 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -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 diff --git a/matlab/endogenous_prior_restrictions.m b/matlab/endogenous_prior_restrictions.m index 71877a30a..efb52b2c6 100644 --- a/matlab/endogenous_prior_restrictions.m +++ b/matlab/endogenous_prior_restrictions.m @@ -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)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 diff --git a/matlab/evaluate_likelihood.m b/matlab/evaluate_likelihood.m index d25bfc927..25b4735b6 100644 --- a/matlab/evaluate_likelihood.m +++ b/matlab/evaluate_likelihood.m @@ -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; \ No newline at end of file diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m index 6162a3b17..be08936b8 100644 --- a/matlab/evaluate_smoother.m +++ b/matlab/evaluate_smoother.m @@ -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 diff --git a/matlab/gsa/map_calibration.m b/matlab/gsa/map_calibration.m index c1755a474..b115fb29c 100644 --- a/matlab/gsa/map_calibration.m +++ b/matlab/gsa/map_calibration.m @@ -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)'; diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m index ea2a5de62..52d68069f 100644 --- a/matlab/gsa/stab_map_.m +++ b/matlab/gsa/stab_map_.m @@ -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; diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index b01235957..516e7f2ab 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -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. % @@ -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 diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m index 1340f7e5e..1aa5d0bd2 100644 --- a/matlab/initial_estimation_checks.m +++ b/matlab/initial_estimation_checks.m @@ -204,7 +204,7 @@ if ~isequal(DynareOptions.mode_compute,11) || ... %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); + [fval,info] = feval(objective_function,xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults.dr,DynareResults.steady_state,DynareResults.exo_steady_state,DynareResults.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') diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_kalman_filter.m index 8f738a734..1f168fc30 100644 --- a/matlab/kalman/likelihood/missing_observations_kalman_filter.m +++ b/matlab/kalman/likelihood/missing_observations_kalman_filter.m @@ -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 diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m index 1acc1829e..23d57a8c4 100644 --- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m @@ -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 && t1 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)); diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m index ecb861fba..3571f56ee 100644 --- a/matlab/non_linear_dsge_likelihood.m +++ b/matlab/non_linear_dsge_likelihood.m @@ -1,4 +1,4 @@ -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,DynareDataset,DatasetInfo,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. % @@ -11,7 +11,10 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,oo % - 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. @@ -24,7 +27,7 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,M_,options_,bayestopt_,oo % - 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. +% - 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 || ... @@ -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; diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m index 37d9ee19e..eb2b55ea8 100644 --- a/matlab/posterior_sampler.m +++ b/matlab/posterior_sampler.m @@ -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') diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m index 405c89c88..551ba4daf 100644 --- a/matlab/posterior_sampler_core.m +++ b/matlab/posterior_sampler_core.m @@ -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() @@ -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; diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index db70a69da..a557cc652 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -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,:)) diff --git a/matlab/prior_sampler.m b/matlab/prior_sampler.m index 4c8b4cbcc..23ffe5f37 100644 --- a/matlab/prior_sampler.m +++ b/matlab/prior_sampler.m @@ -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; diff --git a/matlab/save_display_classical_smoother_results.m b/matlab/save_display_classical_smoother_results.m index 898c1927f..a0becb01f 100644 --- a/matlab/save_display_classical_smoother_results.m +++ b/matlab/save_display_classical_smoother_results.m @@ -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 diff --git a/tests/estimation/system_prior_restriction/Gali_2015_prior_restrictions.m b/tests/estimation/system_prior_restriction/Gali_2015_prior_restrictions.m index b40c60627..84a6e4fd2 100644 --- a/tests/estimation/system_prior_restriction/Gali_2015_prior_restrictions.m +++ b/tests/estimation/system_prior_restriction/Gali_2015_prior_restrictions.m @@ -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. %