Put back pruning in particle filter (implemented in sequential_importance_particle_filter).

matlab/particle/sequential_importance_particle_filter.m

@@ -68,6 +68,9 @@ 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;
@@ -103,7 +106,15 @@ end
 % Get initial condition for the state vector.
 StateVectorMean = ReducedForm.StateVectorMean;
 StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
+if pruning
+    StateVectorMean_ = StateVectorMean;
+    StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
+end
+
+% Get the rank of StateVectorVarianceSquareRoot
 state_variance_rank = size(StateVectorVarianceSquareRoot,2);
+
+% Factorize the covariance matrix of the structural innovations
 Q_lower_triangular_cholesky = chol(Q)';
 
 % Set seed for randn().
@@ -114,35 +125,55 @@ const_lik = log(2*pi)*number_of_observed_variables;
 lik  = NaN(sample_size,1);
 
 % Initialization of the weights across particles.
-nb_obs_resamp = 0 ;
 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
+
+% Loop over observations
 for t=1:sample_size
     yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
     epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
+    if pruning
+        yhat_ = bsxfun(@minus,StateVectors_,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);
+    else
         tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+    end
     PredictedObservedMean = tmp(mf1,:)*transpose(weights);
     PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
     dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
-    PredictedObservedVariance = dPredictedObservedMean*diag(weights)*dPredictedObservedMean' + H;
+    PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H;
     lnw = -.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1));
     dfac = max(lnw);
     wtilde = weights.*exp(lnw-dfac);
     lik(t) = log(sum(wtilde))+dfac;
     weights = wtilde/sum(wtilde);
-    if strcmp(DynareOptions.particle.resampling.status,'generic')
-        Neff = 1/(weights*weights');
+    if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || strcmp(DynareOptions.particle.resampling.status,'systematic')
+        idx = resample(weights,DynareOptions.particle.resampling.method1,DynareOptions.particle.resampling.method2);
+        StateVectors = tmp(mf0,idx);
+        if pruning
+            StateVectors_ = tmp_(mf0,idx);
         end
-    if (strcmp(DynareOptions.particle.resampling.status,'generic') && Neff<DynareOptions.particle.resampling.neff_threshold*sample_size ) || strcmp(DynareOptions.particle.resampling.status,'systematic')
-        nb_obs_resamp = nb_obs_resamp+1 ;
-        StateVectors = tmp(mf0,resample(weights,DynareOptions.particle.resampling.method1,DynareOptions.particle.resampling.method2));
         weights = ones(1,number_of_particles)/number_of_particles;
     elseif strcmp(DynareOptions.particle.resampling.status,'smoothed')
         StateVectors = multivariate_smooth_resampling(weights',tmp(mf0,:)',number_of_particles,DynareOptions.particle.resampling.number_of_partitions)';
+        if pruning
+            StateVectors_ = multivariate_smooth_resampling(weights',tmp_(mf0,:)',number_of_particles,DynareOptions.particle.resampling.number_of_partitions)';
+        end
         weights = ones(1,number_of_particles)/number_of_particles;
     elseif strcmp(DynareOptions.particle.resampling.status,'none')
         StateVectors = tmp(mf0,:);
+        if pruning
+            StateVectors_ = tmp_(mf0,:)
+        end
     end
 end
 
 LIK = -sum(lik(start:end));
+
+
+
+function n = neff(w)
+    n = dot(w,w);
\ No newline at end of file