diff --git a/nonlinear-filters/src/conditional_filter_proposal.m b/nonlinear-filters/src/conditional_filter_proposal.m index f593ad863..e692c66b5 100644 --- a/nonlinear-filters/src/conditional_filter_proposal.m +++ b/nonlinear-filters/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) ; diff --git a/nonlinear-filters/src/conditional_particle_filter.m b/nonlinear-filters/src/conditional_particle_filter.m index 57ee0f1e9..89c37b954 100644 --- a/nonlinear-filters/src/conditional_particle_filter.m +++ b/nonlinear-filters/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