Commit 93039cd9 authored by Frédéric Karamé's avatar Frédéric Karamé

fix a bug on format matrix.

parent 04ad104b
......@@ -67,7 +67,7 @@ end
% Get initial conditions for the state particles
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
......@@ -81,9 +81,9 @@ small_a = sqrt(1-h_square) ;
% Initialization of parameter particles
xparam = zeros(number_of_parameters,number_of_particles) ;
stderr = sqrt(bsxfun(@power,bounds.ub+bounds.lb,2)/12)/100 ;
stderr = sqrt(bsxfun(@power,bounds.ub+bounds.lb,2)/12)/50 ;
stderr = sqrt(bsxfun(@power,bounds.ub+bounds.lb,2)/12)/20 ;
stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/100 ;
stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/50 ;
%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/20 ;
i = 1 ;
while i<=number_of_particles
%candidate = start_param + .001*liu_west_chol_sigma_bar*randn(number_of_parameters,1) ;
......@@ -122,7 +122,7 @@ for t=1:sample_size
chol_sigma_bar = chol(h_square*sigma_bar)' ;
end
% Prediction (without shocks)
ObservedVariables = zeros(number_of_observed_variables,number_of_particles) ;
wtilde = zeros(1,number_of_particles) ;
for i=1:number_of_particles
% model resolution
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
......@@ -144,13 +144,9 @@ for t=1:sample_size
else
tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,1),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
end
ObservedVariables(:,i) = tmp(mf1,:) ;
PredictionError = bsxfun(@minus,Y(t,:)',tmp(mf1,:));
wtilde(i) = exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError),1))) ;
end
PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ;
PredictionError = bsxfun(@minus,Y(t,:)',ObservedVariables);
dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean);
PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ;
wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
% unormalized weights and observation likelihood contribution
tau_tilde = weights.*wtilde ;
sum_tau_tilde = sum(tau_tilde) ;
......@@ -161,10 +157,10 @@ for t=1:sample_size
if pruning
StateVectors_ = StateVectors_(:,indx) ;
end
xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ;
xparam = xparam(:,indx) ;
xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam(:,indx)) ;
wtilde = wtilde(indx) ;
% draw in the new distributions
lnw = zeros(1,number_of_particles) ;
i = 1 ;
while i<=number_of_particles
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters,1) ;
......@@ -194,15 +190,11 @@ for t=1:sample_size
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
end
StateVectors(:,i) = tmp(mf0,:) ;
ObservedVariables(:,i) = tmp(mf1,:) ;
PredictionError = bsxfun(@minus,Y(t,:)',tmp(mf1,:));
lnw(i) = exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError),1)));
i = i+1 ;
end
end
PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ;
PredictionError = bsxfun(@minus,Y(t,:)',ObservedVariables);
dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean);
PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ;
lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1)));
% importance ratio
wtilde = lnw./wtilde ;
% normalization
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment