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

Updated online filter routines. 

- Rewrote doc headers
 - Changed function signatures
 - Removed persistent variables
 - Compute the mode of the particle weights
 - Return the covariance matrix of the particles in the last period
 - Various cosmetic changes
remove-submodule^2
Stéphane Adjemian (Charybdis) 2019-07-11 12:45:26 +02:00
parent 3e6841fd38
commit f0dd928b15
2 changed files with 238 additions and 402 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

159
nonlinear-filters/src/online_auxiliary_filter.m
View File

@ -1,27 +1,27 @@
function [xparam,std_param,lb_95,ub_95,median_param] = online_auxiliary_filter(xparam1,DynareDataset,dataset_info,DynareOptions,Model,EstimatedParameters,BayesInfo,bounds,DynareResults)
function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, DynareResults)
% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
%
% INPUTS
% ReducedForm [structure] Matlab's structure describing the reduced form model.
% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors.
% ReducedForm.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% mf [integer] pp*1 vector of indices.
% number_of_particles [integer] scalar.
% - xparam1 [double] n×1 vector, Initial condition for the estimated parameters.
% - DynareDataset [dseries] Sample used for estimation.
% - dataset_info [struct] Description of the sample.
% - DynareOptions [struct] Option values (options_).
% - Model [struct] Description of the model (M_).
% - EstimatedParameters [struct] Description of the estimated parameters (estim_params_).
% - BayesInfo [struct] Prior definition (bayestopt_).
% - DynareResults [struct] Results (oo_).
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters).
% - pmode [double] n×1 vector, mode of the particles at the end of the sample (for the parameters).
% - pmedian [double] n×1 vector, median of the particles at the end of the sample (for the parameters).
% - pstdev [double] n×1 vector, st. dev. of the particles at the end of the sample (for the parameters).
% - p025 [double] n×1 vector, 2.5 percent of the particles are below p025(i) for i=1,…,n.
% - p975 [double] n×1 vector, 97.5 percent of the particles are below p975(i) for i=1,…,n.
% - covariance [double] n×n matrix, covariance of the particles at the end of the sample.
% Copyright (C) 2013-2017 Dynare Team
% Copyright (C) 2013-2019 Dynare Team
%
% This file is part of Dynare.
%
@ -37,21 +37,18 @@ function [xparam,std_param,lb_95,ub_95,median_param] = online_auxiliary_filter(x
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent Y init_flag mf0 mf1 number_of_particles number_of_parameters liu_west_delta liu_west_chol_sigma_bar
persistent start_param sample_size number_of_observed_variables number_of_structural_innovations
% Set seed for randn().
set_dynare_seed('default');
pruning = DynareOptions.particle.pruning;
second_resample = DynareOptions.particle.resampling.status.systematic;
variance_update = 1 ;
variance_update = true;
bounds = prior_bounds(BayesInfo, DynareOptions.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
% initialization of state particles
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
solve_model_for_online_filter(1,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
[~, Model, DynareOptions, DynareResults, ReducedForm] = solve_model_for_online_filter(true, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
% Set persistent variables.
if isempty(init_flag)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_particles = DynareOptions.particle.number_of_particles;
@ -61,9 +58,6 @@ if isempty(init_flag)
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
liu_west_delta = DynareOptions.particle.liu_west_delta;
start_param = xparam1 ;
init_flag = 1;
end
% Get initial conditions for the state particles
StateVectorMean = ReducedForm.StateVectorMean;
@ -80,27 +74,17 @@ b_square = 1-small_a*small_a ;
% Initialization of parameter particles
xparam = zeros(number_of_parameters,number_of_particles);
%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/100 ;
%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/50 ;
%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/20 ;
bounds = prior_bounds(BayesInfo,DynareOptions.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
prior_draw(BayesInfo,DynareOptions.prior_trunc);
for i=1:number_of_particles
info = 1;
while info==1
%candidate = start_param + .001*liu_west_chol_sigma_bar*randn(number_of_parameters,1) ;
%candidate = start_param + bsxfun(@times,stderr,randn(number_of_parameters,1)) ;
info = 12042009;
while info
candidate = prior_draw()';
if all(candidate(:) >= bounds.lb) && all(candidate(:) <= bounds.ub)
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
solve_model_for_online_filter(1,candidate(:),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
if info==0
[info, Model, DynareOptions, DynareResults] = solve_model_for_online_filter(false, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
if ~info
xparam(:,i) = candidate(:);
end
end
end
end
%xparam = bsxfun(@plus,bounds(:,1),bsxfun(@times,(bounds(:,2)-bounds(:,1)),rand(number_of_parameters,number_of_particles))) ;
% Initialization of the weights of particles.
weights = ones(1,number_of_particles)/number_of_particles;
@ -108,6 +92,7 @@ weights = ones(1,number_of_particles)/number_of_particles ;
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables;
mean_xparam = zeros(number_of_parameters,sample_size);
mode_xparam = zeros(number_of_parameters,sample_size);
median_xparam = zeros(number_of_parameters,sample_size);
std_xparam = zeros(number_of_parameters,sample_size);
lb95_xparam = zeros(number_of_parameters,sample_size);
@ -124,7 +109,7 @@ for t=1:sample_size
m_bar = xparam*(weights');
temp = bsxfun(@minus,xparam,m_bar);
sigma_bar = (bsxfun(@times,weights,temp))*(temp');
if variance_update==1
if variance_update
chol_sigma_bar = chol(b_square*sigma_bar)';
end
% Prediction (without shocks)
@ -132,9 +117,9 @@ for t=1:sample_size
tau_tilde = zeros(1,number_of_particles);
for i=1:number_of_particles
% model resolution
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
solve_model_for_online_filter(t,fore_xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
if info==0
[info, Model, DynareOptions, DynareResults, ReducedForm] = ...
solve_model_for_online_filter(false, fore_xparam(:,i), DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
if ~info
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set local state space model (second-order approximation).
@ -148,16 +133,14 @@ for t=1:sample_size
yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
if pruning
yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
[tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,1),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
[tmp, ~] = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2);
else
tmp = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2);
end
PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
% Replace Gaussian density with a Student density with 3 degrees of
% freedom for fat tails.
% Replace Gaussian density with a Student density with 3 degrees of freedom for fat tails.
z = sum(PredictionError.*(ReducedForm.H\PredictionError), 1) ;
tau_tilde(i) = weights(i).*(tpdf(z, 3*ones(size(z)))+1e-99) ;
%tau_tilde(i) = weights(i).*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError),1))) ;
end
end
% particles selection
@ -172,14 +155,14 @@ for t=1:sample_size
% draw in the new distributions
wtilde = zeros(1, number_of_particles);
for i=1:number_of_particles
info=1 ;
while info==1
info = 12042009;
while info
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
% model resolution for new parameters particles
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
solve_model_for_online_filter(t,candidate,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
if info==0
[info, Model, DynareOptions, DynareResults, ReducedForm] = ...
solve_model_for_online_filter(false, candidate, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults) ;
if ~info
xparam(:,i) = candidate ;
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
@ -210,11 +193,13 @@ for t=1:sample_size
end
% normalization
weights = wtilde/sum(wtilde);
if (variance_update==1) && (neff(weights)<DynareOptions.particle.resampling.threshold*sample_size)
variance_update = 0 ;
if variance_update && (neff(weights)<DynareOptions.particle.resampling.threshold*sample_size)
variance_update = false;
end
% final resampling (not advised)
if second_resample==1
if second_resample
[~, idmode] = max(weights);
mode_xparam(:,t) = xparam(:,idmode);
indx = resample(0, weights,DynareOptions.particle);
StateVectors = StateVectors(:,indx) ;
if pruning
@ -233,7 +218,9 @@ for t=1:sample_size
ub95_xparam(i,t) = temp(0.975*number_of_particles);
end
end
if second_resample==0
if second_resample
[~, idmode] = max(weights);
mode_xparam(:,t) = xparam(:,idmode);
mean_xparam(:,t) = xparam*(weights');
mat_var_cov = bsxfun(@minus, xparam,mean_xparam(:,t));
mat_var_cov = mat_var_cov*(bsxfun(@times, mat_var_cov, weights)');
@ -241,21 +228,21 @@ for t=1:sample_size
for i=1:number_of_parameters
temp = sortrows([xparam(i,:)' weights'], 1);
cumulated_weights = cumsum(temp(:,2));
pass1=1 ;
pass2=1 ;
pass3=1 ;
pass1 = false;
pass2 = false;
pass3 = false;
for j=1:number_of_particles
if cumulated_weights(j)>=0.025 && pass1==1
if ~pass1 && cumulated_weights(j)>=0.025
lb95_xparam(i,t) = temp(j,1);
pass1 = 2 ;
pass1 = true;
end
if cumulated_weights(j)>=0.5 && pass2==1
if ~pass2 && cumulated_weights(j)>=0.5
median_xparam(i,t) = temp(j,1);
pass2 = 2 ;
pass2 = true;
end
if cumulated_weights(j)>=0.975 && pass3==1
if ~pass3 && cumulated_weights(j)>=0.975
ub95_xparam(i,t) = temp(j,1);
pass3 = 2 ;
pass3 = true;
end
end
end
@ -267,22 +254,22 @@ for t=1:sample_size
disp([str])
disp('')
end
distrib_param = xparam ;
xparam = mean_xparam(:,sample_size) ;
std_param = std_xparam(:,sample_size) ;
lb_95 = lb95_xparam(:,sample_size) ;
ub_95 = ub95_xparam(:,sample_size) ;
median_param = median_xparam(:,sample_size) ;
pmean = xparam(:,sample_size);
pmode = mode_xparam(:,sample_size);
pstdev = std_xparam(:,sample_size) ;
p025 = lb95_xparam(:,sample_size) ;
p975 = ub95_xparam(:,sample_size) ;
pmedian = median_xparam(:,sample_size) ;
covariance = mat_var_cov;
%% Plot parameters trajectory
TeX = DynareOptions.TeX;
[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters);
nr = ceil(sqrt(number_of_parameters)) ;
nc = floor(sqrt(number_of_parameters));
nbplt = 1 ;
if TeX
fidTeX = fopen([Model.fname '_param_traj.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
@ -290,15 +277,13 @@ if TeX
fprintf(fidTeX,' \n');
end
z = 1:1:sample_size ;
for plt = 1:nbplt,
for plt = 1:nbplt
if TeX
NAMES = [];
TeXNAMES = [];
end
hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories');
for k=1:length(xparam)
for k=1:length(pmean)
subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
if TeX
@ -310,9 +295,11 @@ for plt = 1:nbplt,
TeXNAMES = char(TeXNAMES,texname);
end
end
y = [mean_xparam(k,:)' median_xparam(k,:)' lb95_xparam(k,:)' ub95_xparam(k,:)' xparam(k)*ones(sample_size,1)] ;
plot(z,y);
% Draw the surface for an interval containing 95% of the particles.
shade(1:sample_size, ub95_xparam(k,:)', 1:sample_size, lb95_xparam(k,:)', 'FillType',[1 2], 'LineStyle', 'none', 'Marker', 'none')
hold on
% Draw the mean of particles.
plot(1:sample_size, mean_xparam(k,:), '-k', 'linewidth', 2)
title(name,'interpreter','none')
hold off
axis tight
@ -334,17 +321,17 @@ for plt = 1:nbplt,
end
end
%% Plot Parameter Densities
% Plot Parameter Densities
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation.
for plt = 1:nbplt,
for plt = 1:nbplt
if TeX
NAMES = [];
TeXNAMES = [];
end
hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities');
for k=1:length(xparam)
for k=1:length(pmean)
subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
if TeX
@ -356,8 +343,8 @@ for plt = 1:nbplt,
TeXNAMES = char(TeXNAMES,texname);
end
end
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',number_of_particles,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ...
number_of_particles, optimal_bandwidth, kernel_function);
plot(density(:,1), density(:,2));
hold on

273
nonlinear-filters/src/solve_model_for_online_filter.m
View File

@ -1,107 +1,26 @@
function [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = solve_model_for_online_filter(observation_number,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% solve the dsge model for an particular parameters set.
function [info, Model, DynareOptions, DynareResults, ReducedForm] = ...
solve_model_for_online_filter(setinitialcondition, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults)
%@info:
%! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults}] =} non_linear_dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults})
%! @anchor{dsge_likelihood}
%! @sp 1
%! Evaluates the posterior kernel of a dsge model using a non linear filter.
%! @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_}).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item fval
%! Double scalar, value of (minus) the likelihood.
%! @item exit_flag
%! Integer scalar, equal to zero if the routine return with a penalty (one otherwise).
%! @item ys
%! Vector of doubles, steady state level for the endogenous variables.
%! @item trend_coeffs
%! Matrix of doubles, coefficients of the deterministic trend in the measurement equation.
%! @item info
%! Integer scalar, 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 thrown an exception (inconsistent deep parameters).
%! @item info==20
%! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations).
%! @item info==21
%! The steady state is complex, info(2) 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==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==45
%! Likelihood is a complex valued number.
%! @end table
%! @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{priordens}
%! @end deftypefn
%@eod:
% Solves the dsge model for an particular parameters set.
%
% INPUTS
% - setinitialcondition [logical] return initial condition if true.
% - xparam1 [double] n×1 vector, parameter values.
% - DynareDataset [struct] Dataset for estimation (dataset_).
% - DynareOptions [struct] Dynare options (options_).
% - Model [struct] Model description (M_).
% - EstimatedParameters [struct] Estimated parameters (estim_params_).
% - BayesInfo [struct] Prior definition (bayestopt_).
% - DynareResults [struct] Dynare results (oo_).
%
% OUTPUTS
% - info [integer] scalar, nonzero if any problem occur when computing the reduced form.
% - Model [struct] Model description (M_).
% - DynareOptions [struct] Dynare options (options_).
% - DynareResults [struct] Dynare results (oo_).
% - ReducedForm [struct] Reduced form model.
% Copyright (C) 2013-2017 Dynare Team
% Copyright (C) 2013-2019 Dynare Team
%
% This file is part of Dynare.
%
@ -118,46 +37,25 @@ function [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResu
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
% frederic DOT karame AT univ DASH lemans DOT fr
persistent init_flag restrict_variables_idx state_variables_idx mf0 mf1 number_of_state_variables
%global objective_function_penalty_base
% Declaration of the penalty as a persistent variable.
persistent init_flag
persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1
persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
info = 0;
% Initialization of the returned arguments.
fval = [];
ys = [];
trend_coeff = [];
exit_flag = 1;
% Set the number of observed variables
%nvobs = DynareDataset.info.nvobs;
nvobs = size(DynareDataset.data,1) ;
%------------------------------------------------------------------------------
%----------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
%----------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
%if (DynareOptions.mode_compute~=1) && any(xparam1<BayesInfo.lb)
% k = find(xparam1(:) < BayesInfo.lb);
% fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
% exit_flag = 0;
% info = 41;
% return
%end
% Test if some parameters are smaller than the lower bound of the prior domain.
if any(xparam1<bounds.lb)
info = 41;
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
%if (DynareOptions.mode_compute~=1) && any(xparam1>BayesInfo.ub)
% k = find(xparam1(:)>BayesInfo.ub);
% fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
% exit_flag = 0;
% info = 42;
% return
%end
% Test if some parameters are greater than the upper bound of the prior domain.
if any(xparam1>bounds.ub)
info = 42;
return
end
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
Q = Model.Sigma_e;
@ -173,7 +71,7 @@ if EstimatedParameters.nvn
end
offset = offset+EstimatedParameters.nvn;
else
H = zeros(nvobs);
H = zeros(size(DynareDataset.data, 1));
end
% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
@ -185,18 +83,12 @@ if EstimatedParameters.ncx
Q(k2,k1) = Q(k1,k2);
end
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
% [CholQ,testQ] = chol(Q);
% if testQ
% The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
% a = diag(eig(Q));
% k = find(a < 0);
% if k > 0
% fval = objective_function_penalty_base+sum(-a(k));
% exit_flag = 0;
% info = 43;
% return
% end
% end
[~, testQ] = chol(Q);
if testQ
% The variance-covariance matrix of the structural innovations is not definite positive.
info = 43;
return
end
offset = offset+EstimatedParameters.ncx;
end
@ -210,18 +102,12 @@ if EstimatedParameters.ncn
H(k2,k1) = H(k1,k2);
end
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
% [CholH,testH] = chol(H);
% if testH
% The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
% a = diag(eig(H));
% k = find(a < 0);
% if k > 0
% fval = objective_function_penalty_base+sum(-a(k));
% exit_flag = 0;
% info = 44;
% return
% end
% end
[~, testH] = chol(H);
if testH
% The variance-covariance matrix of the measurement errors is not definite positive.
info = 44;
return
end
offset = offset+EstimatedParameters.ncn;
end
@ -238,55 +124,18 @@ Model.H = H;
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
%disp(info)
warning('off', 'MATLAB:nearlySingularMatrix')
[~, ~, ~, info, Model, DynareOptions, DynareResults] = ...
dynare_resolve(Model, DynareOptions, DynareResults, 'restrict');
warning('on', 'MATLAB:nearlySingularMatrix')
if info(1)~=0
if nargout==5
ReducedForm = 0;
exit_flag = 55;
end
return
end
% Define a vector of indices for the observed variables. Is this really usefull?...
BayesInfo.mf = BayesInfo.mf1;
% Define the deterministic linear trend of the measurement equation.
if DynareOptions.noconstant
constant = zeros(nvobs,1);
else
if DynareOptions.loglinear
constant = log(SteadyState(BayesInfo.mfys));
else
constant = SteadyState(BayesInfo.mfys);
end
end
% Define the deterministic linear trend of the measurement equation.
%if BayesInfo.with_trend
% trend_coeff = zeros(DynareDataset.info.nvobs,1);
% t = DynareOptions.trend_coeffs;
% for i=1:length(t)
% if ~isempty(t{i})
% trend_coeff(i) = evalin('base',t{i});
% end
% end
% trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs];
%else
% trend = repmat(constant,1,DynareDataset.info.ntobs);
%end
% Get needed informations for kalman filter routines.
start = DynareOptions.presample+1;
np = size(T,1);
mf = BayesInfo.mf;
Y = transpose(DynareDataset.data);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
% Get decision rules and transition equations.
dr = DynareResults.dr;
@ -295,15 +144,14 @@ if isempty(init_flag)
mf0 = BayesInfo.mf0;
mf1 = BayesInfo.mf1;
restrict_variables_idx = dr.restrict_var_list;
observed_variables_idx = restrict_variables_idx(mf1);
state_variables_idx = restrict_variables_idx(mf0);
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(Q);
init_flag = 1;
init_flag = true;
end
% Return reduced form model.
if nargout>4
ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
@ -323,9 +171,10 @@ ReducedForm.Q = Q;
ReducedForm.H = H;
ReducedForm.mf0 = mf0;
ReducedForm.mf1 = mf1;
end
% Set initial condition for t=1
if observation_number==1
% Set initial condition
if setinitialcondition
switch DynareOptions.particle.initialization
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
StateVectorMean = ReducedForm.constant(mf0);