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

Replaced DynareOptions by ParticleOptions and ThreadsOptions.

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

4
nonlinear-filters/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
%
@ -41,7 +41,7 @@ proposal = probability2(mut_t,sqr_Pss_t_t,particles) ;
% prior density
prior = probability2(st_t_1,sqr_Pss_t_t_1,particles) ;
% 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)' ;
yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;

28
nonlinear-filters/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
% predictive (prior) and filtered (posterior) densities for state variables
% by gaussian distributions.
@ -67,25 +67,25 @@ if isempty(init_flag)
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_particles = DynareOptions.particle.number_of_particles;
number_of_particles = ParticleOptions.number_of_particles;
init_flag = 1;
end
% compute gaussian quadrature nodes and weights on states and shocks
if isempty(nodes2)
if DynareOptions.particle.distribution_approximation.cubature
if ParticleOptions.distribution_approximation.cubature
[nodes2,weights2] = spherical_radial_sigma_points(number_of_state_variables);
weights_c2 = weights2;
elseif DynareOptions.particle.distribution_approximation.unscented
[nodes2,weights2,weights_c2] = unscented_sigma_points(number_of_state_variables,DynareOptions);
elseif ParticleOptions.distribution_approximation.unscented
[nodes2,weights2,weights_c2] = unscented_sigma_points(number_of_state_variables,ParticleOptions);
else
if ~DynareOptions.particle.distribution_approximation.montecarlo
if ~ParticleOptions.distribution_approximation.montecarlo
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
end
end
if DynareOptions.particle.distribution_approximation.montecarlo
if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default');
end
@ -117,16 +117,16 @@ ks = 0 ;
for t=1:sample_size
% build the proposal
[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 ;
%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') ;
IncrementalWeights = ...
gaussian_densities(Y(:,t),StateVectorMean,...
StateVectorVarianceSquareRoot,PredictedStateMean,...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,...
weights2,weights_c2,ReducedForm,DynareOptions) ;
weights2,weights_c2,ReducedForm,ThreadsOptions) ;
SampleWeights = weights2.*IncrementalWeights ;
SumSampleWeights = sum(SampleWeights) ;
lik(t) = log(SumSampleWeights) ;
@ -137,7 +137,7 @@ for t=1:sample_size
gaussian_densities(Y(:,t),StateVectorMean,...
StateVectorVarianceSquareRoot,PredictedStateMean,...
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 ;
SumSampleWeights = sum(SampleWeights) ;
%VarSampleWeights = IncrementalWeights-SumSampleWeights ;
@ -145,13 +145,13 @@ for t=1:sample_size
lik(t) = log(SumSampleWeights) ; %+ .5*VarSampleWeights/(number_of_particles*(SumSampleWeights*SumSampleWeights)) ;
SampleWeights = SampleWeights./SumSampleWeights ;
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 ;
StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ;
StateParticles = resample(StateParticles',SampleWeights,ParticleOptions)' ;
StateVectorMean = mean(StateParticles,2) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )';
SampleWeights = 1/number_of_particles ;
elseif DynareOptions.particle.resampling.status.none
elseif ParticleOptions.resampling.status.none
StateVectorMean = (sampleWeights*StateParticles)' ;
temp = sqrt(SampleWeights').*StateParticles ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp'*temp - StateVectorMean*(StateVectorMean') )';

12
nonlinear-filters/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
% sampling
@ -71,11 +71,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_state_variables+number_of_structural_innovations) ;
weights_c = weights ;
elseif DynareOptions.particle.proposal_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,DynareOptions);
elseif ParticleOptions.proposal_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,ParticleOptions);
else
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
@ -87,11 +87,11 @@ sigma_points = bsxfun(@plus,xbar,sqr_Px*(nodes'));
StateVectors = sigma_points(1:number_of_state_variables,:);
epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
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 ;
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);

5
nonlinear-filters/src/gaussian_mixture_densities.m
View File

@ -1,7 +1,6 @@
function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
StateMuPost,StateSqrtPPost,StateWeightsPost,...
StateParticles,H,normconst,weigths1,weigths2,ReducedForm,DynareOptions)
StateParticles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions)
%
% Elements to calculate the importance sampling ratio
%
@ -46,7 +45,7 @@ prior = prior' ;
[ras,ras,proposal] = probability(StateMuPost,StateSqrtPPost,StateWeightsPost,StateParticles) ;
proposal = proposal' ;
% 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)' ;
yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;

