From 24980007195733ccc19975724ad5eda63801b032 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Sat, 21 Dec 2019 12:06:41 +0100 Subject: [PATCH] Allow k order approximation in Gaussian Filter (gf). Ref. dynare#1673 --- nonlinear-filters/src/gaussian_densities.m | 25 +++-- nonlinear-filters/src/gaussian_filter.m | 95 +++++++++---------- nonlinear-filters/src/gaussian_filter_bank.m | 93 +++++++++--------- nonlinear-filters/src/measurement_equations.m | 27 ++++-- 4 files changed, 115 insertions(+), 125 deletions(-) diff --git a/nonlinear-filters/src/gaussian_densities.m b/nonlinear-filters/src/gaussian_densities.m index f1c595100..64128cf58 100644 --- a/nonlinear-filters/src/gaussian_densities.m +++ b/nonlinear-filters/src/gaussian_densities.m @@ -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 . % 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; diff --git a/nonlinear-filters/src/gaussian_filter.m b/nonlinear-filters/src/gaussian_filter.m index e885cce72..29cc021e0 100644 --- a/nonlinear-filters/src/gaussian_filter.m +++ b/nonlinear-filters/src/gaussian_filter.m @@ -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 . -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). -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 \ No newline at end of file diff --git a/nonlinear-filters/src/measurement_equations.m b/nonlinear-filters/src/measurement_equations.m index 1c386b6fd..e450d1c11 100644 --- a/nonlinear-filters/src/measurement_equations.m +++ b/nonlinear-filters/src/measurement_equations.m @@ -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 . 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 \ No newline at end of file