diff --git a/nonlinear-filters/src/online_auxiliary_filter.m b/nonlinear-filters/src/online_auxiliary_filter.m
index 8598a1318..18f762ec9 100644
--- a/nonlinear-filters/src/online_auxiliary_filter.m
+++ b/nonlinear-filters/src/online_auxiliary_filter.m
@@ -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,33 +37,27 @@ 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 .
-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') ;
+set_dynare_seed('default');
pruning = DynareOptions.particle.pruning;
-second_resample = DynareOptions.particle.resampling.status.systematic ;
-variance_update = 1 ;
+second_resample = DynareOptions.particle.resampling.status.systematic;
+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;
- number_of_parameters = size(xparam1,1) ;
- Y = DynareDataset.data ;
- sample_size = size(Y,1);
- 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
+mf0 = ReducedForm.mf0;
+mf1 = ReducedForm.mf1;
+number_of_particles = DynareOptions.particle.number_of_particles;
+number_of_parameters = size(xparam1,1);
+Y = DynareDataset.data;
+sample_size = size(Y,1);
+number_of_observed_variables = length(mf1);
+number_of_structural_innovations = length(ReducedForm.Q);
+liu_west_delta = DynareOptions.particle.liu_west_delta;
% Get initial conditions for the state particles
StateVectorMean = ReducedForm.StateVectorMean;
@@ -75,43 +69,34 @@ if pruning
end
% parameters for the Liu & West filter
-small_a = (3*liu_west_delta-1)/(2*liu_west_delta) ;
-b_square = 1-small_a*small_a ;
+small_a = (3*liu_west_delta-1)/(2*liu_west_delta);
+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
+xparam = zeros(number_of_parameters,number_of_particles);
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
- xparam(:,i) = candidate(:) ;
- end
- end
+ [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
-%xparam = bsxfun(@plus,bounds(:,1),bsxfun(@times,(bounds(:,2)-bounds(:,1)),rand(number_of_parameters,number_of_particles))) ;
+end
% Initialization of the weights of particles.
-weights = ones(1,number_of_particles)/number_of_particles ;
+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) ;
-median_xparam = zeros(number_of_parameters,sample_size) ;
-std_xparam = zeros(number_of_parameters,sample_size) ;
-lb95_xparam = zeros(number_of_parameters,sample_size) ;
-ub95_xparam = zeros(number_of_parameters,sample_size) ;
+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);
+ub95_xparam = zeros(number_of_parameters,sample_size);
%% The Online filter
for t=1:sample_size
@@ -121,20 +106,20 @@ for t=1:sample_size
fprintf('\nSubsample with only the first observation.\n\n', int2str(t))
end
% Moments of parameters particles distribution
- m_bar = xparam*(weights') ;
- temp = bsxfun(@minus,xparam,m_bar) ;
- sigma_bar = (bsxfun(@times,weights,temp))*(temp') ;
- if variance_update==1
- chol_sigma_bar = chol(b_square*sigma_bar)' ;
+ m_bar = xparam*(weights');
+ temp = bsxfun(@minus,xparam,m_bar);
+ sigma_bar = (bsxfun(@times,weights,temp))*(temp');
+ if variance_update
+ chol_sigma_bar = chol(b_square*sigma_bar)';
end
% Prediction (without shocks)
- fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ;
- tau_tilde = zeros(1,number_of_particles) ;
+ fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam);
+ 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,38 +133,36 @@ 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);
+ 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.
- 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
+ PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
+ % 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) ;
+ end
end
% particles selection
- tau_tilde = tau_tilde/sum(tau_tilde) ;
- indx = resample(0,tau_tilde',DynareOptions.particle);
- StateVectors = StateVectors(:,indx) ;
- xparam = fore_xparam(:,indx);
+ tau_tilde = tau_tilde/sum(tau_tilde);
+ indx = resample(0, tau_tilde', DynareOptions.particle);
+ StateVectors = StateVectors(:,indx);
+ xparam = fore_xparam(:,indx);
if pruning
- StateVectors_ = StateVectors_(:,indx) ;
+ StateVectors_ = StateVectors_(:,indx);
end
- w_stage1 = weights(indx)./tau_tilde(indx) ;
+ w_stage1 = weights(indx)./tau_tilde(indx);
% draw in the new distributions
- wtilde = zeros(1,number_of_particles) ;
+ wtilde = zeros(1, number_of_particles);
for i=1:number_of_particles
- info=1 ;
- while info==1
- 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 = 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
+ [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;
@@ -191,71 +174,75 @@ for t=1:sample_size
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
% Get covariance matrices and structural shocks
- epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations,1) ;
+ epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations, 1);
% compute particles likelihood contribution
- yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
+ 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,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
- StateVectors_(:,i) = tmp_(mf0,:) ;
+ yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state);
+ [tmp, tmp_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2);
+ StateVectors_(:,i) = tmp_(mf0,:);
else
- tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+ tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2);
end
- StateVectors(:,i) = tmp(mf0,:) ;
- PredictionError = bsxfun(@minus,Y(t,:)',tmp(mf1,:));
- wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError),1)));
- end
+ StateVectors(:,i) = tmp(mf0,:);
+ PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
+ wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1)));
+ end
end
end
end
% normalization
weights = wtilde/sum(wtilde);
- if (variance_update==1) && (neff(weights)=0.025 && pass1==1
- lb95_xparam(i,t) = temp(j,1) ;
- pass1 = 2 ;
+ if ~pass1 && cumulated_weights(j)>=0.025
+ lb95_xparam(i,t) = temp(j,1);
+ pass1 = true;
end
- if cumulated_weights(j)>=0.5 && pass2==1
- median_xparam(i,t) = temp(j,1) ;
- pass2 = 2 ;
+ if ~pass2 && cumulated_weights(j)>=0.5
+ median_xparam(i,t) = temp(j,1);
+ pass2 = true;
end
- if cumulated_weights(j)>=0.975 && pass3==1
- ub95_xparam(i,t) = temp(j,1) ;
- pass3 = 2 ;
+ if ~pass3 && cumulated_weights(j)>=0.975
+ ub95_xparam(i,t) = temp(j,1);
+ 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,15 +295,17 @@ 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
drawnow
end
- dyn_saveas(hh,[ Model.fname '_param_traj' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
+ dyn_saveas(hh, [Model.fname '_param_traj' int2str(plt)], DynareOptions.nodisplay, DynareOptions.graph_format);
if TeX
% TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n');
@@ -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,12 +343,12 @@ 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,...
- number_of_particles,optimal_bandwidth,kernel_function);
- plot(density(:,1),density(:,2));
+ 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
- title(name,'interpreter','none')
+ title(name, 'interpreter', 'none')
hold off
axis tight
drawnow
@@ -369,9 +356,9 @@ for plt = 1:nbplt,
dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
if TeX
% TeX eps loader file
- fprintf(fidTeX,'\\begin{figure}[H]\n');
+ fprintf(fidTeX, '\\begin{figure}[H]\n');
for jj = 1:length(x)
- fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
+ fprintf(fidTeX, '\\psfrag{%s}[1][][0.5][0]{%s}\n', deblank(NAMES(jj,:)), deblank(TeXNAMES(jj,:)));
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParametersDensities%s}\n',Model.fname,int2str(plt));
diff --git a/nonlinear-filters/src/solve_model_for_online_filter.m b/nonlinear-filters/src/solve_model_for_online_filter.m
index 74d2fd925..4df219d87 100644
--- a/nonlinear-filters/src/solve_model_for_online_filter.m
+++ b/nonlinear-filters/src/solve_model_for_online_filter.m
@@ -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 .
-% 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(xparam1BayesInfo.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');
+warning('off', 'MATLAB:nearlySingularMatrix')
+[~, ~, ~, info, Model, DynareOptions, DynareResults] = ...
+ dynare_resolve(Model, DynareOptions, DynareResults, 'restrict');
+warning('on', 'MATLAB:nearlySingularMatrix')
-%disp(info)
-
-if info(1) ~= 0
- ReducedForm = 0 ;
- exit_flag = 55;
+if info(1)~=0
+ if nargout==5
+ ReducedForm = 0;
+ 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,37 +144,37 @@ 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);
+ state_variables_idx = restrict_variables_idx(mf0);
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
-ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
-ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
-ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
-if DynareOptions.order>1
- ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
- ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
- ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
- ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
-else
- ReducedForm.ghxx = zeros(size(restrict_variables_idx,1),size(dr.kstate,2));
- ReducedForm.ghuu = zeros(size(restrict_variables_idx,1),size(dr.ghu,2));
- ReducedForm.ghxu = zeros(size(restrict_variables_idx,1),size(dr.ghx,2));
- ReducedForm.constant = ReducedForm.steadystate ;
-end
-ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
-ReducedForm.Q = Q;
-ReducedForm.H = H;
-ReducedForm.mf0 = mf0;
-ReducedForm.mf1 = mf1;
-% Set initial condition for t=1
-if observation_number==1
+% 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));
+ if DynareOptions.order>1
+ ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
+ ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
+ ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
+ ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
+ else
+ ReducedForm.ghxx = zeros(size(restrict_variables_idx,1),size(dr.kstate,2));
+ ReducedForm.ghuu = zeros(size(restrict_variables_idx,1),size(dr.ghu,2));
+ ReducedForm.ghxu = zeros(size(restrict_variables_idx,1),size(dr.ghx,2));
+ ReducedForm.constant = ReducedForm.steadystate ;
+ end
+ ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
+ ReducedForm.Q = Q;
+ ReducedForm.H = H;
+ ReducedForm.mf0 = mf0;
+ ReducedForm.mf1 = mf1;
+end
+
+% 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);
@@ -347,4 +196,4 @@ if observation_number==1
end
ReducedForm.StateVectorMean = StateVectorMean;
ReducedForm.StateVectorVariance = StateVectorVariance;
-end
+end
\ No newline at end of file