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

factorize some codes across options.

remove-submodule^2
Frédéric Karamé 2015-01-23 10:05:16 +01:00
parent 9f6a300cf0
commit ffd62a5923
1 changed files with 18 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
nonlinear-filters/src/gaussian_mixture_filter.m
View File

@ -3,7 +3,6 @@ function [LIK,lik] = gaussian_mixture_filter(ReducedForm,Y,start,ParticleOptions
% variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing
% a wide variety of generalized distributions (when multimodal for instance).
% Each gaussian distribution is obtained whether
% - with a Smolyak quadrature à la Kronrod & Paterson (Heiss & Winschel 2010, Winschel & Kratzig 2010).
% - with a radial-spherical cubature
% - with scaled unscented sigma-points
% A Sparse grid Kalman Filter is implemented on each component of the mixture,
@ -72,8 +71,8 @@ if isempty(init_flag)
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
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
I = 1 ; %ParticleOptions.mixture_structural_shocks ; % number of GM components in structural noise
J = 1 ; %ParticleOptions.mixture_measurement_shocks ; % number of GM components in observation noise
Gprime = G*I ;
Gsecond = G*I*J ;
number_of_particles = ParticleOptions.number_of_particles;
@ -116,8 +115,9 @@ Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)';
StateWeights = ones(1,G)/G ;
StateMu = ReducedForm.StateVectorMean*ones(1,G) ;
StateSqrtP = zeros(number_of_state_variables,number_of_state_variables,G) ;
temp = reduced_rank_cholesky(ReducedForm.StateVectorVariance)' ;
for g=1:G
StateSqrtP(:,:,g) = reduced_rank_cholesky(ReducedForm.StateVectorVariance)' ;
StateSqrtP(:,:,g) = temp ;
end
StructuralShocksWeights = ones(1,I)/I ;
@ -197,25 +197,26 @@ for t=1:sample_size
StateParticles,H,const_lik,1/number_of_particles,...
1/number_of_particles,ReducedForm,ThreadsOptions) ;
% calculate importance weights of particles
SampleWeights = SampleWeights.*IncrementalWeights ;
% SampleWeights = SampleWeights.*IncrementalWeights ;
SampleWeights = IncrementalWeights/number_of_particles ;
SumSampleWeights = sum(SampleWeights,1) ;
SampleWeights = SampleWeights./SumSampleWeights ;
lik(t) = log(SumSampleWeights) ;
% First possible state point estimates
%estimate(t,:,1) = SampleWeights*StateParticles' ;
% Resampling if needed of required
Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
if (ParticleOptions.resampling.status.generic && Neff<.5*sample_size) || ParticleOptions.resampling.status.systematic
ks = ks + 1 ;
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 ParticleOptions.resampling.status.none
StateVectorMean = StateParticles*sampleWeights ;
temp = sqrt(SampleWeights').*StateParticles ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )';
end
% Neff = neff(SampleWeights) ;
% if (ParticleOptions.resampling.status.generic && Neff<.5*sample_size) || ParticleOptions.resampling.status.systematic
% ks = ks + 1 ;
% 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 ParticleOptions.resampling.status.none
% StateVectorMean = StateParticles*sampleWeights ;
% temp = sqrt(SampleWeights').*StateParticles ;
% StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )';
% end
% Use the information from particles to update the gaussian mixture on state variables
[StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(StateParticles,StateMu,StateSqrtP,StateWeights,0.001,10,1) ;
%estimate(t,:,3) = StateWeights*StateMu' ;