Allow k order approximation in Gaussian Filter (gf).

Ref. dynare#1673
rm-particles^2
Stéphane Adjemian (Charybdis) 2019-12-21 12:06:41 +01:00
parent 89a941605f
commit 472d755d98
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
4 changed files with 115 additions and 125 deletions

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,ThreadsOptions)
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,DynareOptions, Model)
%
% Elements to calculate the importance sampling ratio
%
@ -19,7 +19,8 @@ function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sq
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2017 Dynare Team
% Copyright (C) 2009-2019 Dynare Team
%
% This file is part of Dynare.
%
@ -37,17 +38,13 @@ function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sq
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% proposal density
proposal = probability2(mut_t,sqr_Pss_t_t,particles) ;
proposal = probability2(mut_t, sqr_Pss_t_t, particles);
% 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
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) ;
%Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ;
%sqr_det = sqrt(det(Pyy)) ;
%foo = (eta_t_i/Pyy).*eta_t_i ;
%likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ;
likelihood = probability2(obs,sqrt(H),yt_t_1_i) ;
IncrementalWeights = likelihood.*prior./proposal ;
yt_t_1_i = measurement_equations(particles, ReducedForm, ThreadsOptions, DynareOptions, Model);
likelihood = probability2(obs, sqrt(H), yt_t_1_i);
IncrementalWeights = likelihood.*prior./proposal;

View File

