Skip to content
GitLab
Commit 8d401a22 authored by Frédéric Karamé's avatar Frédéric Karamé
Browse files

Fix the calculation of the likelihood on the APF.

parent 5922f881
Hide whitespace changes
Inline Side-by-side
src/auxiliary_particle_filter.m
View file @ 8d401a22
function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions)
% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling.
% Copyright (C) 2011-2014 Dynare Team
% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter
% allowing eventually resampling.
%
% Copyright (C) 2011-2015 Dynare Team
%
% This file is part of Dynare (particles module).
%
......@@ -20,7 +21,7 @@ function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptio
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent init_flag mf0 mf1 number_of_particles
persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
persistent sample_size number_of_observed_variables number_of_structural_innovations
% Set default
if isempty(start)
......@@ -40,7 +41,6 @@ if isempty(init_flag)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
number_of_particles = ParticleOptions.number_of_particles;
......@@ -58,55 +58,44 @@ ghxu = ReducedForm.ghxu;
% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
if isempty(H)
H = 0;
end
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';%reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = chol(Q)';
if pruning
StateVectorMean_ = StateVectorMean;
StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
end
% Set seed for randn().
set_dynare_seed('default');
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables +log(det(H));
const_lik = log(2*pi)*number_of_observed_variables+log(det(H));
lik = NaN(sample_size,1);
LIK = NaN;
% Initialization of the weights across particles.
weights = ones(1,number_of_particles)/number_of_particles ;
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
%StateVectors = bsxfun(@plus,zeros(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
StateVectors_ = StateVectors;
end
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
yhat = zeros(2,number_of_particles) ;
if pruning
yhat_ = zeros(2,number_of_particles) ;
[tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
StateVectors_ = tmp_(mf0,:);
else
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
end
StateVectors = tmp(mf0,:) ;
if ParticleOptions.proposal_approximation.cubature
[nodes,nodes_weights] = spherical_radial_sigma_points(number_of_structural_innovations) ;
nodes_weights = ones(size(nodes,1),1)*nodes_weights ;
elseif ParticleOptions.proposal_approximation.unscented
[nodes,nodes_weights,nodes_weights_c] = unscented_sigma_points(number_of_structural_innovations,ParticleOptions);
else
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
nodes = Q_lower_triangular_cholesky*nodes ;
% Uncomment for building the mean average predictions based on a sparse
% grids of structural shocks. Otherwise, all shocks are set to 0 in the
% prediction.
%if ParticleOptions.proposal_approximation.cubature
% [nodes,nodes_weights] = spherical_radial_sigma_points(number_of_structural_innovations) ;
% nodes_weights = ones(size(nodes,1),1)*nodes_weights ;
%elseif ParticleOptions.proposal_approximation.unscented
% [nodes,nodes_weights,nodes_weights_c] = unscented_sigma_points(number_of_structural_innovations,ParticleOptions);
%else
% error('Estimation: This approximation for the proposal is not implemented or unknown!')
%end
%nodes = Q_lower_triangular_cholesky*nodes ;
nodes = zeros(1,number_of_structural_innovations) ;
nodes_weights = 1 ;
for t=1:sample_size
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
......@@ -125,21 +114,19 @@ for t=1:sample_size
tmp = tmp + nodes_weights(i)*local_state_space_iteration_2(yhat,nodes(i,:)*ones(1,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
end
end
%PredictedObservedMean = weights*(tmp(mf1,:)');
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
%dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean');
%PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' +H;
wtilde = exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) ;
tau_tilde = weights.*wtilde ;
sum_tau_tilde = sum(tau_tilde) ;
lik(t) = log(sum_tau_tilde) ;
tau_tilde = tau_tilde/sum_tau_tilde;
%tau_tilde = weights.*(exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) + 1e-99) ;
% Replace Gaussian density with a Student density with 3 degrees of
% freedom for fat tails.
z = sum(PredictionError.*(H\PredictionError),1) ;
tau_tilde = weights.*(tpdf(z,3*ones(size(z)))+1e-99) ;
tau_tilde = tau_tilde/sum(tau_tilde) ;
indx = resample(0,tau_tilde',ParticleOptions);
if pruning
yhat_ = yhat_(:,indx) ;
end
yhat = yhat(:,indx) ;
factor = weights(indx)./tau_tilde(indx) ;
weights_stage_1 = weights(indx)./tau_tilde(indx) ;
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
if pruning
[tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
......@@ -148,13 +135,21 @@ for t=1:sample_size
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
end
StateVectors = tmp(mf0,:);
%PredictedObservedMean = mean(tmp(mf1,:),2);
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
%dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
%PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles + H;
lnw = exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1)));
wtilde = lnw.*factor ;
weights = wtilde/sum(wtilde);
weights_stage_2 = weights_stage_1.*(exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) + 1e-99) ;
lik(t) = log(mean(weights_stage_2)) ;
weights = weights_stage_2/sum(weights_stage_2);
if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
if pruning
temp = resample([StateVectors' StateVectors_'],weights',ParticleOptions);
StateVectors = temp(:,1:number_of_state_variables)';
StateVectors_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)';
else
StateVectors = resample(StateVectors',weights',ParticleOptions)';
end
weights = ones(1,number_of_particles)/number_of_particles;
end
end
%plot(lik) ;
LIK = -sum(lik(start:end));
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment