Commit f66b6b6e authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Replaced DynareOptions by ParticleOptions and ThreadsOptions.

parent 7ad84f3f
function [ProposalStateVector,Weights] = conditional_filter_proposal(ReducedForm,obs,StateVectors,SampleWeights,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions,normconst2)
function [ProposalStateVector,Weights] = conditional_filter_proposal(ReducedForm,obs,StateVectors,SampleWeights,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2)
%
% Computes the proposal for each past particle using Gaussian approximations
% for the state errors and the Kalman filter
......@@ -73,11 +73,11 @@ if isempty(init_flag2)
init_flag2 = 1;
end
if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.proposal_approximation.montecarlo
if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_approximation.montecarlo
[nodes,weights] = spherical_radial_sigma_points(number_of_structural_innovations);
weights_c = weights ;
elseif DynareOptions.particle.proposal_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_structural_innovations,DynareOptions);
elseif ParticleOptions.proposal_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_structural_innovations,ParticleOptions);
else
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
......@@ -85,12 +85,12 @@ end
epsilon = Q_lower_triangular_cholesky*(nodes') ;
yhat = repmat(StateVectors-state_variables_steady_state,1,size(epsilon,2)) ;
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
PredictedStateMean = tmp(mf0,:)*weights ;
PredictedObservedMean = tmp(mf1,:)*weights;
if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.proposal_approximation.montecarlo
if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_approximation.montecarlo
PredictedStateMean = sum(PredictedStateMean,2) ;
PredictedObservedMean = sum(PredictedObservedMean,2) ;
dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean)'.*sqrt(weights) ;
......@@ -117,7 +117,7 @@ else
end
ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ;
ypred = measurement_equations(ProposalStateVector,ReducedForm,DynareOptions) ;
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;
function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,DynareOptions)
function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions)
%
% Evaluates the likelihood of a non-linear model with a particle filter
% - the proposal is built using the Kalman updating step for each particle.
......@@ -74,7 +74,7 @@ if isempty(init_flag)
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
init_flag = 1;
number_of_particles = DynareOptions.particle.number_of_particles ;
number_of_particles = ParticleOptions.number_of_particles ;
end
% Get covariance matrices
......@@ -108,14 +108,14 @@ SampleWeights = ones(1,number_of_particles)/number_of_particles ;
for t=1:sample_size
for i=1:number_of_particles
[StateParticles(:,i),SampleWeights(i)] = ...
conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions,normconst2) ;
conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2) ;
end
SumSampleWeights = sum(SampleWeights) ;
lik(t) = log(SumSampleWeights) ;
SampleWeights = SampleWeights./SumSampleWeights ;
if (DynareOptions.particle.resampling.status.generic && neff(SampleWeights)<DynareOptions.particle.resampling.threshold*sample_size) || DynareOptions.particle.resampling.status.systematic
if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
ks = ks + 1 ;
StateParticles = resample(StateParticles',SampleWeights',DynareOptions)';
StateParticles = resample(StateParticles',SampleWeights',ParticleOptions)';
SampleWeights = ones(1,number_of_particles)/number_of_particles ;
end
end
......
function measure = measurement_equations(StateVectors,ReducedForm,DynareOptions)
function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions)
% Copyright (C) 2013 Dynare Team
%
......@@ -27,7 +27,7 @@ constant = ReducedForm.constant(mf1,:);
state_variables_steady_state = ReducedForm.state_variables_steady_state;
number_of_structural_innovations = length(ReducedForm.Q);
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state) ;
measure = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,size(yhat,2)),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
measure = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,size(yhat,2)),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
function [nodes,W_m,W_c] = unscented_sigma_points(n,DynareOptions)
function [nodes,W_m,W_c] = unscented_sigma_points(n,ParticleOptions)
%
% Computes nodes and weigths for a scaled unscented transform cubature
% INPUTS
......@@ -31,10 +31,10 @@ function [nodes,W_m,W_c] = unscented_sigma_points(n,DynareOptions)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
lambda = (DynareOptions.particle.unscented.alpha^2)*(n+DynareOptions.particle.unscented.kappa) - n ;
lambda = (ParticleOptions.unscented.alpha^2)*(n+ParticleOptions.unscented.kappa) - n ;
nodes = [ zeros(n,1) ( sqrt(n+lambda).*([ eye(n) -eye(n)]) ) ]' ;
W_m = lambda/(n+lambda) ;
W_c = W_m + (1-DynareOptions.particle.unscented.alpha^2+DynareOptions.particle.unscented.beta) ;
W_c = W_m + (1-ParticleOptions.unscented.alpha^2+ParticleOptions.unscented.beta) ;
temp = ones(2*n,1)/(2*(n+lambda)) ;
W_m = [W_m ; temp] ;
W_c = [W_c ; temp] ;
......
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