Correction on a bug and extension of the original paper for the calculation of incremental weights.
parent
3ffcb7b759
commit
35624147df
|
@ -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) ;
|
||||
|
|
|
@ -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));
|
Loading…
Reference in New Issue