stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Replaced DynareOptions by ParticleOptions and ThreadsOptions.

rm-particles^2
Stéphane Adjemian (Charybdis) 2014-12-15 15:27:28 +01:00
parent f66b6b6e8c
commit f512e39453
7 changed files with 53 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/gaussian_densities.m
View File

@ -1,4 +1,4 @@
function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,DynareOptions) function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions)
% %
% Elements to calculate the importance sampling ratio % Elements to calculate the importance sampling ratio
% %
@ -41,7 +41,7 @@ proposal = probability2(mut_t,sqr_Pss_t_t,particles) ;
% prior density % prior density
prior = probability2(st_t_1,sqr_Pss_t_t_1,particles) ; prior = probability2(st_t_1,sqr_Pss_t_t_1,particles) ;
% likelihood % likelihood
yt_t_1_i = measurement_equations(particles,ReducedForm,DynareOptions) ; yt_t_1_i = measurement_equations(particles,ReducedForm,ThreadsOptions) ;
eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ; eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
yt_t_1 = sum(yt_t_1_i*weigths1,2) ; yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ; tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;

28
src/gaussian_filter.m
View File

@ -1,4 +1,4 @@
function [LIK,lik] = gaussian_filter(ReducedForm,Y,start,DynareOptions) function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions)
% Evaluates the likelihood of a non-linear model approximating the % Evaluates the likelihood of a non-linear model approximating the
% predictive (prior) and filtered (posterior) densities for state variables % predictive (prior) and filtered (posterior) densities for state variables
% by gaussian distributions. % by gaussian distributions.
@ -67,25 +67,25 @@ if isempty(init_flag)
sample_size = size(Y,2); sample_size = size(Y,2);
number_of_state_variables = length(mf0); number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1); number_of_observed_variables = length(mf1);
number_of_particles = DynareOptions.particle.number_of_particles; number_of_particles = ParticleOptions.number_of_particles;
init_flag = 1; init_flag = 1;
end end
% compute gaussian quadrature nodes and weights on states and shocks % compute gaussian quadrature nodes and weights on states and shocks
if isempty(nodes2) if isempty(nodes2)
if DynareOptions.particle.distribution_approximation.cubature if ParticleOptions.distribution_approximation.cubature
[nodes2,weights2] = spherical_radial_sigma_points(number_of_state_variables); [nodes2,weights2] = spherical_radial_sigma_points(number_of_state_variables);
weights_c2 = weights2; weights_c2 = weights2;
elseif DynareOptions.particle.distribution_approximation.unscented elseif ParticleOptions.distribution_approximation.unscented
[nodes2,weights2,weights_c2] = unscented_sigma_points(number_of_state_variables,DynareOptions); [nodes2,weights2,weights_c2] = unscented_sigma_points(number_of_state_variables,ParticleOptions);
else else
if ~DynareOptions.particle.distribution_approximation.montecarlo if ~ParticleOptions.distribution_approximation.montecarlo
error('Estimation: This approximation for the proposal is not implemented or unknown!') error('Estimation: This approximation for the proposal is not implemented or unknown!')
end end
end end
end end
if DynareOptions.particle.distribution_approximation.montecarlo if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default'); set_dynare_seed('default');
end end
@ -117,16 +117,16 @@ ks = 0 ;
for t=1:sample_size for t=1:sample_size
% build the proposal % build the proposal
[PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = ... [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm,Y(:,t),StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions) ; gaussian_filter_bank(ReducedForm,Y(:,t),StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions) ;
%Estimate(:,t) = PredictedStateMean ; %Estimate(:,t) = PredictedStateMean ;
%V_Estimate(:,:,t) = PredictedStateVarianceSquareRoot ; %V_Estimate(:,:,t) = PredictedStateVarianceSquareRoot ;
if DynareOptions.particle.distribution_approximation.cubature || DynareOptions.particle.distribution_approximation.unscented if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented
StateParticles = bsxfun(@plus,StateVectorMean,StateVectorVarianceSquareRoot*nodes2') ; StateParticles = bsxfun(@plus,StateVectorMean,StateVectorVarianceSquareRoot*nodes2') ;
IncrementalWeights = ... IncrementalWeights = ...
gaussian_densities(Y(:,t),StateVectorMean,... gaussian_densities(Y(:,t),StateVectorMean,...
StateVectorVarianceSquareRoot,PredictedStateMean,... StateVectorVarianceSquareRoot,PredictedStateMean,...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,... PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,...
weights2,weights_c2,ReducedForm,DynareOptions) ; weights2,weights_c2,ReducedForm,ThreadsOptions) ;
SampleWeights = weights2.*IncrementalWeights ; SampleWeights = weights2.*IncrementalWeights ;
SumSampleWeights = sum(SampleWeights) ; SumSampleWeights = sum(SampleWeights) ;
lik(t) = log(SumSampleWeights) ; lik(t) = log(SumSampleWeights) ;
@ -137,7 +137,7 @@ for t=1:sample_size
gaussian_densities(Y(:,t),StateVectorMean,... gaussian_densities(Y(:,t),StateVectorMean,...
StateVectorVarianceSquareRoot,PredictedStateMean,... StateVectorVarianceSquareRoot,PredictedStateMean,...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,... PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,...
1/number_of_particles,1/number_of_particles,ReducedForm,DynareOptions) ; 1/number_of_particles,1/number_of_particles,ReducedForm,ThreadsOptions) ;
SampleWeights = SampleWeights.*IncrementalWeights ; SampleWeights = SampleWeights.*IncrementalWeights ;
SumSampleWeights = sum(SampleWeights) ; SumSampleWeights = sum(SampleWeights) ;
%VarSampleWeights = IncrementalWeights-SumSampleWeights ; %VarSampleWeights = IncrementalWeights-SumSampleWeights ;
@ -145,13 +145,13 @@ for t=1:sample_size
lik(t) = log(SumSampleWeights) ; %+ .5*VarSampleWeights/(number_of_particles*(SumSampleWeights*SumSampleWeights)) ; lik(t) = log(SumSampleWeights) ; %+ .5*VarSampleWeights/(number_of_particles*(SumSampleWeights*SumSampleWeights)) ;
SampleWeights = SampleWeights./SumSampleWeights ; SampleWeights = SampleWeights./SumSampleWeights ;
Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ; Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
if (Neff<.5*sample_size && DynareOptions.particle.resampling.status.generic) || DynareOptions.particle.resampling.status.systematic if (Neff<.5*sample_size && ParticleOptions.resampling.status.generic) || ParticleOptions.resampling.status.systematic
ks = ks + 1 ; ks = ks + 1 ;
StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ; StateParticles = resample(StateParticles',SampleWeights,ParticleOptions)' ;
StateVectorMean = mean(StateParticles,2) ; StateVectorMean = mean(StateParticles,2) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )'; StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )';
SampleWeights = 1/number_of_particles ; SampleWeights = 1/number_of_particles ;
elseif DynareOptions.particle.resampling.status.none elseif ParticleOptions.resampling.status.none
StateVectorMean = (sampleWeights*StateParticles)' ; StateVectorMean = (sampleWeights*StateParticles)' ;
temp = sqrt(SampleWeights').*StateParticles ; temp = sqrt(SampleWeights').*StateParticles ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp'*temp - StateVectorMean*(StateVectorMean') )'; StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp'*temp - StateVectorMean*(StateVectorMean') )';

12
src/gaussian_filter_bank.m
View File

@ -1,4 +1,4 @@
function [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = gaussian_filter_bank(ReducedForm,obs,StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions) function [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = gaussian_filter_bank(ReducedForm,obs,StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions)
% %
% Computes the proposal with a gaussian approximation for importance % Computes the proposal with a gaussian approximation for importance
% sampling % sampling
@ -71,11 +71,11 @@ if isempty(init_flag2)
init_flag2 = 1; init_flag2 = 1;
end 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_state_variables+number_of_structural_innovations) ; [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations) ;
weights_c = weights ; weights_c = weights ;
elseif DynareOptions.particle.proposal_approximation.unscented elseif ParticleOptions.proposal_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,DynareOptions); [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,ParticleOptions);
else else
error('Estimation: This approximation for the proposal is not implemented or unknown!') error('Estimation: This approximation for the proposal is not implemented or unknown!')
end end
@ -87,11 +87,11 @@ sigma_points = bsxfun(@plus,xbar,sqr_Px*(nodes'));
StateVectors = sigma_points(1:number_of_state_variables,:); StateVectors = sigma_points(1:number_of_state_variables,:);
epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:); epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
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 ; PredictedStateMean = tmp(mf0,:)*weights ;
PredictedObservedMean = tmp(mf1,:)*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); PredictedStateMean = sum(PredictedStateMean,2);
PredictedObservedMean = sum(PredictedObservedMean,2); PredictedObservedMean = sum(PredictedObservedMean,2);
dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean)'.*sqrt(weights); dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean)'.*sqrt(weights);

