Correction on a bug and extension of the original paper for the calculation of incremental weights.

src/conditional_filter_proposal.m

@@ -95,6 +95,7 @@ if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_a
     PredictedObservedMean = sum(PredictedObservedMean,2) ;
     dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean)'.*sqrt(weights) ;
     dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean)'.*sqrt(weights);
+    PredictedStateVariance = dState*dState';
     big_mat = [dObserved  dState; [H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)] ];
     [mat1,mat] = qr2(big_mat,0);
     mat = mat';
@@ -112,12 +113,13 @@ else
     KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance ;
     StateVectorMean = PredictedStateMean + KalmanFilterGain*(obs - PredictedObservedMean);
     StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
-    StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
-    StateVectorVarianceSquareRoot = chol(StateVectorVariance)';%reduced_rank_cholesky(StateVectorVariance)';
+    %StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
+    StateVectorVarianceSquareRoot = chol(StateVectorVariance + 1e-6)' ;
 end
 
+PredictedStateVarianceSquareRoot = chol(PredictedStateVariance + 1e-6)'  ;
 ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ;
-ypred = measurement_equations(ProposalStateVector,ReducedForm,ThreadsOptions) ;
-foo = H_lower_triangular_cholesky \ (obs - ypred) ;
-likelihood = exp(-0.5*(foo)'*foo)/normconst2 + 1e-99 ;
-Weights = SampleWeights.*likelihood;
+Prior = probability2(PredictedStateMean,PredictedStateVarianceSquareRoot,ProposalStateVector) ; 
+Posterior = probability2(StateVectorMean,StateVectorVarianceSquareRoot,ProposalStateVector) ; 
+Likelihood = probability2(obs,H_lower_triangular_cholesky,measurement_equations(ProposalStateVector,ReducedForm,ThreadsOptions)) ; 
+Weights = SampleWeights.*Likelihood.*(Prior./Posterior) ;

src/conditional_particle_filter.m

@@ -68,10 +68,8 @@ end
 
 % Set persistent variables.
 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);
     init_flag = 1;
     number_of_particles = ParticleOptions.number_of_particles ;
@@ -84,25 +82,23 @@ if isempty(H)
     H = 0;
     H_lower_triangular_cholesky = 0;
 else
-    H_lower_triangular_cholesky = chol(H)'; %reduced_rank_cholesky(H)';
+    H_lower_triangular_cholesky = chol(H)'; 
 end
 
 % Get initial condition for the state vector.
 StateVectorMean = ReducedForm.StateVectorMean;
-StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
+StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
 state_variance_rank = size(StateVectorVarianceSquareRoot,2);
-Q_lower_triangular_cholesky = chol(Q)'; %reduced_rank_cholesky(Q)';
+Q_lower_triangular_cholesky = chol(Q)'; 
 
 % Set seed for randn().
 set_dynare_seed('default');
 
 % Initialization of the likelihood.
-normconst2 = log(2*pi)*number_of_observed_variables*prod(diag(H_lower_triangular_cholesky)) ;
+normconst2 = sqrt( ((2*pi)^number_of_observed_variables)*prod(det(H)) ) ;
 lik  = NaN(sample_size,1);
 LIK  = NaN;
-
 ks = 0 ;
-
 StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
 SampleWeights = ones(1,number_of_particles)/number_of_particles ;
 for t=1:sample_size
@@ -120,6 +116,4 @@ for t=1:sample_size
     end
 end
 
-LIK = -sum(lik(start:end));
-
-
+LIK = -sum(lik(start:end));
\ No newline at end of file