sequential_importance_particle_filter

 Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).

0001 function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,DynareOptions)
0002 % Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
0003 
0004 %@info:
0005 %! @deftypefn {Function File} {@var{y}, @var{y_} =} sequential_importance_particle_filter (@var{ReducedForm},@var{Y}, @var{start}, @var{DynareOptions})
0006 %! @anchor{particle/sequential_importance_particle_filter}
0007 %! @sp 1
0008 %! Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
0009 %!
0010 %! @sp 2
0011 %! @strong{Inputs}
0012 %! @sp 1
0013 %! @table @ @var
0014 %! @item ReducedForm
0015 %! Structure describing the state space model (built in @ref{non_linear_dsge_likelihood}).
0016 %! @item Y
0017 %! p*smpl matrix of doubles (p is the number of observed variables), the (detrended) data.
0018 %! @item start
0019 %! Integer scalar, likelihood evaluation starts at observation 'start'.
0020 %! @item DynareOptions
0021 %! Structure specifying Dynare's options.
0022 %! @end table
0023 %! @sp 2
0024 %! @strong{Outputs}
0025 %! @sp 1
0026 %! @table @ @var
0027 %! @item LIK
0028 %! double scalar, value of (minus) the logged likelihood.
0029 %! @item lik
0030 %! smpl*1 vector of doubles, density of the observations at each period.
0031 %! @end table
0032 %! @sp 2
0033 %! @strong{This function is called by:}
0034 %! @ref{non_linear_dsge_likelihood}
0035 %! @sp 2
0036 %! @strong{This function calls:}
0037 %!
0038 %! @end deftypefn
0039 %@eod:
0040 
0041 % Copyright (C) 2011, 2012 Dynare Team
0042 %
0043 % This file is part of Dynare.
0044 %
0045 % Dynare is free software: you can redistribute it and/or modify
0046 % it under the terms of the GNU General Public License as published by
0047 % the Free Software Foundation, either version 3 of the License, or
0048 % (at your option) any later version.
0049 %
0050 % Dynare is distributed in the hope that it will be useful,
0051 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0052 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0053 % GNU General Public License for more details.
0054 %
0055 % You should have received a copy of the GNU General Public License
0056 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0057 
0058 % AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr
0059 %           stephane DOT adjemian AT univ DASH lemans DOT fr
0060 
0061 persistent init_flag
0062 persistent mf0 mf1
0063 persistent number_of_particles
0064 persistent sample_size number_of_observed_variables number_of_structural_innovations
0065 
0066 % Set default value for start
0067 if isempty(start)
0068     start = 1;
0069 end
0070 
0071 % Get steady state and mean.
0072 steadystate = ReducedForm.steadystate;
0073 constant = ReducedForm.constant;
0074 state_variables_steady_state = ReducedForm.state_variables_steady_state;
0075 
0076 % Set persistent variables (if needed).
0077 if isempty(init_flag)
0078     mf0 = ReducedForm.mf0;
0079     mf1 = ReducedForm.mf1;
0080     sample_size = size(Y,2);
0081     number_of_observed_variables = length(mf1);
0082     number_of_structural_innovations = length(ReducedForm.Q);
0083     number_of_particles = DynareOptions.particle.number_of_particles;
0084     init_flag = 1;
0085 end
0086 
0087 % Set local state space model (first order approximation).
0088 ghx  = ReducedForm.ghx;
0089 ghu  = ReducedForm.ghu;
0090 
0091 % Set local state space model (second order approximation).
0092 ghxx = ReducedForm.ghxx;
0093 ghuu = ReducedForm.ghuu;
0094 ghxu = ReducedForm.ghxu;
0095 
0096 % Get covariance matrices.
0097 Q = ReducedForm.Q;
0098 H = ReducedForm.H;
0099 if isempty(H)
0100     H = 0;
0101 end
0102 
0103 % Get initial condition for the state vector.
0104 StateVectorMean = ReducedForm.StateVectorMean;
0105 StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
0106 state_variance_rank = size(StateVectorVarianceSquareRoot,2);
0107 Q_lower_triangular_cholesky = chol(Q)';
0108 
0109 % Set seed for randn().
0110 set_dynare_seed('default');
0111 
0112 % Initialization of the likelihood.
0113 const_lik = log(2*pi)*number_of_observed_variables;
0114 lik  = NaN(sample_size,1);
0115 
0116 % Initialization of the weights across particles.
0117 nb_obs_resamp = 0 ;
0118 weights = ones(1,number_of_particles)/number_of_particles ;
0119 StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
0120 for t=1:sample_size
0121     yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
0122     epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
0123     tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
0124     PredictedObservedMean = mean(tmp(mf1,:),2);
0125     PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
0126     dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
0127     PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles+H;
0128     lnw = -.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1));
0129     dfac = max(lnw);
0130     wtilde = weights.*exp(lnw-dfac);
0131     lik(t) = log(sum(wtilde))+dfac;
0132     weights = wtilde/sum(wtilde);
0133     if strcmp(DynareOptions.particle.resampling.status,'generic')
0134         Neff = 1/(weights*weights');
0135     end
0136     if (strcmp(DynareOptions.particle.resampling.status,'generic') && Neff<DynareOptions.particle.resampling.neff_threshold*sample_size ) || strcmp(DynareOptions.particle.resampling.status,'systematic')
0137         nb_obs_resamp = nb_obs_resamp+1 ;
0138         StateVectors = tmp(mf0,resample(weights,DynareOptions.particle.resampling.method1,DynareOptions.particle.resampling.method2));
0139         weights = ones(1,number_of_particles)/number_of_particles;
0140     elseif strcmp(DynareOptions.particle.resampling.status,'smoothed')
0141         StateVectors = multivariate_smooth_resampling(weights',tmp(mf0,:)',number_of_particles,DynareOptions.particle.resampling.number_of_partitions)';
0142         weights = ones(1,number_of_particles)/number_of_particles;
0143     elseif strcmp(DynareOptions.particle.resampling.status,'none')
0144         StateVectors = tmp(mf0,:);
0145     end
0146 end
0147 
0148 LIK = -sum(lik(start:end));