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 changesrm-particles^2
parent
22be3797a8
commit
a0fb5c7348
|
@ -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.
|
% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% ReducedForm [structure] Matlab's structure describing the reduced form model.
|
% - xparam1 [double] n×1 vector, Initial condition for the estimated parameters.
|
||||||
% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors.
|
% - DynareDataset [dseries] Sample used for estimation.
|
||||||
% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors.
|
% - dataset_info [struct] Description of the sample.
|
||||||
% ReducedForm.state.dr [structure] output of resol.m.
|
% - DynareOptions [struct] Option values (options_).
|
||||||
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
|
% - Model [struct] Description of the model (M_).
|
||||||
% start [integer] scalar, likelihood evaluation starts at 'start'.
|
% - EstimatedParameters [struct] Description of the estimated parameters (estim_params_).
|
||||||
% mf [integer] pp*1 vector of indices.
|
% - BayesInfo [struct] Prior definition (bayestopt_).
|
||||||
% number_of_particles [integer] scalar.
|
% - DynareResults [struct] Results (oo_).
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% LIK [double] scalar, likelihood
|
% - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters).
|
||||||
% lik [double] vector, density of observations in each period.
|
% - 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).
|
||||||
% REFERENCES
|
% - 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.
|
||||||
% NOTES
|
% - p975 [double] n×1 vector, 97.5 percent of the particles are below p975(i) for i=1,…,n.
|
||||||
% The vector "lik" is used to evaluate the jacobian of the likelihood.
|
% - 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.
|
% 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
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
% 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 seed for randn().
|
||||||
set_dynare_seed('default') ;
|
set_dynare_seed('default');
|
||||||
pruning = DynareOptions.particle.pruning;
|
pruning = DynareOptions.particle.pruning;
|
||||||
second_resample = DynareOptions.particle.resampling.status.systematic ;
|
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
|
% initialization of state particles
|
||||||
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
|
[~, Model, DynareOptions, DynareResults, ReducedForm] = solve_model_for_online_filter(true, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
|
||||||
solve_model_for_online_filter(1,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
|
|
||||||
|
|
||||||
% Set persistent variables.
|
mf0 = ReducedForm.mf0;
|
||||||
if isempty(init_flag)
|
mf1 = ReducedForm.mf1;
|
||||||
mf0 = ReducedForm.mf0;
|
number_of_particles = DynareOptions.particle.number_of_particles;
|
||||||
mf1 = ReducedForm.mf1;
|
number_of_parameters = size(xparam1,1);
|
||||||
number_of_particles = DynareOptions.particle.number_of_particles;
|
Y = DynareDataset.data;
|
||||||
number_of_parameters = size(xparam1,1) ;
|
sample_size = size(Y,1);
|
||||||
Y = DynareDataset.data ;
|
number_of_observed_variables = length(mf1);
|
||||||
sample_size = size(Y,1);
|
number_of_structural_innovations = length(ReducedForm.Q);
|
||||||
number_of_observed_variables = length(mf1);
|
liu_west_delta = DynareOptions.particle.liu_west_delta;
|
||||||
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
|
% Get initial conditions for the state particles
|
||||||
StateVectorMean = ReducedForm.StateVectorMean;
|
StateVectorMean = ReducedForm.StateVectorMean;
|
||||||
|
@ -75,43 +69,34 @@ if pruning
|
||||||
end
|
end
|
||||||
|
|
||||||
% parameters for the Liu & West filter
|
% parameters for the Liu & West filter
|
||||||
small_a = (3*liu_west_delta-1)/(2*liu_west_delta) ;
|
small_a = (3*liu_west_delta-1)/(2*liu_west_delta);
|
||||||
b_square = 1-small_a*small_a ;
|
b_square = 1-small_a*small_a;
|
||||||
|
|
||||||
% Initialization of parameter particles
|
% Initialization of parameter particles
|
||||||
xparam = zeros(number_of_parameters,number_of_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);
|
prior_draw(BayesInfo,DynareOptions.prior_trunc);
|
||||||
for i=1:number_of_particles
|
for i=1:number_of_particles
|
||||||
info = 1;
|
info = 12042009;
|
||||||
while info==1
|
while info
|
||||||
%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)) ;
|
|
||||||
candidate = prior_draw()';
|
candidate = prior_draw()';
|
||||||
if all(candidate(:) >= bounds.lb) && all(candidate(:) <= bounds.ub)
|
[info, Model, DynareOptions, DynareResults] = solve_model_for_online_filter(false, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
|
||||||
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
|
if ~info
|
||||||
solve_model_for_online_filter(1,candidate(:),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
|
xparam(:,i) = candidate(:);
|
||||||
if info==0
|
end
|
||||||
xparam(:,i) = candidate(:) ;
|
|
||||||
end
|
|
||||||
end
|
|
||||||
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.
|
% 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.
|
% Initialization of the likelihood.
|
||||||
const_lik = log(2*pi)*number_of_observed_variables;
|
const_lik = log(2*pi)*number_of_observed_variables;
|
||||||
mean_xparam = zeros(number_of_parameters,sample_size) ;
|
mean_xparam = zeros(number_of_parameters,sample_size);
|
||||||
median_xparam = zeros(number_of_parameters,sample_size) ;
|
mode_xparam = zeros(number_of_parameters,sample_size);
|
||||||
std_xparam = zeros(number_of_parameters,sample_size) ;
|
median_xparam = zeros(number_of_parameters,sample_size);
|
||||||
lb95_xparam = zeros(number_of_parameters,sample_size) ;
|
std_xparam = zeros(number_of_parameters,sample_size);
|
||||||
ub95_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
|
%% The Online filter
|
||||||
for t=1:sample_size
|
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))
|
fprintf('\nSubsample with only the first observation.\n\n', int2str(t))
|
||||||
end
|
end
|
||||||
% Moments of parameters particles distribution
|
% Moments of parameters particles distribution
|
||||||
m_bar = xparam*(weights') ;
|
m_bar = xparam*(weights');
|
||||||
temp = bsxfun(@minus,xparam,m_bar) ;
|
temp = bsxfun(@minus,xparam,m_bar);
|
||||||
sigma_bar = (bsxfun(@times,weights,temp))*(temp') ;
|
sigma_bar = (bsxfun(@times,weights,temp))*(temp');
|
||||||
if variance_update==1
|
if variance_update
|
||||||
chol_sigma_bar = chol(b_square*sigma_bar)' ;
|
chol_sigma_bar = chol(b_square*sigma_bar)';
|
||||||
end
|
end
|
||||||
% Prediction (without shocks)
|
% Prediction (without shocks)
|
||||||
fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ;
|
fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam);
|
||||||
tau_tilde = zeros(1,number_of_particles) ;
|
tau_tilde = zeros(1,number_of_particles);
|
||||||
for i=1:number_of_particles
|
for i=1:number_of_particles
|
||||||
% model resolution
|
% model resolution
|
||||||
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
|
[info, Model, DynareOptions, DynareResults, ReducedForm] = ...
|
||||||
solve_model_for_online_filter(t,fore_xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
|
solve_model_for_online_filter(false, fore_xparam(:,i), DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
|
||||||
if info==0
|
if ~info
|
||||||
steadystate = ReducedForm.steadystate;
|
steadystate = ReducedForm.steadystate;
|
||||||
state_variables_steady_state = ReducedForm.state_variables_steady_state;
|
state_variables_steady_state = ReducedForm.state_variables_steady_state;
|
||||||
% Set local state space model (second-order approximation).
|
% 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);
|
yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
|
||||||
if pruning
|
if pruning
|
||||||
yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
|
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
|
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
|
end
|
||||||
PredictionError = bsxfun(@minus,Y(t,:)',tmp(mf1,:));
|
PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
|
||||||
% Replace Gaussian density with a Student density with 3 degrees of
|
% Replace Gaussian density with a Student density with 3 degrees of freedom for fat tails.
|
||||||
% freedom for fat tails.
|
z = sum(PredictionError.*(ReducedForm.H\PredictionError), 1) ;
|
||||||
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).*(tpdf(z,3*ones(size(z)))+1e-99) ;
|
end
|
||||||
%tau_tilde(i) = weights(i).*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError),1))) ;
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
% particles selection
|
% particles selection
|
||||||
tau_tilde = tau_tilde/sum(tau_tilde) ;
|
tau_tilde = tau_tilde/sum(tau_tilde);
|
||||||
indx = resample(0,tau_tilde',DynareOptions.particle);
|
indx = resample(0, tau_tilde', DynareOptions.particle);
|
||||||
StateVectors = StateVectors(:,indx) ;
|
StateVectors = StateVectors(:,indx);
|
||||||
xparam = fore_xparam(:,indx);
|
xparam = fore_xparam(:,indx);
|
||||||
if pruning
|
if pruning
|
||||||
StateVectors_ = StateVectors_(:,indx) ;
|
StateVectors_ = StateVectors_(:,indx);
|
||||||
end
|
end
|
||||||
w_stage1 = weights(indx)./tau_tilde(indx) ;
|
w_stage1 = weights(indx)./tau_tilde(indx);
|
||||||
% draw in the new distributions
|
% draw in the new distributions
|
||||||
wtilde = zeros(1,number_of_particles) ;
|
wtilde = zeros(1, number_of_particles);
|
||||||
for i=1:number_of_particles
|
for i=1:number_of_particles
|
||||||
info=1 ;
|
info = 12042009;
|
||||||
while info==1
|
while info
|
||||||
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters,1) ;
|
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
|
||||||
if all(candidate >= bounds.lb) && all(candidate <= bounds.ub)
|
if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
|
||||||
% model resolution for new parameters particles
|
% model resolution for new parameters particles
|
||||||
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
|
[info, Model, DynareOptions, DynareResults, ReducedForm] = ...
|
||||||
solve_model_for_online_filter(t,candidate,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
|
solve_model_for_online_filter(false, candidate, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults) ;
|
||||||
if info==0
|
if ~info
|
||||||
xparam(:,i) = candidate ;
|
xparam(:,i) = candidate ;
|
||||||
steadystate = ReducedForm.steadystate;
|
steadystate = ReducedForm.steadystate;
|
||||||
state_variables_steady_state = ReducedForm.state_variables_steady_state;
|
state_variables_steady_state = ReducedForm.state_variables_steady_state;
|
||||||
|
@ -191,71 +174,75 @@ for t=1:sample_size
|
||||||
ghuu = ReducedForm.ghuu;
|
ghuu = ReducedForm.ghuu;
|
||||||
ghxu = ReducedForm.ghxu;
|
ghxu = ReducedForm.ghxu;
|
||||||
% Get covariance matrices and structural shocks
|
% 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
|
% compute particles likelihood contribution
|
||||||
yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
|
yhat = bsxfun(@minus,StateVectors(:,i), state_variables_steady_state);
|
||||||
if pruning
|
if pruning
|
||||||
yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
|
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);
|
[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,:) ;
|
StateVectors_(:,i) = tmp_(mf0,:);
|
||||||
else
|
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
|
end
|
||||||
StateVectors(:,i) = tmp(mf0,:) ;
|
StateVectors(:,i) = tmp(mf0,:);
|
||||||
PredictionError = bsxfun(@minus,Y(t,:)',tmp(mf1,:));
|
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)));
|
wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1)));
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
% normalization
|
% normalization
|
||||||
weights = wtilde/sum(wtilde);
|
weights = wtilde/sum(wtilde);
|
||||||
if (variance_update==1) && (neff(weights)<DynareOptions.particle.resampling.threshold*sample_size)
|
if variance_update && (neff(weights)<DynareOptions.particle.resampling.threshold*sample_size)
|
||||||
variance_update = 0 ;
|
variance_update = false;
|
||||||
end
|
end
|
||||||
% final resampling (not advised)
|
% final resampling (not advised)
|
||||||
if second_resample==1
|
if second_resample
|
||||||
indx = resample(0,weights,DynareOptions.particle);
|
[~, idmode] = max(weights);
|
||||||
|
mode_xparam(:,t) = xparam(:,idmode);
|
||||||
|
indx = resample(0, weights,DynareOptions.particle);
|
||||||
StateVectors = StateVectors(:,indx) ;
|
StateVectors = StateVectors(:,indx) ;
|
||||||
if pruning
|
if pruning
|
||||||
StateVectors_ = StateVectors_(:,indx) ;
|
StateVectors_ = StateVectors_(:,indx);
|
||||||
end
|
end
|
||||||
xparam = xparam(:,indx) ;
|
xparam = xparam(:,indx);
|
||||||
weights = ones(1,number_of_particles)/number_of_particles ;
|
weights = ones(1, number_of_particles)/number_of_particles;
|
||||||
mean_xparam(:,t) = mean(xparam,2);
|
mean_xparam(:,t) = mean(xparam, 2);
|
||||||
mat_var_cov = bsxfun(@minus,xparam,mean_xparam(:,t)) ;
|
mat_var_cov = bsxfun(@minus, xparam, mean_xparam(:,t));
|
||||||
mat_var_cov = (mat_var_cov*mat_var_cov')/(number_of_particles-1) ;
|
mat_var_cov = (mat_var_cov*mat_var_cov')/(number_of_particles-1);
|
||||||
std_xparam(:,t) = sqrt(diag(mat_var_cov)) ;
|
std_xparam(:,t) = sqrt(diag(mat_var_cov));
|
||||||
for i=1:number_of_parameters
|
for i=1:number_of_parameters
|
||||||
temp = sortrows(xparam(i,:)') ;
|
temp = sortrows(xparam(i,:)');
|
||||||
lb95_xparam(i,t) = temp(0.025*number_of_particles) ;
|
lb95_xparam(i,t) = temp(0.025*number_of_particles);
|
||||||
median_xparam(i,t) = temp(0.5*number_of_particles) ;
|
median_xparam(i,t) = temp(0.5*number_of_particles);
|
||||||
ub95_xparam(i,t) = temp(0.975*number_of_particles) ;
|
ub95_xparam(i,t) = temp(0.975*number_of_particles);
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
if second_resample==0
|
if second_resample
|
||||||
mean_xparam(:,t) = xparam*(weights') ;
|
[~, idmode] = max(weights);
|
||||||
mat_var_cov = bsxfun(@minus,xparam,mean_xparam(:,t)) ;
|
mode_xparam(:,t) = xparam(:,idmode);
|
||||||
mat_var_cov = mat_var_cov*(bsxfun(@times,mat_var_cov,weights)') ;
|
mean_xparam(:,t) = xparam*(weights');
|
||||||
std_xparam(:,t) = sqrt(diag(mat_var_cov)) ;
|
mat_var_cov = bsxfun(@minus, xparam,mean_xparam(:,t));
|
||||||
|
mat_var_cov = mat_var_cov*(bsxfun(@times, mat_var_cov, weights)');
|
||||||
|
std_xparam(:,t) = sqrt(diag(mat_var_cov));
|
||||||
for i=1:number_of_parameters
|
for i=1:number_of_parameters
|
||||||
temp = sortrows([xparam(i,:)' weights'],1) ;
|
temp = sortrows([xparam(i,:)' weights'], 1);
|
||||||
cumulated_weights = cumsum(temp(:,2)) ;
|
cumulated_weights = cumsum(temp(:,2));
|
||||||
pass1=1 ;
|
pass1 = false;
|
||||||
pass2=1 ;
|
pass2 = false;
|
||||||
pass3=1 ;
|
pass3 = false;
|
||||||
for j=1:number_of_particles
|
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) ;
|
lb95_xparam(i,t) = temp(j,1);
|
||||||
pass1 = 2 ;
|
pass1 = true;
|
||||||
end
|
end
|
||||||
if cumulated_weights(j)>=0.5 && pass2==1
|
if ~pass2 && cumulated_weights(j)>=0.5
|
||||||
median_xparam(i,t) = temp(j,1) ;
|
median_xparam(i,t) = temp(j,1);
|
||||||
pass2 = 2 ;
|
pass2 = true;
|
||||||
end
|
end
|
||||||
if cumulated_weights(j)>=0.975 && pass3==1
|
if ~pass3 && cumulated_weights(j)>=0.975
|
||||||
ub95_xparam(i,t) = temp(j,1) ;
|
ub95_xparam(i,t) = temp(j,1);
|
||||||
pass3 = 2 ;
|
pass3 = true;
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
@ -267,22 +254,22 @@ for t=1:sample_size
|
||||||
disp([str])
|
disp([str])
|
||||||
disp('')
|
disp('')
|
||||||
end
|
end
|
||||||
distrib_param = xparam ;
|
|
||||||
xparam = mean_xparam(:,sample_size) ;
|
pmean = xparam(:,sample_size);
|
||||||
std_param = std_xparam(:,sample_size) ;
|
pmode = mode_xparam(:,sample_size);
|
||||||
lb_95 = lb95_xparam(:,sample_size) ;
|
pstdev = std_xparam(:,sample_size) ;
|
||||||
ub_95 = ub95_xparam(:,sample_size) ;
|
p025 = lb95_xparam(:,sample_size) ;
|
||||||
median_param = median_xparam(:,sample_size) ;
|
p975 = ub95_xparam(:,sample_size) ;
|
||||||
|
pmedian = median_xparam(:,sample_size) ;
|
||||||
|
covariance = mat_var_cov;
|
||||||
|
|
||||||
%% Plot parameters trajectory
|
%% Plot parameters trajectory
|
||||||
TeX = DynareOptions.TeX;
|
TeX = DynareOptions.TeX;
|
||||||
|
|
||||||
[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters);
|
|
||||||
nr = ceil(sqrt(number_of_parameters)) ;
|
nr = ceil(sqrt(number_of_parameters)) ;
|
||||||
nc = floor(sqrt(number_of_parameters));
|
nc = floor(sqrt(number_of_parameters));
|
||||||
nbplt = 1 ;
|
nbplt = 1 ;
|
||||||
|
|
||||||
|
|
||||||
if TeX
|
if TeX
|
||||||
fidTeX = fopen([Model.fname '_param_traj.tex'],'w');
|
fidTeX = fopen([Model.fname '_param_traj.tex'],'w');
|
||||||
fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
|
fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
|
||||||
|
@ -290,15 +277,13 @@ if TeX
|
||||||
fprintf(fidTeX,' \n');
|
fprintf(fidTeX,' \n');
|
||||||
end
|
end
|
||||||
|
|
||||||
z = 1:1:sample_size ;
|
for plt = 1:nbplt
|
||||||
|
|
||||||
for plt = 1:nbplt,
|
|
||||||
if TeX
|
if TeX
|
||||||
NAMES = [];
|
NAMES = [];
|
||||||
TeXNAMES = [];
|
TeXNAMES = [];
|
||||||
end
|
end
|
||||||
hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories');
|
hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories');
|
||||||
for k=1:length(xparam)
|
for k=1:length(pmean)
|
||||||
subplot(nr,nc,k)
|
subplot(nr,nc,k)
|
||||||
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
|
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
|
||||||
if TeX
|
if TeX
|
||||||
|
@ -310,15 +295,17 @@ for plt = 1:nbplt,
|
||||||
TeXNAMES = char(TeXNAMES,texname);
|
TeXNAMES = char(TeXNAMES,texname);
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
y = [mean_xparam(k,:)' median_xparam(k,:)' lb95_xparam(k,:)' ub95_xparam(k,:)' xparam(k)*ones(sample_size,1)] ;
|
% Draw the surface for an interval containing 95% of the particles.
|
||||||
plot(z,y);
|
shade(1:sample_size, ub95_xparam(k,:)', 1:sample_size, lb95_xparam(k,:)', 'FillType',[1 2], 'LineStyle', 'none', 'Marker', 'none')
|
||||||
hold on
|
hold on
|
||||||
|
% Draw the mean of particles.
|
||||||
|
plot(1:sample_size, mean_xparam(k,:), '-k', 'linewidth', 2)
|
||||||
title(name,'interpreter','none')
|
title(name,'interpreter','none')
|
||||||
hold off
|
hold off
|
||||||
axis tight
|
axis tight
|
||||||
drawnow
|
drawnow
|
||||||
end
|
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
|
if TeX
|
||||||
% TeX eps loader file
|
% TeX eps loader file
|
||||||
fprintf(fidTeX,'\\begin{figure}[H]\n');
|
fprintf(fidTeX,'\\begin{figure}[H]\n');
|
||||||
|
@ -334,17 +321,17 @@ for plt = 1:nbplt,
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
%% Plot Parameter Densities
|
% Plot Parameter Densities
|
||||||
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
|
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
|
||||||
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
|
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
|
||||||
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation.
|
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation.
|
||||||
for plt = 1:nbplt,
|
for plt = 1:nbplt
|
||||||
if TeX
|
if TeX
|
||||||
NAMES = [];
|
NAMES = [];
|
||||||
TeXNAMES = [];
|
TeXNAMES = [];
|
||||||
end
|
end
|
||||||
hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities');
|
hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities');
|
||||||
for k=1:length(xparam)
|
for k=1:length(pmean)
|
||||||
subplot(nr,nc,k)
|
subplot(nr,nc,k)
|
||||||
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
|
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
|
||||||
if TeX
|
if TeX
|
||||||
|
@ -356,12 +343,12 @@ for plt = 1:nbplt,
|
||||||
TeXNAMES = char(TeXNAMES,texname);
|
TeXNAMES = char(TeXNAMES,texname);
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',number_of_particles,bandwidth,kernel_function);
|
optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function);
|
||||||
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
|
[density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ...
|
||||||
number_of_particles,optimal_bandwidth,kernel_function);
|
number_of_particles, optimal_bandwidth, kernel_function);
|
||||||
plot(density(:,1),density(:,2));
|
plot(density(:,1), density(:,2));
|
||||||
hold on
|
hold on
|
||||||
title(name,'interpreter','none')
|
title(name, 'interpreter', 'none')
|
||||||
hold off
|
hold off
|
||||||
axis tight
|
axis tight
|
||||||
drawnow
|
drawnow
|
||||||
|
@ -369,9 +356,9 @@ for plt = 1:nbplt,
|
||||||
dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
|
dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
|
||||||
if TeX
|
if TeX
|
||||||
% TeX eps loader file
|
% TeX eps loader file
|
||||||
fprintf(fidTeX,'\\begin{figure}[H]\n');
|
fprintf(fidTeX, '\\begin{figure}[H]\n');
|
||||||
for jj = 1:length(x)
|
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
|
end
|
||||||
fprintf(fidTeX,'\\centering \n');
|
fprintf(fidTeX,'\\centering \n');
|
||||||
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParametersDensities%s}\n',Model.fname,int2str(plt));
|
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParametersDensities%s}\n',Model.fname,int2str(plt));
|
||||||
|
|
|
@ -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)
|
function [info, Model, DynareOptions, DynareResults, ReducedForm] = ...
|
||||||
% solve the dsge model for an particular parameters set.
|
solve_model_for_online_filter(setinitialcondition, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults)
|
||||||
|
|
||||||
%@info:
|
% Solves the dsge model for an particular parameters set.
|
||||||
%! @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}
|
% INPUTS
|
||||||
%! @sp 1
|
% - setinitialcondition [logical] return initial condition if true.
|
||||||
%! Evaluates the posterior kernel of a dsge model using a non linear filter.
|
% - xparam1 [double] n×1 vector, parameter values.
|
||||||
%! @sp 2
|
% - DynareDataset [struct] Dataset for estimation (dataset_).
|
||||||
%! @strong{Inputs}
|
% - DynareOptions [struct] Dynare options (options_).
|
||||||
%! @sp 1
|
% - Model [struct] Model description (M_).
|
||||||
%! @table @ @var
|
% - EstimatedParameters [struct] Estimated parameters (estim_params_).
|
||||||
%! @item xparam1
|
% - BayesInfo [struct] Prior definition (bayestopt_).
|
||||||
%! Vector of doubles, current values for the estimated parameters.
|
% - DynareResults [struct] Dynare results (oo_).
|
||||||
%! @item DynareDataset
|
%
|
||||||
%! Matlab's structure describing the dataset (initialized by dynare, see @ref{dataset_}).
|
% OUTPUTS
|
||||||
%! @item DynareOptions
|
% - info [integer] scalar, nonzero if any problem occur when computing the reduced form.
|
||||||
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
|
% - Model [struct] Model description (M_).
|
||||||
%! @item Model
|
% - DynareOptions [struct] Dynare options (options_).
|
||||||
%! Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
|
% - DynareResults [struct] Dynare results (oo_).
|
||||||
%! @item EstimatedParamemeters
|
% - ReducedForm [struct] Reduced form model.
|
||||||
%! 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:
|
|
||||||
|
|
||||||
% Copyright (C) 2013-2017 Dynare Team
|
% Copyright (C) 2013-2019 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% 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
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
|
persistent init_flag restrict_variables_idx state_variables_idx mf0 mf1 number_of_state_variables
|
||||||
% frederic DOT karame AT univ DASH lemans DOT fr
|
|
||||||
|
|
||||||
%global objective_function_penalty_base
|
info = 0;
|
||||||
% 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
|
|
||||||
|
|
||||||
% 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
|
% 1. Get the structural parameters & define penalties
|
||||||
%------------------------------------------------------------------------------
|
%----------------------------------------------------
|
||||||
|
|
||||||
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
|
% Test if some parameters are smaller than the lower bound of the prior domain.
|
||||||
%if (DynareOptions.mode_compute~=1) && any(xparam1<BayesInfo.lb)
|
if any(xparam1<bounds.lb)
|
||||||
% k = find(xparam1(:) < BayesInfo.lb);
|
info = 41;
|
||||||
% fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
|
return
|
||||||
% exit_flag = 0;
|
end
|
||||||
% info = 41;
|
|
||||||
% return
|
|
||||||
%end
|
|
||||||
|
|
||||||
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
|
% Test if some parameters are greater than the upper bound of the prior domain.
|
||||||
%if (DynareOptions.mode_compute~=1) && any(xparam1>BayesInfo.ub)
|
if any(xparam1>bounds.ub)
|
||||||
% k = find(xparam1(:)>BayesInfo.ub);
|
info = 42;
|
||||||
% fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
|
return
|
||||||
% exit_flag = 0;
|
end
|
||||||
% info = 42;
|
|
||||||
% return
|
|
||||||
%end
|
|
||||||
|
|
||||||
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
|
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
|
||||||
Q = Model.Sigma_e;
|
Q = Model.Sigma_e;
|
||||||
|
@ -173,7 +71,7 @@ if EstimatedParameters.nvn
|
||||||
end
|
end
|
||||||
offset = offset+EstimatedParameters.nvn;
|
offset = offset+EstimatedParameters.nvn;
|
||||||
else
|
else
|
||||||
H = zeros(nvobs);
|
H = zeros(size(DynareDataset.data, 1));
|
||||||
end
|
end
|
||||||
|
|
||||||
% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
|
% 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);
|
Q(k2,k1) = Q(k1,k2);
|
||||||
end
|
end
|
||||||
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
|
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
|
||||||
% [CholQ,testQ] = chol(Q);
|
[~, testQ] = chol(Q);
|
||||||
% if testQ
|
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.
|
% The variance-covariance matrix of the structural innovations is not definite positive.
|
||||||
% a = diag(eig(Q));
|
info = 43;
|
||||||
% k = find(a < 0);
|
return
|
||||||
% if k > 0
|
end
|
||||||
% fval = objective_function_penalty_base+sum(-a(k));
|
|
||||||
% exit_flag = 0;
|
|
||||||
% info = 43;
|
|
||||||
% return
|
|
||||||
% end
|
|
||||||
% end
|
|
||||||
offset = offset+EstimatedParameters.ncx;
|
offset = offset+EstimatedParameters.ncx;
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -210,18 +102,12 @@ if EstimatedParameters.ncn
|
||||||
H(k2,k1) = H(k1,k2);
|
H(k2,k1) = H(k1,k2);
|
||||||
end
|
end
|
||||||
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
|
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
|
||||||
% [CholH,testH] = chol(H);
|
[~, testH] = chol(H);
|
||||||
% if testH
|
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.
|
% The variance-covariance matrix of the measurement errors is not definite positive.
|
||||||
% a = diag(eig(H));
|
info = 44;
|
||||||
% k = find(a < 0);
|
return
|
||||||
% if k > 0
|
end
|
||||||
% fval = objective_function_penalty_base+sum(-a(k));
|
|
||||||
% exit_flag = 0;
|
|
||||||
% info = 44;
|
|
||||||
% return
|
|
||||||
% end
|
|
||||||
% end
|
|
||||||
offset = offset+EstimatedParameters.ncn;
|
offset = offset+EstimatedParameters.ncn;
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -238,55 +124,18 @@ Model.H = H;
|
||||||
% 2. call model setup & reduction program
|
% 2. call model setup & reduction program
|
||||||
%------------------------------------------------------------------------------
|
%------------------------------------------------------------------------------
|
||||||
|
|
||||||
% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
|
warning('off', 'MATLAB:nearlySingularMatrix')
|
||||||
[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
|
[~, ~, ~, info, Model, DynareOptions, DynareResults] = ...
|
||||||
|
dynare_resolve(Model, DynareOptions, DynareResults, 'restrict');
|
||||||
|
warning('on', 'MATLAB:nearlySingularMatrix')
|
||||||
|
|
||||||
%disp(info)
|
if info(1)~=0
|
||||||
|
if nargout==5
|
||||||
if info(1) ~= 0
|
ReducedForm = 0;
|
||||||
ReducedForm = 0 ;
|
end
|
||||||
exit_flag = 55;
|
|
||||||
return
|
return
|
||||||
end
|
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.
|
% Get decision rules and transition equations.
|
||||||
dr = DynareResults.dr;
|
dr = DynareResults.dr;
|
||||||
|
|
||||||
|
@ -295,37 +144,37 @@ if isempty(init_flag)
|
||||||
mf0 = BayesInfo.mf0;
|
mf0 = BayesInfo.mf0;
|
||||||
mf1 = BayesInfo.mf1;
|
mf1 = BayesInfo.mf1;
|
||||||
restrict_variables_idx = dr.restrict_var_list;
|
restrict_variables_idx = dr.restrict_var_list;
|
||||||
observed_variables_idx = restrict_variables_idx(mf1);
|
state_variables_idx = restrict_variables_idx(mf0);
|
||||||
state_variables_idx = restrict_variables_idx(mf0);
|
|
||||||
sample_size = size(Y,2);
|
|
||||||
number_of_state_variables = length(mf0);
|
number_of_state_variables = length(mf0);
|
||||||
number_of_observed_variables = length(mf1);
|
init_flag = true;
|
||||||
number_of_structural_innovations = length(Q);
|
|
||||||
init_flag = 1;
|
|
||||||
end
|
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
|
% Return reduced form model.
|
||||||
if observation_number==1
|
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
|
switch DynareOptions.particle.initialization
|
||||||
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
|
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
|
||||||
StateVectorMean = ReducedForm.constant(mf0);
|
StateVectorMean = ReducedForm.constant(mf0);
|
||||||
|
@ -347,4 +196,4 @@ if observation_number==1
|
||||||
end
|
end
|
||||||
ReducedForm.StateVectorMean = StateVectorMean;
|
ReducedForm.StateVectorMean = StateVectorMean;
|
||||||
ReducedForm.StateVectorVariance = StateVectorVariance;
|
ReducedForm.StateVectorVariance = StateVectorVariance;
|
||||||
end
|
end
|
Loading…
Reference in New Issue