Commit 155f6f3b authored by Frédéric Karamé's avatar Frédéric Karamé
Browse files

modification in the call of resampling: only one call now

parent 9df38579
function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,DynareOptions)
% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling.
%
% 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.
%
% 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.
% Copyright (C) 2011, 2012 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent init_flag mf0 mf1 number_of_particles
persistent number_of_observed_variables number_of_structural_innovations
% Set default
if isempty(start)
start = 1;
end
% Set flag for prunning
%pruning = DynareOptions.particle.pruning;
% Get steady state and mean.
%steadystate = ReducedForm.steadystate;
constant = ReducedForm.constant;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set persistent variables.
if isempty(init_flag)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
number_of_particles = DynareOptions.particle.number_of_particles;
init_flag = 1;
end
% Set local state space model (first order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
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 = reduced_rank_cholesky(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;
% 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);
%if pruning
% StateVectors_ = StateVectors;
%end
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
%if pruning
% yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state);
% [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),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,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
%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+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
tau_tilde = weights.*wtilde ;
tau_tilde = tau_tilde/sum(tau_tilde);
initial_distribution = resample(StateVectors',tau_tilde',DynareOptions)' ;
\ No newline at end of file
...@@ -37,7 +37,7 @@ function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,DynareOptions ...@@ -37,7 +37,7 @@ function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,DynareOptions
% 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 init_flag mf0 mf1 number_of_particles persistent init_flag mf0 mf1 number_of_particles
persistent sample_size number_of_observed_variables number_of_structural_innovations persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
% Set default % Set default
if isempty(start) if isempty(start)
...@@ -57,6 +57,7 @@ if isempty(init_flag) ...@@ -57,6 +57,7 @@ if isempty(init_flag)
mf0 = ReducedForm.mf0; mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1; mf1 = ReducedForm.mf1;
sample_size = size(Y,2); sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1); number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q); number_of_structural_innovations = length(ReducedForm.Q);
number_of_particles = DynareOptions.particle.number_of_particles; number_of_particles = DynareOptions.particle.number_of_particles;
...@@ -121,12 +122,25 @@ for t=1:sample_size ...@@ -121,12 +122,25 @@ for t=1:sample_size
%var_wtilde = var_wtilde'*var_wtilde/(number_of_particles-1) ; %var_wtilde = var_wtilde'*var_wtilde/(number_of_particles-1) ;
lik(t) = log(sum_tau_tilde) ; %+ .5*var_wtilde/(number_of_particles*(sum_tau_tilde*sum_tau_tilde)) ; lik(t) = log(sum_tau_tilde) ; %+ .5*var_wtilde/(number_of_particles*(sum_tau_tilde*sum_tau_tilde)) ;
tau_tilde = tau_tilde/sum_tau_tilde; tau_tilde = tau_tilde/sum_tau_tilde;
indx_resmpl = resample(tau_tilde,DynareOptions.particle.resampling.method1,DynareOptions.particle.resampling.method2); if pruning
yhat = yhat(:,indx_resmpl); temp = resample([yhat' yhat_'],tau_tilde',DynareOptions);
wtilde = wtilde(indx_resmpl); yhat = temp(:,1:number_of_state_variables)' ;
yhat_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)' ;
else
yhat = resample(yhat',tau_tilde',DynareOptions)' ;
end
if pruning
[tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),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,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
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+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles); epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
if pruning if pruning
yhat_ = yhat_(:,indx_resmpl);
[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_ = tmp_(mf0,:); StateVectors_ = tmp_(mf0,:);
else else
......
...@@ -149,12 +149,7 @@ for t=1:sample_size ...@@ -149,12 +149,7 @@ for t=1:sample_size
if (Neff<.5*sample_size && strcmpi(DynareOptions.particle.resampling.status,'generic')) || ... if (Neff<.5*sample_size && strcmpi(DynareOptions.particle.resampling.status,'generic')) || ...
strcmpi(DynareOptions.particle.resampling.status,'systematic') strcmpi(DynareOptions.particle.resampling.status,'systematic')
ks = ks + 1 ; ks = ks + 1 ;
StateParticles = StateParticles(:,resample(SampleWeights',DynareOptions.particle.resampling.method1,DynareOptions.particle.resampling.method2)) ; StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ;
StateVectorMean = mean(StateParticles,2) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )';
SampleWeights = 1/number_of_particles ;
elseif strcmp(DynareOptions.particle.resampling.status,'smoothed')
StateParticles = multivariate_smooth_resampling(SampleWeights,StateParticles',number_of_particles,DynareOptions.particle.resampling.number_of_partitions)';
StateVectorMean = mean(StateParticles,2) ; StateVectorMean = mean(StateParticles,2) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )'; StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )';
SampleWeights = 1/number_of_particles ; SampleWeights = 1/number_of_particles ;
......
...@@ -60,7 +60,7 @@ function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,D ...@@ -60,7 +60,7 @@ function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,D
persistent init_flag persistent init_flag
persistent mf0 mf1 persistent mf0 mf1
persistent number_of_particles persistent number_of_particles number_of_state_variables
persistent sample_size number_of_observed_variables number_of_structural_innovations persistent sample_size number_of_observed_variables number_of_structural_innovations
% Set default value for start % Set default value for start
...@@ -81,6 +81,7 @@ if isempty(init_flag) ...@@ -81,6 +81,7 @@ if isempty(init_flag)
mf0 = ReducedForm.mf0; mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1; mf1 = ReducedForm.mf1;
sample_size = size(Y,2); sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1); number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q); number_of_structural_innovations = length(ReducedForm.Q);
number_of_particles = DynareOptions.particle.number_of_particles; number_of_particles = DynareOptions.particle.number_of_particles;
...@@ -151,17 +152,14 @@ for t=1:sample_size ...@@ -151,17 +152,14 @@ for t=1:sample_size
wtilde = weights.*exp(lnw-dfac); wtilde = weights.*exp(lnw-dfac);
lik(t) = log(sum(wtilde))+dfac; lik(t) = log(sum(wtilde))+dfac;
weights = wtilde/sum(wtilde); weights = wtilde/sum(wtilde);
if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || strcmp(DynareOptions.particle.resampling.status,'systematic') if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || ...
idx = resample(weights,DynareOptions.particle.resampling.method1,DynareOptions.particle.resampling.method2); strcmp(DynareOptions.particle.resampling.status,'systematic')
StateVectors = tmp(mf0,idx);
if pruning if pruning
StateVectors_ = tmp_(mf0,idx); temp = resample([tmp(mf0,:)' tmp_(mf0,:)'],weights,DynareOptions);
end StateVectors = temp(:,1:number_of_state_variables)' ;
weights = ones(1,number_of_particles)/number_of_particles; StateVectors_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)';
elseif strcmp(DynareOptions.particle.resampling.status,'smoothed') else
StateVectors = multivariate_smooth_resampling(weights',tmp(mf0,:)',number_of_particles,DynareOptions.particle.resampling.number_of_partitions)'; StateVectors = resample(tmp(mf0,:)',weights,DynareOptions)';
if pruning
StateVectors_ = multivariate_smooth_resampling(weights',tmp_(mf0,:)',number_of_particles,DynareOptions.particle.resampling.number_of_partitions)';
end end
weights = ones(1,number_of_particles)/number_of_particles; weights = ones(1,number_of_particles)/number_of_particles;
elseif strcmp(DynareOptions.particle.resampling.status,'none') elseif strcmp(DynareOptions.particle.resampling.status,'none')
......
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