5
src/gaussian_mixture_densities.m
View File

@ -1,7 +1,6 @@
function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
StateMuPost,StateSqrtPPost,StateWeightsPost,... StateMuPost,StateSqrtPPost,StateWeightsPost,...
StateParticles,H,normconst,weigths1,weigths2,ReducedForm,DynareOptions) StateParticles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions)
% %
% Elements to calculate the importance sampling ratio % Elements to calculate the importance sampling ratio
% %
@ -46,7 +45,7 @@ prior = prior' ;
[ras,ras,proposal] = probability(StateMuPost,StateSqrtPPost,StateWeightsPost,StateParticles) ; [ras,ras,proposal] = probability(StateMuPost,StateSqrtPPost,StateWeightsPost,StateParticles) ;
proposal = proposal' ; proposal = proposal' ;
% Compute the density of the current observation conditionally to each particle % Compute the density of the current observation conditionally to each particle
yt_t_1_i = measurement_equations(StateParticles,ReducedForm,DynareOptions) ; yt_t_1_i = measurement_equations(StateParticles,ReducedForm,ThreadsOptions) ;
eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ; eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
yt_t_1 = sum(yt_t_1_i*weigths1,2) ; yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ; tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;

36
src/gaussian_mixture_filter.m
View File

@ -1,4 +1,4 @@
function [LIK,lik] = gaussian_mixture_filter(ReducedForm,Y,start,DynareOptions) function [LIK,lik] = gaussian_mixture_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions)
% Evaluates the likelihood of a non-linear model approximating the state % Evaluates the likelihood of a non-linear model approximating the state
% variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing % variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing
% a wide variety of generalized distributions (when multimodal for instance). % a wide variety of generalized distributions (when multimodal for instance).
@ -71,12 +71,12 @@ if isempty(init_flag)
number_of_state_variables = length(mf0); number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1); number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q); number_of_structural_innovations = length(ReducedForm.Q);
G = DynareOptions.particle.mixture_state_variables; % number of GM components in state G = ParticleOptions.mixture_state_variables; % number of GM components in state
I = DynareOptions.particle.mixture_structural_shocks ; % number of GM components in structural noise I = ParticleOptions.mixture_structural_shocks ; % number of GM components in structural noise
J = DynareOptions.particle.mixture_measurement_shocks ; % number of GM components in observation noise J = ParticleOptions.mixture_measurement_shocks ; % number of GM components in observation noise
Gprime = G*I ; Gprime = G*I ;
Gsecond = G*I*J ; Gsecond = G*I*J ;
number_of_particles = DynareOptions.particle.number_of_particles; number_of_particles = ParticleOptions.number_of_particles;
init_flag = 1; init_flag = 1;
end end
@ -84,19 +84,19 @@ SampleWeights = ones(Gsecond,1)/Gsecond ;
% compute gaussian quadrature nodes and weights on states and shocks % compute gaussian quadrature nodes and weights on states and shocks
if isempty(nodes) if isempty(nodes)
if DynareOptions.particle.distribution_approximation.cubature if ParticleOptions.distribution_approximation.cubature
[nodes,weights] = spherical_radial_sigma_points(number_of_state_variables); [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables);
weights_c = weights; weights_c = weights;
elseif DynareOptions.particle.distribution_approximation.unscented elseif ParticleOptions.distribution_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables,DynareOptions); [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables,ParticleOptions);
else else
if ~DynareOptions.particle.distribution_approximation.montecarlo if ~ParticleOptions.distribution_approximation.montecarlo
error('Estimation: This approximation for the proposal is not implemented or unknown!') error('Estimation: This approximation for the proposal is not implemented or unknown!')
end end
end end
end end
if DynareOptions.particle.distribution_approximation.montecarlo if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default'); set_dynare_seed('default');
SampleWeights = 1/number_of_particles ; SampleWeights = 1/number_of_particles ;
end end
@ -161,7 +161,7 @@ for t=1:sample_size
gaussian_mixture_filter_bank(ReducedForm,Y(:,t),StateMu(:,g),StateSqrtP(:,:,g),StateWeights(1,g),... gaussian_mixture_filter_bank(ReducedForm,Y(:,t),StateMu(:,g),StateSqrtP(:,:,g),StateWeights(1,g),...
StructuralShocksMu(:,i),StructuralShocksSqrtP(:,:,i),StructuralShocksWeights(1,i),... StructuralShocksMu(:,i),StructuralShocksSqrtP(:,:,i),StructuralShocksWeights(1,i),...
ObservationShocksMu(:,j),ObservationShocksSqrtP(:,:,j),ObservationShocksWeights(1,j),... ObservationShocksMu(:,j),ObservationShocksSqrtP(:,:,j),ObservationShocksWeights(1,j),...
H,H_lower_triangular_cholesky,const_lik,DynareOptions) ; H,H_lower_triangular_cholesky,const_lik,ParticleOptions,ThreadsOptions) ;
end end
end end
end end
@ -170,12 +170,12 @@ for t=1:sample_size
StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior,2) ; StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior,2) ;
StateWeightsPost = StateWeightsPost/sum(StateWeightsPost,2) ; StateWeightsPost = StateWeightsPost/sum(StateWeightsPost,2) ;
if DynareOptions.particle.distribution_approximation.cubature || DynareOptions.particle.distribution_approximation.unscented if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented
for i=1:Gsecond for i=1:Gsecond
StateParticles = bsxfun(@plus,StateMuPost(:,i),StateSqrtPPost(:,:,i)*nodes') ; StateParticles = bsxfun(@plus,StateMuPost(:,i),StateSqrtPPost(:,:,i)*nodes') ;
IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
StateMuPost,StateSqrtPPost,StateWeightsPost,... StateMuPost,StateSqrtPPost,StateWeightsPost,...
StateParticles,H,const_lik,weights,weights_c,ReducedForm,DynareOptions) ; StateParticles,H,const_lik,weights,weights_c,ReducedForm,ThreadsOptions) ;
SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights) ; SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights) ;
end end
SumSampleWeights = sum(SampleWeights) ; SumSampleWeights = sum(SampleWeights) ;
@ -183,7 +183,7 @@ for t=1:sample_size
SampleWeights = SampleWeights./SumSampleWeights ; SampleWeights = SampleWeights./SumSampleWeights ;
[ras,SortedRandomIndx] = sort(rand(1,Gsecond)); [ras,SortedRandomIndx] = sort(rand(1,Gsecond));
SortedRandomIndx = SortedRandomIndx(1:G); SortedRandomIndx = SortedRandomIndx(1:G);
indx = index_resample(0,SampleWeights,DynareOptions) ; indx = index_resample(0,SampleWeights,ParticleOptions) ;
indx = indx(SortedRandomIndx) ; indx = indx(SortedRandomIndx) ;
StateMu = StateMuPost(:,indx); StateMu = StateMuPost(:,indx);
StateSqrtP = StateSqrtPPost(:,:,indx); StateSqrtP = StateSqrtPPost(:,:,indx);
@ -195,7 +195,7 @@ for t=1:sample_size
IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
StateMuPost,StateSqrtPPost,StateWeightsPost,... StateMuPost,StateSqrtPPost,StateWeightsPost,...
StateParticles,H,const_lik,1/number_of_particles,... StateParticles,H,const_lik,1/number_of_particles,...
1/number_of_particles,ReducedForm,DynareOptions) ; 1/number_of_particles,ReducedForm,ThreadsOptions) ;
% calculate importance weights of particles % calculate importance weights of particles
SampleWeights = SampleWeights.*IncrementalWeights ; SampleWeights = SampleWeights.*IncrementalWeights ;
SumSampleWeights = sum(SampleWeights,1) ; SumSampleWeights = sum(SampleWeights,1) ;
@ -205,13 +205,13 @@ for t=1:sample_size
%estimate(t,:,1) = SampleWeights*StateParticles' ; %estimate(t,:,1) = SampleWeights*StateParticles' ;
% Resampling if needed of required % Resampling if needed of required
Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ; Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
if (DynareOptions.particle.resampling.status.generic && Neff<.5*sample_size) || DynareOptions.particle.resampling.status.systematic if (ParticleOptions.resampling.status.generic && Neff<.5*sample_size) || ParticleOptions.resampling.status.systematic
ks = ks + 1 ; ks = ks + 1 ;
StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ; StateParticles = resample(StateParticles',SampleWeights,ParticleOptions)' ;
StateVectorMean = mean(StateParticles,2) ; StateVectorMean = mean(StateParticles,2) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/number_of_particles - StateVectorMean*(StateVectorMean') )'; StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/number_of_particles - StateVectorMean*(StateVectorMean') )';
SampleWeights = 1/number_of_particles ; SampleWeights = 1/number_of_particles ;
elseif DynareOptions.particle.resampling.status.none elseif ParticleOptions.resampling.status.none
StateVectorMean = StateParticles*sampleWeights ; StateVectorMean = StateParticles*sampleWeights ;
temp = sqrt(SampleWeights').*StateParticles ; temp = sqrt(SampleWeights').*StateParticles ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )'; StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )';

12
src/gaussian_mixture_filter_bank.m
View File

@ -2,7 +2,7 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPP
gaussian_mixture_filter_bank(ReducedForm,obs,StateMu,StateSqrtP,StateWeights,... gaussian_mixture_filter_bank(ReducedForm,obs,StateMu,StateSqrtP,StateWeights,...
StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,... StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,...
ObservationShocksMu,ObservationShocksSqrtP,ObservationShocksWeights,... ObservationShocksMu,ObservationShocksSqrtP,ObservationShocksWeights,...
H,H_lower_triangular_cholesky,normfactO,DynareOptions) H,H_lower_triangular_cholesky,normfactO,ParticleOptions,ThreadsOptions)
% %
% Computes the proposal with a gaussian approximation for importance % Computes the proposal with a gaussian approximation for importance
% sampling % sampling
@ -79,11 +79,11 @@ end
numb = number_of_state_variables+number_of_structural_innovations ; numb = number_of_state_variables+number_of_structural_innovations ;
if DynareOptions.particle.proposal_approximation.cubature if ParticleOptions.proposal_approximation.cubature
[nodes3,weights3] = spherical_radial_sigma_points(numb); [nodes3,weights3] = spherical_radial_sigma_points(numb);
weights_c3 = weights3; weights_c3 = weights3;
elseif DynareOptions.particle.proposal_approximation.unscented elseif ParticleOptions.proposal_approximation.unscented
[nodes3,weights3,weights_c3] = unscented_sigma_points(numb,DynareOptions); [nodes3,weights3,weights_c3] = unscented_sigma_points(numb,ParticleOptions);
else else
error('Estimation: This approximation for the proposal is not implemented or unknown!') error('Estimation: This approximation for the proposal is not implemented or unknown!')
end end
@ -91,11 +91,11 @@ end
epsilon = bsxfun(@plus,StructuralShocksSqrtP*nodes3(:,number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations)',StructuralShocksMu) ; epsilon = bsxfun(@plus,StructuralShocksSqrtP*nodes3(:,number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations)',StructuralShocksMu) ;
StateVectors = bsxfun(@plus,StateSqrtP*nodes3(:,1:number_of_state_variables)',StateMu); StateVectors = bsxfun(@plus,StateSqrtP*nodes3(:,1:number_of_state_variables)',StateMu);
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
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,:)*weights3; PredictedStateMean = tmp(mf0,:)*weights3;
PredictedObservedMean = tmp(mf1,:)*weights3; PredictedObservedMean = tmp(mf1,:)*weights3;
if DynareOptions.particle.proposal_approximation.cubature if ParticleOptions.proposal_approximation.cubature
PredictedStateMean = sum(PredictedStateMean,2); PredictedStateMean = sum(PredictedStateMean,2);
PredictedObservedMean = sum(PredictedObservedMean,2); PredictedObservedMean = sum(PredictedObservedMean,2);
dState = (bsxfun(@minus,tmp(mf0,:),PredictedStateMean)').*sqrt(weights3); dState = (bsxfun(@minus,tmp(mf0,:),PredictedStateMean)').*sqrt(weights3);

10
src/index_resample.m
View File

@ -1,4 +1,4 @@
function resampling_index = index_resample(particles,weights,DynareOptions) function resampling_index = index_resample(particles,weights,ParticleOptions)
% Resamples particles. % Resamples particles.
%@info: %@info:
@ -54,17 +54,17 @@ function resampling_index = index_resample(particles,weights,DynareOptions)
defaultmethod = 1; % For residual based method set this variable equal to 0. defaultmethod = 1; % For residual based method set this variable equal to 0.
if defaultmethod if defaultmethod
if DynareOptions.particle.resampling.method.kitagawa if ParticleOptions.resampling.method.kitagawa
resampling_index = traditional_resampling(particles,weights,rand); resampling_index = traditional_resampling(particles,weights,rand);
elseif DynareOptions.particle.resampling.method.stratified elseif ParticleOptions.resampling.method.stratified
resampling_index = traditional_resampling(particles,weights,rand(size(weights))); resampling_index = traditional_resampling(particles,weights,rand(size(weights)));
else else
error('Unknow sampling method!') error('Unknow sampling method!')
end end
else else
if DynareOptions.particle.resampling.method.kitagawa if ParticleOptions.resampling.method.kitagawa
resampled_particles = residual_resampling(particles,weights,rand); resampled_particles = residual_resampling(particles,weights,rand);
elseif DynareOptions.particle.resampling.method.stratified elseif ParticleOptions.resampling.method.stratified
resampled_particles = residual_resampling(particles,weights,rand(size(weights))); resampled_particles = residual_resampling(particles,weights,rand(size(weights)));
else else
error('Unknown sampling method!') error('Unknown sampling method!')