diff --git a/matlab/particle/sequential_importance_particle_filter.m b/matlab/particle/sequential_importance_particle_filter.m
index ee13ec2e0..129f54db5 100644
--- a/matlab/particle/sequential_importance_particle_filter.m
+++ b/matlab/particle/sequential_importance_particle_filter.m
@@ -121,10 +121,10 @@ 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);
     tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
-    PredictedObservedMean = mean(tmp(mf1,:),2);
+    PredictedObservedMean = tmp(mf1,:)*transpose(weights);
     PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
     dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
-    PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles+H;
+    PredictedObservedVariance = dPredictedObservedMean*diag(weights)*dPredictedObservedMean' + H;
     lnw = -.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1));
     dfac = max(lnw);
     wtilde = weights.*exp(lnw-dfac);