36
nonlinear-filters/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
% variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing
% 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_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
G = DynareOptions.particle.mixture_state_variables; % number of GM components in state
I = DynareOptions.particle.mixture_structural_shocks ; % number of GM components in structural noise
J = DynareOptions.particle.mixture_measurement_shocks ; % number of GM components in observation noise
G = ParticleOptions.mixture_state_variables; % number of GM components in state
I = ParticleOptions.mixture_structural_shocks ; % number of GM components in structural noise
J = ParticleOptions.mixture_measurement_shocks ; % number of GM components in observation noise
Gprime = G*I ;
Gsecond = G*I*J ;
number_of_particles = DynareOptions.particle.number_of_particles;
number_of_particles = ParticleOptions.number_of_particles;
init_flag = 1;
end
@ -84,19 +84,19 @@ SampleWeights = ones(Gsecond,1)/Gsecond ;
% compute gaussian quadrature nodes and weights on states and shocks
if isempty(nodes)
if DynareOptions.particle.distribution_approximation.cubature
if ParticleOptions.distribution_approximation.cubature
[nodes,weights] = spherical_radial_sigma_points(number_of_state_variables);
weights_c = weights;
elseif DynareOptions.particle.distribution_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables,DynareOptions);
elseif ParticleOptions.distribution_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables,ParticleOptions);
else
if ~DynareOptions.particle.distribution_approximation.montecarlo
if ~ParticleOptions.distribution_approximation.montecarlo
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
end
end
if DynareOptions.particle.distribution_approximation.montecarlo
if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default');
SampleWeights = 1/number_of_particles ;
end
@ -161,7 +161,7 @@ for t=1:sample_size
gaussian_mixture_filter_bank(ReducedForm,Y(:,t),StateMu(:,g),StateSqrtP(:,:,g),StateWeights(1,g),...
StructuralShocksMu(:,i),StructuralShocksSqrtP(:,:,i),StructuralShocksWeights(1,i),...
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
@ -170,12 +170,12 @@ for t=1:sample_size
StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior,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
StateParticles = bsxfun(@plus,StateMuPost(:,i),StateSqrtPPost(:,:,i)*nodes') ;
IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
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) ;
end
SumSampleWeights = sum(SampleWeights) ;
@ -183,7 +183,7 @@ for t=1:sample_size
SampleWeights = SampleWeights./SumSampleWeights ;
[ras,SortedRandomIndx] = sort(rand(1,Gsecond));
SortedRandomIndx = SortedRandomIndx(1:G);
indx = index_resample(0,SampleWeights,DynareOptions) ;
indx = index_resample(0,SampleWeights,ParticleOptions) ;
indx = indx(SortedRandomIndx) ;
StateMu = StateMuPost(:,indx);
StateSqrtP = StateSqrtPPost(:,:,indx);
@ -195,7 +195,7 @@ for t=1:sample_size
IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
StateMuPost,StateSqrtPPost,StateWeightsPost,...
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
SampleWeights = SampleWeights.*IncrementalWeights ;
SumSampleWeights = sum(SampleWeights,1) ;
@ -205,13 +205,13 @@ for t=1:sample_size
%estimate(t,:,1) = SampleWeights*StateParticles' ;
% Resampling if needed of required
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 ;
StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ;
StateParticles = resample(StateParticles',SampleWeights,ParticleOptions)' ;
StateVectorMean = mean(StateParticles,2) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/number_of_particles - StateVectorMean*(StateVectorMean') )';
SampleWeights = 1/number_of_particles ;
elseif DynareOptions.particle.resampling.status.none
elseif ParticleOptions.resampling.status.none
StateVectorMean = StateParticles*sampleWeights ;
temp = sqrt(SampleWeights').*StateParticles ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )';

12
nonlinear-filters/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,...
StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,...
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
% sampling
@ -79,11 +79,11 @@ end
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);
weights_c3 = weights3;
elseif DynareOptions.particle.proposal_approximation.unscented
[nodes3,weights3,weights_c3] = unscented_sigma_points(numb,DynareOptions);
elseif ParticleOptions.proposal_approximation.unscented
[nodes3,weights3,weights_c3] = unscented_sigma_points(numb,ParticleOptions);
else
error('Estimation: This approximation for the proposal is not implemented or unknown!')
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) ;
StateVectors = bsxfun(@plus,StateSqrtP*nodes3(:,1:number_of_state_variables)',StateMu);
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;
PredictedObservedMean = tmp(mf1,:)*weights3;
if DynareOptions.particle.proposal_approximation.cubature
if ParticleOptions.proposal_approximation.cubature
PredictedStateMean = sum(PredictedStateMean,2);
PredictedObservedMean = sum(PredictedObservedMean,2);
dState = (bsxfun(@minus,tmp(mf0,:),PredictedStateMean)').*sqrt(weights3);

10
nonlinear-filters/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.
%@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.
if defaultmethod
if DynareOptions.particle.resampling.method.kitagawa
if ParticleOptions.resampling.method.kitagawa
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)));
else
error('Unknow sampling method!')
end
else
if DynareOptions.particle.resampling.method.kitagawa
if ParticleOptions.resampling.method.kitagawa
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)));
else
error('Unknown sampling method!')