@ -1,4 +1,5 @@
function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions)
function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model)
% Evaluates the likelihood of a non-linear model approximating the
% predictive (prior) and filtered (posterior) densities for state variables
% by gaussian distributions.
@ -30,7 +31,8 @@ function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, Thr
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2017 Dynare Team
% Copyright (C) 2009-2019 Dynare Team
%
% This file is part of Dynare.
%
@ -47,37 +49,27 @@ function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, Thr
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent init_flag mf0 mf1
persistent nodes2 weights2 weights_c2 number_of_particles
persistent sample_size number_of_state_variables number_of_observed_variables
% Set default
if isempty(start)
start = 1;
end
% Set persistent variables.
if isempty(init_flag)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_particles = ParticleOptions.number_of_particles;
init_flag = 1;
end
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_particles = ParticleOptions.number_of_particles;
% compute gaussian quadrature nodes and weights on states and shocks
if isempty(nodes2)
if ParticleOptions.distribution_approximation.cubature
[nodes2,weights2] = spherical_radial_sigma_points(number_of_state_variables);
weights_c2 = weights2;
elseif ParticleOptions.distribution_approximation.unscented
[nodes2,weights2,weights_c2] = unscented_sigma_points(number_of_state_variables,ParticleOptions);
else
if ~ParticleOptions.distribution_approximation.montecarlo
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
if ParticleOptions.distribution_approximation.cubature
[nodes2, weights2] = spherical_radial_sigma_points(number_of_state_variables);
weights_c2 = weights2;
elseif ParticleOptions.distribution_approximation.unscented
[nodes2, weights2, weights_c2] = unscented_sigma_points(number_of_state_variables,ParticleOptions);
else
if ~ParticleOptions.distribution_approximation.montecarlo
error('This approximation for the proposal is unknown!')
end
end
@ -107,38 +99,37 @@ lik = NaN(sample_size,1);
LIK = NaN;
for t=1:sample_size
[PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm,Y(:,t),StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions) ;
[PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm, Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ...
H, ParticleOptions, ThreadsOptions, DynareOptions, Model);
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,ThreadsOptions) ;
SampleWeights = weights2.*IncrementalWeights ;
StateParticles = bsxfun(@plus, StateVectorMean, StateVectorVarianceSquareRoot*nodes2');
IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ...
PredictedStateVarianceSquareRoot, StateParticles, H, const_lik, ...
weights2, weights_c2, ReducedForm, ThreadsOptions, ...
DynareOptions, Model);
SampleWeights = weights2.*IncrementalWeights;
else
StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean) ;
IncrementalWeights = ...
gaussian_densities(Y(:,t),StateVectorMean,...
StateVectorVarianceSquareRoot,PredictedStateMean,...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,...
1/number_of_particles,1/number_of_particles,ReducedForm,ThreadsOptions) ;
SampleWeights = IncrementalWeights/number_of_particles ;
StateParticles = bsxfun(@plus, StateVectorVarianceSquareRoot*randn(state_variance_rank, number_of_particles), StateVectorMean) ;
IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik, ...
1/number_of_particles,1/number_of_particles,ReducedForm,ThreadsOptions, ...
DynareOptions, Model);
SampleWeights = IncrementalWeights/number_of_particles;
end
SampleWeights = SampleWeights + 1e-6*ones(size(SampleWeights,1),1) ;
SumSampleWeights = sum(SampleWeights) ;
lik(t) = log(SumSampleWeights) ;
SampleWeights = SampleWeights./SumSampleWeights ;
SampleWeights = SampleWeights + 1e-6*ones(size(SampleWeights, 1), 1);
SumSampleWeights = sum(SampleWeights);
lik(t) = log(SumSampleWeights);
SampleWeights = SampleWeights./SumSampleWeights;
if not(ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented)
if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
StateParticles = resample(StateParticles',SampleWeights,ParticleOptions)' ;
SampleWeights = ones(number_of_particles,1)/number_of_particles;
StateParticles = resample(StateParticles', SampleWeights, ParticleOptions)';
SampleWeights = ones(number_of_particles, 1)/number_of_particles;
end
end
StateVectorMean = StateParticles*SampleWeights ;
temp = bsxfun(@minus,StateParticles,StateVectorMean) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( bsxfun(@times,SampleWeights',temp)*temp' )';
StateVectorMean = StateParticles*SampleWeights;
temp = bsxfun(@minus, StateParticles, StateVectorMean);
StateVectorVarianceSquareRoot = reduced_rank_cholesky(bsxfun(@times,SampleWeights',temp)*temp')';
end
LIK = -sum(lik(start:end));
LIK = -sum(lik(start:end));

View File

@ -1,4 +1,6 @@
function [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = gaussian_filter_bank(ReducedForm,obs,StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions)
function [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm, obs, StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ...
ParticleOptions, ThreadsOptions, DynareOptions, Model)
%
% Computes the proposal with a gaussian approximation for importance
% sampling
@ -19,6 +21,7 @@ function [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,St
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2017 Dynare Team
%
% This file is part of Dynare.
@ -36,83 +39,73 @@ function [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,St
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent init_flag2 mf0 mf1
persistent number_of_state_variables number_of_observed_variables
persistent number_of_structural_innovations
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ...
any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ...
any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4))
ghx
ghu
ghxx
ghuu
ghxu
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
else
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
end
constant = ReducedForm.constant;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set persistent variables.
if isempty(init_flag2)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
init_flag2 = 1;
end
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
if ParticleOptions.proposal_approximation.montecarlo
nodes = randn(ParticleOptions.number_of_particles,number_of_state_variables+number_of_structural_innovations) ;
nodes = randn(ParticleOptions.number_of_particles, number_of_state_variables+number_of_structural_innovations) ;
weights = 1/ParticleOptions.number_of_particles ;
weights_c = weights ;
elseif ParticleOptions.proposal_approximation.cubature
[nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations) ;
weights_c = weights ;
elseif ParticleOptions.proposal_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,ParticleOptions);
[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!')
error('This approximation for the proposal is not implemented or unknown!')
end
xbar = [StateVectorMean ; zeros(number_of_structural_innovations,1) ] ;
sqr_Px = [ [ StateVectorVarianceSquareRoot zeros(number_of_state_variables,number_of_structural_innovations) ] ;
[ zeros(number_of_structural_innovations,number_of_state_variables) Q_lower_triangular_cholesky ] ];
sigma_points = bsxfun(@plus,xbar,sqr_Px*(nodes'));
xbar = [StateVectorMean ; zeros(number_of_structural_innovations,1)] ;
sqr_Px = [ StateVectorVarianceSquareRoot, zeros(number_of_state_variables, number_of_structural_innovations);
zeros(number_of_structural_innovations, number_of_state_variables) Q_lower_triangular_cholesky];
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,ThreadsOptions.local_state_space_iteration_2);
PredictedStateMean = tmp(mf0,:)*weights ;
yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
else
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
end
PredictedStateMean = tmp(mf0,:)*weights;
PredictedObservedMean = tmp(mf1,:)*weights;
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);
dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean)'.*sqrt(weights);
PredictedStateMean = sum(PredictedStateMean, 2);
PredictedObservedMean = sum(PredictedObservedMean, 2);
dState = bsxfun(@minus,tmp(mf0,:), PredictedStateMean)'.*sqrt(weights);
dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean)'.*sqrt(weights);
PredictedStateVarianceSquareRoot = chol(dState'*dState)';
big_mat = [dObserved dState ; [H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)] ];
[mat1,mat] = qr2(big_mat,0);
big_mat = [dObserved, dState ; H_lower_triangular_cholesky, zeros(number_of_observed_variables,number_of_state_variables)];
[~, mat] = qr2(big_mat, 0);
mat = mat';
clear('mat1');
PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables);
CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables);
StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables));
PredictionError = obs - PredictedObservedMean;
StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*PredictionError;
else
dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean);
dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean);
dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean);
PredictedStateVariance = dState*diag(weights_c)*dState';
PredictedObservedVariance = dObserved*diag(weights_c)*dObserved' + H;
PredictedStateAndObservedCovariance = dState*diag(weights_c)*dObserved';
@ -123,4 +116,4 @@ else
StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';
end
end

View File

@ -1,6 +1,6 @@
function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions)
function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions, DynareOptions, Model)
% Copyright (C) 2013-2017 Dynare Team
% Copyright (C) 2013-2019 Dynare Team
%
% This file is part of Dynare.
%
@ -18,13 +18,22 @@ function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
mf1 = ReducedForm.mf1;
ghx = ReducedForm.ghx(mf1,:);
ghu = ReducedForm.ghu(mf1,:);
ghxx = ReducedForm.ghxx(mf1,:);
ghuu = ReducedForm.ghuu(mf1,:);
ghxu = ReducedForm.ghxu(mf1,:);
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
else
ghx = ReducedForm.ghx(mf1,:);
ghu = ReducedForm.ghu(mf1,:);
ghxx = ReducedForm.ghxx(mf1,:);
ghuu = ReducedForm.ghuu(mf1,:);
ghxu = ReducedForm.ghxu(mf1,:);
end
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,ThreadsOptions.local_state_space_iteration_2);
yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions);
measure = tmp(mf1,:);
else
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);
end