diff --git a/nonlinear-filters/auxiliary_initialization.m b/nonlinear-filters/auxiliary_initialization.m new file mode 100644 index 000000000..e05741052 --- /dev/null +++ b/nonlinear-filters/auxiliary_initialization.m @@ -0,0 +1,117 @@ +function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,ParticleOptions,ThreadsOptions) + +% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling. + +% Copyright © 2011-2022 Dynare Team +% +% This file is part of Dynare (particles module). +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare particles module is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +persistent init_flag mf0 mf1 number_of_particles +persistent number_of_observed_variables number_of_structural_innovations + +% Set default +if isempty(start) + start = 1; +end + +% Set flag for prunning +%pruning = ParticleOptions.pruning; + +% Get steady state and mean. +%steadystate = ReducedForm.steadystate; +constant = ReducedForm.constant; +ss = ReducedForm.ys; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +% Set persistent variables. +if isempty(init_flag) + mf0 = ReducedForm.mf0; + mf1 = ReducedForm.mf1; + number_of_observed_variables = length(mf1); + number_of_structural_innovations = length(ReducedForm.Q); + number_of_particles = ParticleOptions.number_of_particles; + init_flag = 1; +end + +order = DynareOptions.order; + +% 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; +ghs2 = ReducedForm.ghs2; +if (order == 3) + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; +end + +% Get covariance matrices +Q = ReducedForm.Q; +H = ReducedForm.H; +if isempty(H) + H = 0; +end + +% Get initial condition for the state vector. +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)'; +state_variance_rank = size(StateVectorVarianceSquareRoot,2); +%Q_lower_triangular_cholesky = chol(Q)'; +%if pruning +% StateVectorMean_ = StateVectorMean; +% StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot; +%end + +% Set seed for randn(). +set_dynare_seed('default'); + +% Initialization of the likelihood. +const_lik = log(2*pi)*number_of_observed_variables; + +% Initialization of the weights across particles. +weights = ones(1,number_of_particles)/number_of_particles ; +StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); +%if pruning +% StateVectors_ = StateVectors; +%end +yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); +%if pruning +% yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state); +% [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); +%else +if (order == 2) + tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); +elseif (order == 3) + tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ss, options_.threads.local_state_space_iteration_3, false); +else + error('Orders > 3 not allowed'); +end +%end +PredictedObservedMean = weights*(tmp(mf1,:)'); +PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); +dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean'); +PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + H; +wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ; +tau_tilde = weights.*wtilde ; +tau_tilde = tau_tilde/sum(tau_tilde); +initial_distribution = resample(StateVectors',tau_tilde',ParticleOptions)' ; \ No newline at end of file diff --git a/nonlinear-filters/auxiliary_particle_filter.m b/nonlinear-filters/auxiliary_particle_filter.m new file mode 100644 index 000000000..21caa73bd --- /dev/null +++ b/nonlinear-filters/auxiliary_particle_filter.m @@ -0,0 +1,183 @@ +function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model) + +% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter +% allowing eventually resampling. +% +% Copyright © 2011-2022 Dynare Team +% +% This file is part of Dynare (particles module). +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare particles module is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Set default +if isempty(start) + start = 1; +end +% Get perturbation order +order = DynareOptions.order; + +% Set flag for prunning +pruning = ParticleOptions.pruning; + +% Get steady state and mean. +steadystate = ReducedForm.steadystate; +constant = ReducedForm.constant; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +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_structural_innovations = length(ReducedForm.Q); +number_of_particles = ParticleOptions.number_of_particles; + +if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; +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; + ghs2 = ReducedForm.ghs2; + if (order == 3) + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end +end + +% Get covariance matrices +Q = ReducedForm.Q; +H = ReducedForm.H; + +% Get initial condition for the state vector. +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)'; +state_variance_rank = size(StateVectorVarianceSquareRoot,2); +Q_lower_triangular_cholesky = chol(Q)'; + +% Set seed for randn(). +set_dynare_seed('default'); + +% Initialization of the likelihood. +const_lik = log(2*pi)*number_of_observed_variables+log(det(H)); +lik = NaN(sample_size,1); +LIK = NaN; + +% Initialization of the weights across particles. +weights = ones(1,number_of_particles)/number_of_particles ; +StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); +%StateVectors = bsxfun(@plus,zeros(state_variance_rank,number_of_particles),StateVectorMean); +if pruning + if order == 2 + StateVectors_ = StateVectors; + state_variables_steady_state_ = state_variables_steady_state; + mf0_ = mf0; + elseif order == 3 + StateVectors_ = repmat(StateVectors,3,1); + state_variables_steady_state_ = repmat(state_variables_steady_state,3,1); + mf0_ = repmat(mf0,1,3); + mask2 = number_of_state_variables+1:2*number_of_state_variables; + mask3 = 2*number_of_state_variables+1:3*number_of_state_variables; + mf0_(mask2) = mf0_(mask2)+size(ghx,1); + mf0_(mask3) = mf0_(mask3)+2*size(ghx,1); + else + error('Pruning is not available for orders > 3'); + end +end + +for t=1:sample_size + yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); + if pruning + yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state_); + if order == 2 + [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat_, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning); + else + error('Pruning is not available for orders > 3'); + end + else + if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations,number_of_particles), dr, Model, DynareOptions, udr); + else + if order == 2 + tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning); + else + error('Order > 3: use_k_order_solver should be set to true'); + end + end + end + PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); + z = sum(PredictionError.*(H\PredictionError),1) ; + tau_tilde = weights.*(tpdf(z,3*ones(size(z)))+1e-99) ; + tau_tilde = tau_tilde/sum(tau_tilde) ; + indx = resample(0,tau_tilde',ParticleOptions); + if pruning + yhat_ = yhat_(:,indx) ; + end + yhat = yhat(:,indx) ; + weights_stage_1 = weights(indx)./tau_tilde(indx) ; + epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles); + if pruning + if order == 2 + [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning); + else + error('Pruning is not available for orders > 3'); + end + StateVectors_ = tmp_(mf0_,:); + else + if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); + else + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning); + else + error('Order > 3: use_k_order_solver should be set to true'); + end + end + end + StateVectors = tmp(mf0,:); + PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); + weights_stage_2 = weights_stage_1.*(exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) + 1e-99) ; + lik(t) = log(mean(weights_stage_2)) ; + weights = weights_stage_2/sum(weights_stage_2); + if (ParticleOptions.resampling.status.generic && neff(weights). + +flag = false; + +order = DynareOptions.order; + +if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; +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; + ghs2 = ReducedForm.ghs2; + if order == 3 + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end +end + +constant = ReducedForm.constant; +steadystate = ReducedForm.steadystate; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +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/10, number_of_structural_innovations); + weights = 1.0/ParticleOptions.number_of_particles; + weights_c = weights; +elseif ParticleOptions.proposal_approximation.cubature + [nodes, weights] = spherical_radial_sigma_points(number_of_structural_innovations); + weights_c = weights; +elseif ParticleOptions.proposal_approximation.unscented + [nodes, weights, weights_c] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions); +else + error('Estimation: This approximation for the proposal is not implemented or unknown!') +end + +epsilon = Q_lower_triangular_cholesky*nodes'; +yhat = repmat(StateVectors-state_variables_steady_state, 1, size(epsilon, 2)); + +if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); +else + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false); + else + error('Order > 3: use_k_order_solver should be set to true'); + end +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); + PredictedStateVariance = dState*dState'; + big_mat = [dObserved dState; H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)]; + [~, mat] = qr2(big_mat,0); + mat = mat'; + 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)); + Error = y-PredictedObservedMean; + StateVectorMean = PredictedStateMean+(CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*Error; + if ParticleOptions.cpf_weights_method.amisanotristani + Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1), PredictedObservedVarianceSquareRoot, Error); + end +else + 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'; + KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance; + Error = y-PredictedObservedMean; + StateVectorMean = PredictedStateMean+KalmanFilterGain*Error; + StateVectorVariance = PredictedStateVariance-KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; + StateVectorVariance = 0.5*(StateVectorVariance+StateVectorVariance'); + [StateVectorVarianceSquareRoot, p] = chol(StateVectorVariance, 'lower') ; + if p + flag = true; + ProposalStateVector = zeros(number_of_state_variables, 1); + Weights = 0.0; + return + end + if ParticleOptions.cpf_weights_method.amisanotristani + Weights = SampleWeights.*probability2(zeros(number_of_observed_variables, 1), chol(PredictedObservedVariance)', Error); + end +end + +ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot, 2), 1)+StateVectorMean; +if ParticleOptions.cpf_weights_method.murrayjonesparslow + PredictedStateVariance = 0.5*(PredictedStateVariance+PredictedStateVariance'); + [PredictedStateVarianceSquareRoot, p] = chol(PredictedStateVariance, 'lower'); + if p + flag = true; + ProposalStateVector = zeros(number_of_state_variables,1); + Weights = 0.0; + return + end + Prior = probability2(PredictedStateMean, PredictedStateVarianceSquareRoot, ProposalStateVector); + Posterior = probability2(StateVectorMean, StateVectorVarianceSquareRoot, ProposalStateVector); + Likelihood = probability2(y, H_lower_triangular_cholesky, measurement_equations(ProposalStateVector, ReducedForm, ThreadsOptions, DynareOptions, Model)); + Weights = SampleWeights.*Likelihood.*(Prior./Posterior); +end diff --git a/nonlinear-filters/conditional_particle_filter.m b/nonlinear-filters/conditional_particle_filter.m new file mode 100644 index 000000000..db204dc85 --- /dev/null +++ b/nonlinear-filters/conditional_particle_filter.m @@ -0,0 +1,110 @@ +function [LIK,lik] = conditional_particle_filter(ReducedForm, Y, s, ParticleOptions, ThreadsOptions, DynareOptions, Model) + +% Evaluates the likelihood of a non-linear model with a particle filter +% +% INPUTS +% - ReducedForm [structure] Matlab's structure describing the reduced form model. +% - Y [double] p×T matrix of (detrended) data, where p is the number of observed variables. +% - s [integer] scalar, likelihood evaluation starts at s (has to be smaller than T, the sample length provided in Y). +% - ParticlesOptions [struct] +% - ThreadsOptions [struct] +% - DynareOptions [struct] +% - Model [struct] +% +% OUTPUTS +% - LIK [double] scalar, likelihood +% - lik [double] (T-s+1)×1 vector, density of observations in each period. +% +% REMARKS +% - The proposal is built using the Kalman updating step for each particle. +% - we need draws in the errors distributions +% Whether we use Monte-Carlo draws from a multivariate gaussian distribution +% as in Amisano & Tristani (JEDC 2010). +% Whether we use multidimensional Gaussian sparse grids approximations: +% - a univariate Kronrod-Paterson Gaussian quadrature combined by the Smolyak +% operator (ref: Winschel & Kratzig, 2010). +% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009a,2009b). +% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1997, van der +% Merwe & Wan 2003). +% +% Pros: +% - Allows using current observable information in the proposal +% - The use of sparse grids Gaussian approximation is much faster than the Monte-Carlo approach +% Cons: +% - The use of the Kalman updating step may biais the proposal distribution since +% it has been derived in a linear context and is implemented in a nonlinear +% context. That is why particle resampling is performed. + +% Copyright © 2009-2020 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Set default for third input argument. +if isempty(s) + s = 1; +end + +T = size(Y,2); +p = length(ReducedForm.mf1); +n = ParticleOptions.number_of_particles; + +% Get covariance matrices +Q = ReducedForm.Q; +H = ReducedForm.H; +if isempty(H) + H = 0; + H_lower_triangular_cholesky = 0; +else + H_lower_triangular_cholesky = chol(H)'; +end + +% Get initial condition for the state vector. +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)'; +state_variance_rank = size(StateVectorVarianceSquareRoot, 2); +Q_lower_triangular_cholesky = chol(Q)'; + +% Set seed for randn(). +set_dynare_seed('default'); + +% Initialization of the likelihood. +lik = NaN(T, 1); +ks = 0; +StateParticles = bsxfun(@plus, StateVectorVarianceSquareRoot*randn(state_variance_rank, n), StateVectorMean); +SampleWeights = ones(1, n)/n; + +for t=1:T + flags = false(n, 1); + for i=1:n + [StateParticles(:,i), SampleWeights(i), flags(i)] = ... + conditional_filter_proposal(ReducedForm, Y(:,t), StateParticles(:,i), SampleWeights(i), Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ParticleOptions, ThreadsOptions, DynareOptions, Model); + end + if any(flags) + LIK = -Inf; + lik(t) = -Inf; + return + end + SumSampleWeights = sum(SampleWeights); + lik(t) = log(SumSampleWeights); + SampleWeights = SampleWeights./SumSampleWeights; + if (ParticleOptions.resampling.status.generic && neff(SampleWeights). + +[dim,Ndata] = size(X); +M = size(StateMu,2) ; +if check % Ensure that covariances don't collapse + MIN_COVAR_SQRT = sqrt(eps); + init_covars = StateSqrtP; +end +eold = -Inf; +for n=1:niters + % Calculate posteriors based on old parameters + [prior,likelihood,marginal,posterior] = probability3(StateMu,StateSqrtP,StateWeights,X,X_weights); + e = sum(log(marginal)); + if (n > 1 && abs((e - eold)/eold) < crit) + return; + else + eold = e; + end + new_pr = (sum(posterior,2))'; + StateWeights = new_pr/Ndata; + StateMu = bsxfun(@rdivide,(posterior*X')',new_pr); + for j=1:M + diffs = bsxfun(@minus,X,StateMu(:,j)); + tpost = (1/sqrt(new_pr(j)))*sqrt(posterior(j,:)); + diffs = bsxfun(@times,diffs,tpost); + [foo,tcov] = qr2(diffs',0); + StateSqrtP(:,:,j) = tcov'; + if check + if min(abs(diag(StateSqrtP(:,:,j)))) < MIN_COVAR_SQRT + StateSqrtP(:,:,j) = init_covars(:,:,j); + end + end + end +end diff --git a/nonlinear-filters/gaussian_densities.m b/nonlinear-filters/gaussian_densities.m new file mode 100644 index 000000000..727d083c4 --- /dev/null +++ b/nonlinear-filters/gaussian_densities.m @@ -0,0 +1,50 @@ +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 +% +% INPUTS +% reduced_form_model [structure] Matlab's structure describing the reduced form model. +% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors. +% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors. +% reduced_form_model.state.dr [structure] output of resol.m. +% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables. +% start [integer] scalar, likelihood evaluation starts at 'start'. +% smolyak_accuracy [integer] scalar. +% +% OUTPUTS +% LIK [double] scalar, likelihood +% lik [double] vector, density of observations in each period. +% +% REFERENCES +% +% NOTES +% The vector "lik" is used to evaluate the jacobian of the likelihood. + +% Copyright © 2009-2019 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% proposal density +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, ThreadsOptions, DynareOptions, Model); +likelihood = probability2(obs, sqrt(H), yt_t_1_i); + +IncrementalWeights = likelihood.*prior./proposal; diff --git a/nonlinear-filters/gaussian_filter.m b/nonlinear-filters/gaussian_filter.m new file mode 100644 index 000000000..9b1b29e4c --- /dev/null +++ b/nonlinear-filters/gaussian_filter.m @@ -0,0 +1,135 @@ +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. +% Gaussian approximation is done by: +% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009). +% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1995) +% - Monte-Carlo draws from a multivariate gaussian distribution. +% First and second moments of prior and posterior state densities are computed +% from the resulting nodes/particles and allows to generate new distributions at the +% following observation. +% Pros: The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles +% filters since it treats a lesser number of particles. Furthermore, in all cases, there is no need +% of resampling. +% Cons: estimations may be biaised if the model is truly non-gaussian +% since predictive and filtered densities are unimodal. +% +% INPUTS +% Reduced_Form [structure] Matlab's structure describing the reduced form model. +% Y [double] matrix of original observed variables. +% start [double] structural parameters. +% ParticleOptions [structure] Matlab's structure describing options concerning particle filtering. +% ThreadsOptions [structure] Matlab's structure. +% +% OUTPUTS +% LIK [double] scalar, likelihood +% lik [double] vector, density of observations in each period. +% +% REFERENCES +% +% NOTES +% The vector "lik" is used to evaluate the jacobian of the likelihood. + +% Copyright © 2009-2019 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Set default +if isempty(start) + start = 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 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 + +if ParticleOptions.distribution_approximation.montecarlo + set_dynare_seed('default'); +end + +% Get covariance matrices +Q = ReducedForm.Q; +H = ReducedForm.H; +if isempty(H) + H = 0; + H_lower_triangular_cholesky = 0; +else + H_lower_triangular_cholesky = reduced_rank_cholesky(H)'; +end + +% Get initial condition for the state vector. +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)'; +state_variance_rank = size(StateVectorVarianceSquareRoot,2); +Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)'; + +% Initialization of the likelihood. +const_lik = (2*pi)^(number_of_observed_variables/2) ; +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, 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, ... + 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, ... + 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; + if not(ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented) + if (ParticleOptions.resampling.status.generic && neff(SampleWeights). + +order = DynareOptions.order; + +if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; +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; + ghs2 = ReducedForm.ghs2; + if order == 3 + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end +end + +constant = ReducedForm.constant; +steadystate = ReducedForm.steadystate; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +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) ; + 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); +else + 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')); +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); +if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); +else + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false); + else + error('Order > 3: use_k_order_solver should be set to true'); + end +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); + PredictedStateVarianceSquareRoot = chol(dState'*dState)'; + big_mat = [dObserved, dState ; H_lower_triangular_cholesky, zeros(number_of_observed_variables,number_of_state_variables)]; + [~, mat] = qr2(big_mat, 0); + mat = mat'; + 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); + PredictedStateVariance = dState*diag(weights_c)*dState'; + PredictedObservedVariance = dObserved*diag(weights_c)*dObserved' + H; + PredictedStateAndObservedCovariance = dState*diag(weights_c)*dObserved'; + PredictedStateVarianceSquareRoot = chol(PredictedStateVariance)'; + PredictionError = obs - PredictedObservedMean; + KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance; + StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError; + StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; + StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance'); + StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)'; +end diff --git a/nonlinear-filters/gaussian_mixture_densities.m b/nonlinear-filters/gaussian_mixture_densities.m new file mode 100644 index 000000000..0cf9dd79e --- /dev/null +++ b/nonlinear-filters/gaussian_mixture_densities.m @@ -0,0 +1,55 @@ +function IncrementalWeights = gaussian_mixture_densities(obs, StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ... + StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ... + ReducedForm, ThreadsOptions, DynareOptions, Model) + +% Elements to calculate the importance sampling ratio +% +% INPUTS +% reduced_form_model [structure] Matlab's structure describing the reduced form model. +% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors. +% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors. +% reduced_form_model.state.dr [structure] output of resol.m. +% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables. +% start [integer] scalar, likelihood evaluation starts at 'start'. +% smolyak_accuracy [integer] scalar. +% +% OUTPUTS +% LIK [double] scalar, likelihood +% lik [double] vector, density of observations in each period. +% +% REFERENCES +% +% NOTES +% The vector "lik" is used to evaluate the jacobian of the likelihood. + +% Copyright © 2009-2019 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Compute the density of particles under the prior distribution +[~, ~, prior] = probability(StateMuPrior, StateSqrtPPrior, StateWeightsPrior, StateParticles); +prior = prior'; + +% Compute the density of particles under the proposal distribution +[~, ~, 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, ThreadsOptions, DynareOptions, Model); + +% likelihood +likelihood = probability2(obs, sqrt(H), yt_t_1_i); +IncrementalWeights = likelihood.*prior./proposal; diff --git a/nonlinear-filters/gaussian_mixture_filter.m b/nonlinear-filters/gaussian_mixture_filter.m new file mode 100644 index 000000000..e06461b85 --- /dev/null +++ b/nonlinear-filters/gaussian_mixture_filter.m @@ -0,0 +1,226 @@ +function [LIK, lik] = gaussian_mixture_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model) + +% 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). +% Each gaussian distribution is obtained whether +% - with a radial-spherical cubature +% - with scaled unscented sigma-points +% A Sparse grid Kalman Filter is implemented on each component of the mixture, +% which confers it a weight about current information. +% Information on the current observables is then embodied in the proposal +% distribution in which we draw particles, which allows +% - reaching a greater precision relatively to a standard particle filter, +% - reducing the number of particles needed, +% - still being faster. +% +% +% INPUTS +% reduced_form_model [structure] Matlab's structure describing the reduced form model. +% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors. +% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors. +% reduced_form_model.state.dr [structure] output of resol.m. +% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables. +% start [integer] scalar, likelihood evaluation starts at 'start'. +% +% OUTPUTS +% LIK [double] scalar, likelihood +% lik [double] vector, density of observations in each period. +% +% REFERENCES +% +% Van der Meerwe & Wan, Gaussian Mixture Sigma-Point Particle Filters for Sequential +% Probabilistic Inference in Dynamic State-Space Models. +% Heiss & Winschel, 2010, Journal of Applied Economics. +% Winschel & Kratzig, 2010, Econometrica. +% +% NOTES +% The vector "lik" is used to evaluate the jacobian of the likelihood. +% Copyright © 2009-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Set default +if isempty(start) + start = 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_structural_innovations = length(ReducedForm.Q); +G = ParticleOptions.mixture_state_variables; % number of GM components in state +number_of_particles = ParticleOptions.number_of_particles; + +% compute gaussian quadrature nodes and weights on states and shocks +if ParticleOptions.distribution_approximation.cubature + [nodes, weights] = spherical_radial_sigma_points(number_of_state_variables); +elseif ParticleOptions.distribution_approximation.unscented + [nodes, weights] = unscented_sigma_points(number_of_state_variables, ParticleOptions); +else + if ~ParticleOptions.distribution_approximation.montecarlo + error('This approximation for the proposal is unknown!') + end +end + +if ParticleOptions.distribution_approximation.montecarlo + set_dynare_seed('default'); +end + +% Get covariance matrices +Q = ReducedForm.Q; +H = ReducedForm.H; +if isempty(H) + H = 0; + H_lower_triangular_cholesky = 0; +else + H_lower_triangular_cholesky = reduced_rank_cholesky(H)'; +end +Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)'; + +% Initialize mixtures +StateWeights = ones(1, G)/G; +StateMu = ReducedForm.StateVectorMean; +StateSqrtP = zeros(number_of_state_variables, number_of_state_variables, G); +temp = reduced_rank_cholesky(ReducedForm.StateVectorVariance)'; +StateMu = bsxfun(@plus, StateMu, bsxfun(@times,diag(temp), (-(G-1)/2:1:(G-1)/2))/10); +for g=1:G + StateSqrtP(:,:,g) = temp/sqrt(G) ; +end + +if ~ParticleOptions.mixture_structural_shocks + StructuralShocksMu = zeros(1, number_of_structural_innovations); + StructuralShocksWeights = 1; + I = 1; + StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu'; + StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I); + StructuralShocksSqrtP(:,:,1) = Q_lower_triangular_cholesky; +elseif ParticleOptions.mixture_structural_shocks==1 + if ParticleOptions.proposal_approximation.cubature + [StructuralShocksMu, StructuralShocksWeights] = spherical_radial_sigma_points(number_of_structural_innovations); + StructuralShocksWeights = ones(size(StructuralShocksMu, 1), 1)*StructuralShocksWeights; + elseif ParticleOptions.proposal_approximation.unscented + [StructuralShocksMu, StructuralShocksWeights] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions); + else + if ~ParticleOptions.distribution_approximation.montecarlo + error('This approximation for the proposal is unknown!') + end + end + I = size(StructuralShocksWeights, 1); + StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu'; + StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I); + for i=1:I + StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky; + end +else + if ParticleOptions.proposal_approximation.cubature + [StructuralShocksMu, StructuralShocksWeights] = spherical_radial_sigma_points(number_of_structural_innovations); + StructuralShocksWeights = ones(size(StructuralShocksMu, 1), 1)*StructuralShocksWeights ; + elseif ParticleOptions.proposal_approximation.unscented + [StructuralShocksMu, StructuralShocksWeights] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions); + else + if ~ParticleOptions.distribution_approximation.montecarlo + error('This approximation for the proposal is unknown!') + end + end + I = size(StructuralShocksWeights, 1); + StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu'; + StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I); + for i=1:I + StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky/sqrt(StructuralShocksWeights(i)); + end +end + +ObservationShocksWeights = 1; +J = 1 ; + +Gprime = G*I; +Gsecond = G*I*J; +SampleWeights = ones(Gsecond, 1)/Gsecond; + +StateWeightsPrior = zeros(1,Gprime); +StateMuPrior = zeros(number_of_state_variables,Gprime); +StateSqrtPPrior = zeros(number_of_state_variables, number_of_state_variables, Gprime); + +StateWeightsPost = zeros(1, Gsecond); +StateMuPost = zeros(number_of_state_variables, Gsecond); +StateSqrtPPost = zeros(number_of_state_variables, number_of_state_variables, Gsecond); + +const_lik = (2*pi)^(.5*number_of_observed_variables); + +lik = NaN(sample_size, 1); +LIK = NaN; +for t=1:sample_size + % Build the proposal joint quadratures of Gaussian on states, structural + % shocks and observation shocks based on each combination of mixtures + for i=1:I + for j=1:J + for g=1:G + gprime = g + (i-1)*G; + gsecond = gprime + (j-1)*Gprime; + [StateMuPrior(:,gprime), StateSqrtPPrior(:,:,gprime), StateWeightsPrior(1,gprime), ... + StateMuPost(:,gsecond), StateSqrtPPost(:,:,gsecond), StateWeightsPost(1,gsecond)] = ... + gaussian_mixture_filter_bank(ReducedForm,Y(:,t), StateMu(:,g), StateSqrtP(:,:,g), StateWeights(g),... + StructuralShocksMu(:,i), StructuralShocksSqrtP(:,:,i), StructuralShocksWeights(i),... + ObservationShocksWeights(j), H, H_lower_triangular_cholesky, const_lik, ... + ParticleOptions, ThreadsOptions, DynareOptions, Model); + end + end + end + + % Normalize weights + StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior, 2); + StateWeightsPost = StateWeightsPost/sum(StateWeightsPost, 2); + + 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, ... + ReducedForm, ThreadsOptions, DynareOptions, Model); + SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights); + end + SumSampleWeights = sum(SampleWeights); + lik(t) = log(SumSampleWeights); + SampleWeights = SampleWeights./SumSampleWeights; + [~, SortedRandomIndx] = sort(rand(1,Gsecond)); + SortedRandomIndx = SortedRandomIndx(1:G); + indx = resample(0,SampleWeights,ParticleOptions); + indx = indx(SortedRandomIndx); + StateMu = StateMuPost(:,indx); + StateSqrtP = StateSqrtPPost(:,:,indx); + StateWeights = ones(1,G)/G; + else + % Sample particle in the proposal distribution, ie the posterior state GM + StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles); + IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ... + StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ... + ReducedForm, ThreadsOptions, DynareOptions, Model); + SampleWeights = IncrementalWeights/number_of_particles; + SumSampleWeights = sum(SampleWeights,1); + SampleWeights = SampleWeights./SumSampleWeights; + lik(t) = log(SumSampleWeights); + if (ParticleOptions.resampling.status.generic && neff(SampleWeights). + +order = DynareOptions.order; + +if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; +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; + ghs2 = ReducedForm.ghs2; + if order == 3 + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end +end + +constant = ReducedForm.constant; +steadystate = ReducedForm.steadystate; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +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); + +numb = number_of_state_variables+number_of_structural_innovations; + +if ParticleOptions.proposal_approximation.cubature + [nodes3, weights3] = spherical_radial_sigma_points(numb); + weights_c3 = weights3; +elseif ParticleOptions.proposal_approximation.unscented + [nodes3, weights3, weights_c3] = unscented_sigma_points(numb, ParticleOptions); +else + error('This approximation for the proposal is unknown!') +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); +if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); +else + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false); + else + error('Order > 3: use_k_order_solver should be set to true'); + end +end +PredictedStateMean = tmp(mf0,:)*weights3; +PredictedObservedMean = tmp(mf1,:)*weights3; + +if ParticleOptions.proposal_approximation.cubature + PredictedStateMean = sum(PredictedStateMean, 2); + PredictedObservedMean = sum(PredictedObservedMean, 2); + dState = (bsxfun(@minus, tmp(mf0,:), PredictedStateMean)').*sqrt(weights3); + dObserved = (bsxfun(@minus, tmp(mf1,:), PredictedObservedMean)').*sqrt(weights3); + PredictedStateVariance = dState'*dState; + big_mat = [dObserved, dState ; H_lower_triangular_cholesky, zeros(number_of_observed_variables, number_of_state_variables)]; + [~, mat] = qr2(big_mat, 0); + mat = mat'; + 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)); + iPredictedObservedVarianceSquareRoot = inv(PredictedObservedVarianceSquareRoot); + iPredictedObservedVariance = iPredictedObservedVarianceSquareRoot'*iPredictedObservedVarianceSquareRoot; + sqrdet = 1/sqrt(det(iPredictedObservedVariance)); + PredictionError = obs - PredictedObservedMean; + StateVectorMean = PredictedStateMean + CovarianceObservedStateSquareRoot*iPredictedObservedVarianceSquareRoot*PredictionError; +else + dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean); + dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean); + PredictedStateVariance = dState*diag(weights_c3)*dState'; + PredictedObservedVariance = dObserved*diag(weights_c3)*dObserved' + H; + PredictedStateAndObservedCovariance = dState*diag(weights_c3)*dObserved'; + sqrdet = sqrt(det(PredictedObservedVariance)); + iPredictedObservedVariance = inv(PredictedObservedVariance); + PredictionError = obs - PredictedObservedMean; + KalmanFilterGain = PredictedStateAndObservedCovariance*iPredictedObservedVariance; + StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError; + StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; + StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance'); + StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)'; +end + +data_lik_GM_g = exp(-0.5*PredictionError'*iPredictedObservedVariance*PredictionError)/abs(normfactO*sqrdet) + 1e-99; +StateMuPrior = PredictedStateMean; +StateSqrtPPrior = reduced_rank_cholesky(PredictedStateVariance)'; +StateWeightsPrior = StateWeights*StructuralShocksWeights; +StateMuPost = StateVectorMean; +StateSqrtPPost = StateVectorVarianceSquareRoot; +StateWeightsPost = StateWeightsPrior*ObservationShocksWeights*data_lik_GM_g; diff --git a/nonlinear-filters/importance_sampling.m b/nonlinear-filters/importance_sampling.m new file mode 100644 index 000000000..d477f9f50 --- /dev/null +++ b/nonlinear-filters/importance_sampling.m @@ -0,0 +1,28 @@ +function State_Particles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost,numP) + +% Copyright © 2013-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +[Xdim,Gsecond] = size(StateMuPost) ; +u = rand(numP,1); +[Nc,comp] = histc(u, cumsum([0; StateWeightsPost])); +State_Particles = zeros(Xdim,numP); +for k=1:Gsecond + idx = bsxfun(@eq,comp,k*ones(size(comp))) ; + State_Particles(:,idx) = StateSqrtPPost(:,:,k)*randn(Xdim,Nc(k)); +end +State_Particles= State_Particles + StateMuPost(:,comp); diff --git a/nonlinear-filters/measurement_equations.m b/nonlinear-filters/measurement_equations.m new file mode 100644 index 000000000..8ed06349d --- /dev/null +++ b/nonlinear-filters/measurement_equations.m @@ -0,0 +1,57 @@ +function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions, DynareOptions, Model) + +% Copyright © 2013-2022 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +order = DynareOptions.order; +mf1 = ReducedForm.mf1; +if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; +else + ghx = ReducedForm.ghx(mf1,:); + ghu = ReducedForm.ghu(mf1,:); + ghxx = ReducedForm.ghxx(mf1,:); + ghuu = ReducedForm.ghuu(mf1,:); + ghxu = ReducedForm.ghxu(mf1,:); + ghs2 = ReducedForm.ghs2(mf1,:); + if order == 3 + ghxxx = ReducedForm.ghxxx(mf1,:); + ghuuu = ReducedForm.ghuuu(mf1,:); + ghxxu = ReducedForm.ghxxu(mf1,:); + ghxuu = ReducedForm.ghxuu(mf1,:); + ghxss = ReducedForm.ghxss(mf1,:); + ghuss = ReducedForm.ghuss(mf1,:); + end +end +steadystate = ReducedForm.steadystate(mf1,:); +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); +if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions, udr); + measure = tmp(mf1,:); +else + if order == 2 + 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); + elseif order == 3 + measure = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, size(yhat,2)), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false); + else + error('Order > 3: use_k_order_solver should be set to true'); + end +end diff --git a/nonlinear-filters/multivariate_smooth_resampling.m b/nonlinear-filters/multivariate_smooth_resampling.m new file mode 100644 index 000000000..a83d348c2 --- /dev/null +++ b/nonlinear-filters/multivariate_smooth_resampling.m @@ -0,0 +1,72 @@ +function new_particles = multivariate_smooth_resampling(particles,weights) +% Smooth Resampling of the particles. + +%@info: +%! @deftypefn {Function File} {@var{new_particles} =} multivariate_smooth_resampling (@var{weights}, @var{particles}, @var{number_of_new_particles}, @var{number_of_partitions}) +%! @anchor{particle/multivariate_smooth_resampling} +%! @sp 1 +%! Smooth Resampling of the particles (multivariate version). +%! @sp 2 +%! @strong{Inputs} +%! @sp 1 +%! @table @ @var +%! @item weights +%! n*1 vector of doubles, particles' weights. +%! @item particles +%! n*1 vector of doubles, particles. +%! @item number_of_new_particles +%! Integer scalar. +%! @item number_of_partitions +%! Integer scalar. +%! @end table +%! @sp 2 +%! @strong{Outputs} +%! @sp 1 +%! @table @ @var +%! @item indx +%! number_of_new_particles*1 vector of doubles, new particles. +%! @end table +%! @sp 2 +%! @strong{This function is called by:} +%! @sp 1 +%! @ref{particle/sequantial_importance_particle_filter} +%! @sp 2 +%! @strong{This function calls:} +%! @sp 1 +%! @ref{particle/univariate_smooth_resampling} +%! @sp 2 +%! @end deftypefn +%@eod: + +% Copyright © 2012-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr +% stephane DOT adjemian AT univ DASH lemans DOT fr + +number_of_particles = length(weights); +number_of_states = size(particles,2); +[P,D] = eig(particles'*(bsxfun(@times,1/number_of_particles,particles))) ; +D = diag(D) ; +vectors = bsxfun(@times,P,sqrt(D)') ; +orthogonalized_particles = bsxfun(@rdivide,particles*vectors,D') ; +new_particles = zeros(number_of_particles,number_of_states) ; +for j=1:number_of_states + tout = sortrows( [orthogonalized_particles(:,j) weights],1) ; + new_particles(:,j) = univariate_smooth_resampling(tout(:,2),tout(:,1),number_of_particles) ; +end +new_particles = new_particles*(vectors') ; diff --git a/nonlinear-filters/mykmeans.m b/nonlinear-filters/mykmeans.m new file mode 100644 index 000000000..e7c424b00 --- /dev/null +++ b/nonlinear-filters/mykmeans.m @@ -0,0 +1,53 @@ +function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod) + +% Copyright © 2013-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +[n,m] = size(x) ; +indold = zeros(1,m) ; +if cod==0 + d = transpose(sum(bsxfun(@power,bsxfun(@minus,x,mean(x)),2))); + d = sortrows( [transpose(1:m) d],2) ; + d = d((1+(0:1:g-1))*m/g,1) ; + c = x(:,d); +else + c = init ; +end +for iter=1:300 + dist = zeros(g,m) ; + for i=1:g + dist(i,:) = sum(bsxfun(@power,bsxfun(@minus,x,c(:,i)),2)); + end + [rien,ind] = min(dist) ; + if isequal(ind,indold) + break ; + end + indold = ind ; + for i=1:g + lin = bsxfun(@eq,ind,i.*ones(1,m)) ; + h = x(:,lin) ; + c(:,i) = mean(h,2) ; + end +end +SqrtVariance = zeros(n,n,g) ; +Weights = zeros(1,g) ; +for i=1:g + temp = x(:,bsxfun(@eq,ind,i*ones(1,m))) ; + u = bsxfun(@minus,temp,mean(temp,2)); %temp-mean(temp,1)' ; + SqrtVariance(:,:,i) = chol( (u*u')/size(temp,2) )' ; + Weights(i) = size(temp,2)/m ; +end diff --git a/nonlinear-filters/neff.m b/nonlinear-filters/neff.m new file mode 100644 index 000000000..bd428fcd1 --- /dev/null +++ b/nonlinear-filters/neff.m @@ -0,0 +1,21 @@ +function n = neff(w) +% Evaluates the criterion for resampling + +% Copyright © 2013-2014 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +n = 1/dot(w,w); \ No newline at end of file diff --git a/nonlinear-filters/nonlinear_kalman_filter.m b/nonlinear-filters/nonlinear_kalman_filter.m new file mode 100644 index 000000000..ccae0731d --- /dev/null +++ b/nonlinear-filters/nonlinear_kalman_filter.m @@ -0,0 +1,182 @@ +function [LIK,lik] = nonlinear_kalman_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 a Kalman filter. +% Gaussian distribution approximation is done by: +% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009). +% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1995) +% - Monte-Carlo draws from a multivariate gaussian distribution. +% First and second moments of prior and posterior state densities are computed +% from the resulting nodes/particles and allows to generate new distributions at the +% following observation. +% Pros: The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles +% filters since it treats a lesser number of particles. +% Cons: 1. Application a linear projection formulae in a nonlinear context. +% 2. Parameter estimations may be biaised if the model is truly non-gaussian since predictive and +% filtered densities are unimodal. +% +% INPUTS +% Reduced_Form [structure] Matlab's structure describing the reduced form model. +% Y [double] matrix of original observed variables. +% start [double] structural parameters. +% ParticleOptions [structure] Matlab's structure describing options concerning particle filtering. +% ThreadsOptions [structure] Matlab's structure. +% +% OUTPUTS +% LIK [double] scalar, likelihood +% lik [double] vector, density of observations in each period. +% +% REFERENCES +% +% NOTES +% The vector "lik" is used to evaluate the jacobian of the likelihood. + +% Copyright © 2009-2022 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Set default +if isempty(start) + start = 1; +end + +order = DynareOptions.order; + +if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; +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; + ghs2 = ReducedForm.ghs2; + if (order == 3) + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end +end + +constant = ReducedForm.constant; +steadystate = ReducedForm.steadystate; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +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_structural_innovations = length(ReducedForm.Q); + +% compute gaussian quadrature nodes and weights on states and shocks +if ParticleOptions.proposal_approximation.montecarlo + 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); +else + error('Estimation: This approximation for the proposal is not implemented or unknown!') +end + +if ParticleOptions.distribution_approximation.montecarlo + set_dynare_seed('default'); +end + +% Get covariance matrices +H = ReducedForm.H; +H_lower_triangular_cholesky = chol(H)' ; +Q_lower_triangular_cholesky = chol(ReducedForm.Q)'; + +% Get initial condition for the state vector. +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)'; + +% Initialization of the likelihood. +lik = NaN(sample_size,1); +LIK = NaN; + +for t=1:sample_size + 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); + if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); + else + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false); + end + 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); + big_mat = [dObserved dState ; [H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)] ]; + [~, mat] = qr2(big_mat,0); + mat = mat'; + 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 = Y(:,t) - PredictedObservedMean; + StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*PredictionError; + else + 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'; + PredictionError = Y(:,t) - PredictedObservedMean; + KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance; + StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError; + StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; + [StateVectorVarianceSquareRoot, p]= chol(StateVectorVariance,'lower'); + if p + LIK=-Inf; + lik(t)=-Inf; + return + end + [~, p]= chol(PredictedObservedVariance,'lower'); + if p + LIK=-Inf; + lik(t)=-Inf; + return + end + end + lik(t) = log( sum(probability2(Y(:,t),H_lower_triangular_cholesky,tmp(mf1,:)).*weights,1) ) ; +end + +LIK = -sum(lik(start:end)); diff --git a/nonlinear-filters/online_auxiliary_filter.m b/nonlinear-filters/online_auxiliary_filter.m new file mode 100644 index 000000000..bbfb21e3e --- /dev/null +++ b/nonlinear-filters/online_auxiliary_filter.m @@ -0,0 +1,455 @@ +function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, DynareResults) + +% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters. +% +% INPUTS +% - xparam1 [double] n×1 vector, Initial condition for the estimated parameters. +% - DynareDataset [dseries] Sample used for estimation. +% - dataset_info [struct] Description of the sample. +% - DynareOptions [struct] Option values (options_). +% - Model [struct] Description of the model (M_). +% - EstimatedParameters [struct] Description of the estimated parameters (estim_params_). +% - BayesInfo [struct] Prior definition (bayestopt_). +% - DynareResults [struct] Results (oo_). +% +% OUTPUTS +% - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters). +% - pmode [double] n×1 vector, mode of the particles at the end of the sample (for the parameters). +% - pmedian [double] n×1 vector, median of the particles at the end of the sample (for the parameters). +% - pstdev [double] n×1 vector, st. dev. of the particles at the end of the sample (for the parameters). +% - p025 [double] n×1 vector, 2.5 percent of the particles are below p025(i) for i=1,…,n. +% - p975 [double] n×1 vector, 97.5 percent of the particles are below p975(i) for i=1,…,n. +% - covariance [double] n×n matrix, covariance of the particles at the end of the sample. + +% Copyright © 2013-2022 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% Set seed for randn(). +set_dynare_seed('default'); +pruning = DynareOptions.particle.pruning; +second_resample = DynareOptions.particle.resampling.status.systematic; +variance_update = true; + +bounds = prior_bounds(BayesInfo, DynareOptions.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding + +% initialization of state particles +[~, Model, DynareOptions, DynareResults, ReducedForm] = solve_model_for_online_filter(true, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults); + +order = DynareOptions.order; +mf0 = ReducedForm.mf0; +mf1 = ReducedForm.mf1; +number_of_particles = DynareOptions.particle.number_of_particles; +number_of_parameters = size(xparam1,1); +Y = DynareDataset.data; +sample_size = size(Y,1); +number_of_observed_variables = length(mf1); +number_of_structural_innovations = length(ReducedForm.Q); +liu_west_delta = DynareOptions.particle.liu_west_delta; + +% Get initial conditions for the state particles +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)'; +state_variance_rank = size(StateVectorVarianceSquareRoot,2); +StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); +if pruning + if order == 2 + StateVectors_ = StateVectors; + elseif order == 3 + StateVectors_ = repmat(StateVectors,3,1); + else + error('Pruning is not available for orders > 3'); + end +end + +% parameters for the Liu & West filter +small_a = (3*liu_west_delta-1)/(2*liu_west_delta); +b_square = 1-small_a*small_a; + +% Initialization of parameter particles +xparam = zeros(number_of_parameters,number_of_particles); +Prior = dprior(BayesInfo, DynareOptions.prior_trunc); +for i=1:number_of_particles + info = 12042009; + while info + candidate = Prior.draw(); + [info, Model, DynareOptions, DynareResults] = solve_model_for_online_filter(false, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults); + if ~info + xparam(:,i) = candidate(:); + end + end +end + +% Initialization of the weights of particles. +weights = ones(1,number_of_particles)/number_of_particles; + +% Initialization of the likelihood. +const_lik = log(2*pi)*number_of_observed_variables; +mean_xparam = zeros(number_of_parameters,sample_size); +mode_xparam = zeros(number_of_parameters,sample_size); +median_xparam = zeros(number_of_parameters,sample_size); +std_xparam = zeros(number_of_parameters,sample_size); +lb95_xparam = zeros(number_of_parameters,sample_size); +ub95_xparam = zeros(number_of_parameters,sample_size); + +%% The Online filter +for t=1:sample_size + if t>1 + fprintf('\nSubsample with %s first observations.\n\n', int2str(t)) + else + fprintf('\nSubsample with only the first observation.\n\n') + end + % Moments of parameters particles distribution + m_bar = xparam*(weights'); + temp = bsxfun(@minus,xparam,m_bar); + sigma_bar = (bsxfun(@times,weights,temp))*(temp'); + if variance_update + chol_sigma_bar = chol(b_square*sigma_bar)'; + end + % Prediction (without shocks) + fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam); + tau_tilde = zeros(1,number_of_particles); + for i=1:number_of_particles + % model resolution + [info, Model, DynareOptions, DynareResults, ReducedForm] = ... + solve_model_for_online_filter(false, fore_xparam(:,i), DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults); + if ~info(1) + steadystate = ReducedForm.steadystate; + state_variables_steady_state = ReducedForm.state_variables_steady_state; + % Set local state space model (second-order approximation). + if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; + else + steadystate = ReducedForm.steadystate; + constant = ReducedForm.constant; + % 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; + ghs2 = ReducedForm.ghs2; + if (order == 3) + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end + if pruning + if order == 2 + state_variables_steady_state_ = state_variables_steady_state; + elseif order == 3 + state_variables_steady_state_ = repmat(state_variables_steady_state,3,1); + else + error('Pruning is not available for orders > 3'); + end + end + + end + % particle likelihood contribution + yhat = bsxfun(@minus, StateVectors(:,i), state_variables_steady_state); + if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions, udr); + else + if pruning + yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state_); + if order == 2 + [tmp, ~] = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat_, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning); + else + error('Pruning is not available for orders > 3'); + end + else + if order == 2 + tmp = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning); + else + error('Order > 3: use_k_order_solver should be set to true'); + end + end + end + PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:)); + % Replace Gaussian density with a Student density with 3 degrees of freedom for fat tails. + z = sum(PredictionError.*(ReducedForm.H\PredictionError), 1) ; + tau_tilde(i) = weights(i).*(tpdf(z, 3*ones(size(z)))+1e-99) ; + end + end + % particles selection + tau_tilde = tau_tilde/sum(tau_tilde); + indx = resample(0, tau_tilde', DynareOptions.particle); + StateVectors = StateVectors(:,indx); + xparam = fore_xparam(:,indx); + if pruning + StateVectors_ = StateVectors_(:,indx); + end + w_stage1 = weights(indx)./tau_tilde(indx); + % draw in the new distributions + wtilde = zeros(1, number_of_particles); + for i=1:number_of_particles + info = 12042009; + counter=0; + while info(1) && counter =bounds.lb) && all(candidate<=bounds.ub) + % model resolution for new parameters particles + [info, Model, DynareOptions, DynareResults, ReducedForm] = ... + solve_model_for_online_filter(false, candidate, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults) ; + if ~info(1) + xparam(:,i) = candidate ; + steadystate = ReducedForm.steadystate; + state_variables_steady_state = ReducedForm.state_variables_steady_state; + % Set local state space model (second order approximation). + if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; + else + constant = ReducedForm.constant; + % 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; + ghs2 = ReducedForm.ghs2; + if (order == 3) + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end + if pruning + if order == 2 + state_variables_steady_state_ = state_variables_steady_state; + mf0_ = mf0; + elseif order == 3 + state_variables_steady_state_ = repmat(state_variables_steady_state,3,1); + mf0_ = repmat(mf0,1,3); + mask2 = number_of_state_variables+1:2*number_of_state_variables; + mask3 = 2*number_of_state_variables+1:number_of_state_variables; + mf0_(mask2) = mf0_(mask2)+size(ghx,1); + mf0_(mask3) = mf0_(mask3)+2*size(ghx,1); + else + error('Pruning is not available for orders > 3'); + end + end + end + % Get covariance matrices and structural shocks + epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations, 1); + % compute particles likelihood contribution + yhat = bsxfun(@minus,StateVectors(:,i), state_variables_steady_state); + if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); + else + if pruning + yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state_); + if order == 2 + [tmp, tmp_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning); + else + error('Pruning is not available for orders > 3'); + end + StateVectors_(:,i) = tmp_(mf0_,:); + else + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning); + else + error('Order > 3: use_k_order_solver should be set to true'); + end + end + end + StateVectors(:,i) = tmp(mf0,:); + PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:)); + wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1))); + end + end + if counter==DynareOptions.particle.liu_west_max_resampling_tries + fprintf('\nLiu & West particle filter: I haven''t been able to solve the model in %u tries.\n',DynareOptions.particle.liu_west_max_resampling_tries) + fprintf('Liu & West particle filter: The last error message was: %s\n',get_error_message(info)) + fprintf('Liu & West particle filter: You can try to increase liu_west_max_resampling_tries, but most\n') + fprintf('Liu & West particle filter: likely there is an issue with the model.\n') + error('Liu & West particle filter: unable to solve the model.') + end + end + end + % normalization + weights = wtilde/sum(wtilde); + if variance_update && (neff(weights)=0.025 + lb95_xparam(i,t) = temp(j,1); + pass1 = true; + end + if ~pass2 && cumulated_weights(j)>=0.5 + median_xparam(i,t) = temp(j,1); + pass2 = true; + end + if ~pass3 && cumulated_weights(j)>=0.975 + ub95_xparam(i,t) = temp(j,1); + pass3 = true; + end + end + end + end + str = sprintf(' Lower Bound (95%%) \t Mean \t\t\t Upper Bound (95%%)'); + for l=1:size(xparam,1) + str = sprintf('%s\n %5.4f \t\t %7.5f \t\t %5.4f', str, lb95_xparam(l,t), mean_xparam(l,t), ub95_xparam(l,t)); + end + disp(str) + disp('') +end + +pmean = xparam(:,sample_size); +pmode = mode_xparam(:,sample_size); +pstdev = std_xparam(:,sample_size) ; +p025 = lb95_xparam(:,sample_size) ; +p975 = ub95_xparam(:,sample_size) ; +pmedian = median_xparam(:,sample_size) ; +covariance = mat_var_cov; + +%% Plot parameters trajectory +TeX = DynareOptions.TeX; + +nr = ceil(sqrt(number_of_parameters)) ; +nc = floor(sqrt(number_of_parameters)); +nbplt = 1 ; + +if TeX + fidTeX = fopen([Model.fname '_param_traj.tex'],'w'); + fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n'); + fprintf(fidTeX,['%% ' datestr(now,0) '\n']); + fprintf(fidTeX,' \n'); +end + +for plt = 1:nbplt + hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories'); + for k=1:length(pmean) + subplot(nr,nc,k) + [name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions); + % Draw the surface for an interval containing 95% of the particles. + area(1:sample_size, ub95_xparam(k,:), 'FaceColor', [.9 .9 .9], 'BaseValue', min(lb95_xparam(k,:))); + hold on + area(1:sample_size, lb95_xparam(k,:), 'FaceColor', [1 1 1], 'BaseValue', min(lb95_xparam(k,:))); + % Draw the mean of particles. + plot(1:sample_size, mean_xparam(k,:), '-k', 'linewidth', 2) + if TeX + title(texname,'interpreter','latex') + else + title(name,'interpreter','none') + end + hold off + axis tight + drawnow + end + dyn_saveas(hh, [Model.fname '_param_traj' int2str(plt)], DynareOptions.nodisplay, DynareOptions.graph_format); + if TeX + % TeX eps loader file + fprintf(fidTeX,'\\begin{figure}[H]\n'); + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParamTraj%s}\n',Model.fname,int2str(plt)); + fprintf(fidTeX,'\\caption{Parameters trajectories.}'); + fprintf(fidTeX,'\\label{Fig:ParametersPlots:%s}\n',int2str(plt)); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,' \n'); + end +end + +% Plot Parameter Densities +number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two. +bandwidth = 0; % Rule of thumb optimal bandwidth parameter. +kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation. +for plt = 1:nbplt + hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities'); + for k=1:length(pmean) + subplot(nr,nc,k) + [name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions); + optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function); + [density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ... + number_of_particles, optimal_bandwidth, kernel_function); + plot(density(:,1), density(:,2)); + hold on + if TeX + title(texname,'interpreter','latex') + else + title(name,'interpreter','none') + end + hold off + axis tight + drawnow + end + dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format); + if TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format))) + % TeX eps loader file + fprintf(fidTeX, '\\begin{figure}[H]\n'); + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/nc,1),M_.fname,int2str(plt)); + fprintf(fidTeX,'\\caption{Parameter densities based on the Liu/West particle filter.}'); + fprintf(fidTeX,'\\label{Fig:ParameterDensities:%s}\n',int2str(plt)); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,' \n'); + end +end diff --git a/nonlinear-filters/probability.m b/nonlinear-filters/probability.m new file mode 100644 index 000000000..a934b48d8 --- /dev/null +++ b/nonlinear-filters/probability.m @@ -0,0 +1,39 @@ +function [prior,likelihood,C,posterior] = probability(mu,sqrtP,prior,X) + +% Copyright © 2013-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +[dim,nov] = size(X); +M = size(mu,2) ; +if nargout>1 + likelihood = zeros(M,nov); + normfact = (2*pi)^(dim/2); + for k=1:M + XX = bsxfun(@minus,X,mu(:,k)); + S = sqrtP(:,:,k); + foo = S \ XX; + likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S)))); + end +end +likelihood = likelihood + 1e-99; +if nargout>2 + C = prior*likelihood + 1e-99; +end +if nargout>3 + posterior = bsxfun(@rdivide,bsxfun(@times,prior',likelihood),C) + 1e-99 ; + posterior = bsxfun(@rdivide,posterior,sum(posterior,1)); +end diff --git a/nonlinear-filters/probability2.m b/nonlinear-filters/probability2.m new file mode 100644 index 000000000..e37cfe3d1 --- /dev/null +++ b/nonlinear-filters/probability2.m @@ -0,0 +1,37 @@ +function [density] = probability2(mu,S,X) +% +% Multivariate gaussian density +% +% INPUTS +% n [integer] scalar, number of variables. +% +% OUTPUTS +% nodes [double] nodes of the cubature +% weigths [double] associated weigths +% +% REFERENCES +% +% +% NOTES +% +% Copyright © 2009-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +dim = size(X,1) ; +normfact = bsxfun(@power,(2*pi),(dim/2)) ; +foo = S\(bsxfun(@minus,X,mu)) ; +density = exp(-0.5*sum(foo.*foo)')./abs((normfact*prod(diag(S)))) + 1e-99 ; \ No newline at end of file diff --git a/nonlinear-filters/probability3.m b/nonlinear-filters/probability3.m new file mode 100644 index 000000000..18dc442b8 --- /dev/null +++ b/nonlinear-filters/probability3.m @@ -0,0 +1,39 @@ +function [prior,likelihood,C,posterior] = probability3(mu,sqrtP,prior,X,X_weights) + +% Copyright © 2013-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +[dim,nov] = size(X); +M = size(mu,2) ; +if nargout>1 + likelihood = zeros(M,nov); + normfact = (2*pi)^(dim/2); + for k=1:M + XX = bsxfun(@minus,X,mu(:,k)); + S = sqrtP(:,:,k); + foo = S \ XX; + likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S)))); + end +end +wlikelihood = bsxfun(@times,X_weights,likelihood) + 1e-99; +if nargout>2 + C = prior*wlikelihood + 1e-99; +end +if nargout>3 + posterior = bsxfun(@rdivide,bsxfun(@times,prior',wlikelihood),C) + 1e-99 ; + posterior = bsxfun(@rdivide,posterior,sum(posterior,1)); +end diff --git a/nonlinear-filters/resample.m b/nonlinear-filters/resample.m new file mode 100644 index 000000000..fc50173c2 --- /dev/null +++ b/nonlinear-filters/resample.m @@ -0,0 +1,79 @@ +function resampled_output = resample(particles,weights,ParticleOptions) +% Resamples particles. +% if particles = 0, returns the resampling index (except for smooth resampling) +% Otherwise, returns the resampled particles set. + +%@info: +%! @deftypefn {Function File} {@var{indx} =} resample (@var{weights}, @var{method}) +%! @anchor{particle/resample} +%! @sp 1 +%! Resamples particles. +%! @sp 2 +%! @strong{Inputs} +%! @sp 1 +%! @table @ @var +%! @item weights +%! n*1 vector of doubles, particles' weights. +%! @item method +%! string equal to 'residual' or 'traditional'. +%! @end table +%! @sp 2 +%! @strong{Outputs} +%! @sp 1 +%! @table @ @var +%! @item indx +%! n*1 vector of intergers, indices. +%! @end table +%! @sp 2 +%! @strong{This function is called by:} +%! @sp 1 +%! @ref{particle/sequantial_importance_particle_filter} +%! @sp 2 +%! @strong{This function calls:} +%! @sp 1 +%! @ref{residual_resampling}, @ref{traditional_resampling} +%! @sp 2 +%! @end deftypefn +%@eod: + +% Copyright © 2011-2014 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +defaultmethod = 1; % For residual based method set this variable equal to 0. + +if defaultmethod + if ParticleOptions.resampling.method.kitagawa + resampled_output = traditional_resampling(particles,weights,rand); + elseif ParticleOptions.resampling.method.stratified + resampled_output = traditional_resampling(particles,weights,rand(size(weights))); + elseif ParticleOptions.resampling.method.smooth + if particles==0 + error('Particle = 0 is incompatible with this resampling method!') + end + resampled_output = multivariate_smooth_resampling(particles,weights); + else + error('Unknown sampling method!') + end +else + if ParticleOptions.resampling.method.kitagawa + resampled_output = residual_resampling(particles,weights,rand); + elseif ParticleOptions.resampling.method.stratified + resampled_output = residual_resampling(particles,weights,rand(size(weights))); + else + error('Unknown sampling method!') + end +end diff --git a/nonlinear-filters/residual_resampling.m b/nonlinear-filters/residual_resampling.m new file mode 100644 index 000000000..ee0ccf9ee --- /dev/null +++ b/nonlinear-filters/residual_resampling.m @@ -0,0 +1,143 @@ +function return_resample = residual_resampling(particles,weights,noise) +% Resamples particles. + +%@info: +%! @deftypefn {Function File} {@var{indx} =} residual_resampling (@var{weights}) +%! @anchor{particle/residual_resampling} +%! @sp 1 +%! Resamples particles. +%! @sp 2 +%! @strong{Inputs} +%! @sp 1 +%! @table @ @var +%! @item weights +%! n*1 vector of doubles, particles' weights. +%! @end table +%! @sp 2 +%! @strong{Outputs} +%! @sp 1 +%! @table @ @var +%! @item indx +%! n*1 vector of intergers, indices. +%! @end table +%! @sp 2 +%! @strong{This function is called by:} +%! @sp 1 +%! @ref{particle/resample} +%! @sp 2 +%! @strong{This function calls:} +%! @sp 2 +%! @end deftypefn +%@eod: + +% Copyright © 2011-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr +% stephane DOT adjemian AT univ DASH lemans DOT fr + +% What is the number of particles? +number_of_particles = length(weights); + +switch length(noise) + case 1 + kitagawa_resampling = 1; + case number_of_particles + kitagawa_resampling = 0; + otherwise + error(['particle::resampling: Unknown method! The size of the second argument (' inputname(3) ') is wrong.']) +end + +% Set vectors of indices. +jndx = 1:number_of_particles; +indx = zeros(1,number_of_particles); + +% Multiply the weights by the number of particles. +WEIGHTS = number_of_particles*weights; + +% Compute the integer part of the normalized weights. +iWEIGHTS = fix(WEIGHTS); + +% Compute the number of resample +number_of_trials = number_of_particles-sum(iWEIGHTS); + +if number_of_trials + WEIGHTS = (WEIGHTS-iWEIGHTS)/number_of_trials; + EmpiricalCDF = cumsum(WEIGHTS); + if kitagawa_resampling + u = (transpose(1:number_of_trials)-1+noise(:))/number_of_trials; + else + u = fliplr(cumprod(noise(1:number_of_trials).^(1./(number_of_trials:-1:1)))); + end + j=1; + for i=1:number_of_trials + while (u(i)>EmpiricalCDF(j)) + j=j+1; + end + iWEIGHTS(j)=iWEIGHTS(j)+1; + if kitagawa_resampling==0 + j=1; + end + end +end + +k=1; +for i=1:number_of_particles + if (iWEIGHTS(i)>0) + for j=k:k+iWEIGHTS(i)-1 + indx(j) = jndx(i); + end + end + k = k + iWEIGHTS(i); +end + +if particles==0 + return_resample = indx ; +else + return_resample = particles(indx,:) ; +end +%@test:1 +%$ % Define the weights +%$ weights = randn(2000,1).^2; +%$ weights = weights/sum(weights); +%$ % Initialize t. +%$ t = ones(1,1); +%$ +%$ try +%$ indx1 = residual_resampling(weights); +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ T = all(t); +%@eof:1 + +%@test:2 +%$ % Define the weights +%$ weights = exp(randn(2000,1)); +%$ weights = weights/sum(weights); +%$ % Initialize t. +%$ t = ones(1,1); +%$ +%$ try +%$ indx1 = residual_resampling(weights); +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ T = all(t); +%@eof:2 \ No newline at end of file diff --git a/nonlinear-filters/sequential_importance_particle_filter.m b/nonlinear-filters/sequential_importance_particle_filter.m new file mode 100644 index 000000000..c69ded02f --- /dev/null +++ b/nonlinear-filters/sequential_importance_particle_filter.m @@ -0,0 +1,185 @@ +function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model) + +% Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling). + +% Copyright © 2011-2022 Dynare Team +% +% This file is part of Dynare (particles module). +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare particles module is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +persistent init_flag +persistent mf0 mf1 +persistent number_of_particles number_of_state_variables +persistent sample_size number_of_observed_variables number_of_structural_innovations + +% Set default value for start +if isempty(start) + start = 1; +end + +% Set flag for prunning +pruning = ParticleOptions.pruning; + +% Get steady state and mean. +steadystate = ReducedForm.steadystate; +constant = ReducedForm.constant; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +order = DynareOptions.order; + +% Set persistent variables (if needed). +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_structural_innovations = length(ReducedForm.Q); + number_of_particles = ParticleOptions.number_of_particles; + init_flag = 1; +end + +if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; + udr = ReducedForm.udr; +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; + ghs2 = ReducedForm.ghs2; + if order == 3 + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end +end + +% Get covariance matrices. +Q = ReducedForm.Q; % Covariance matrix of the structural innovations. +H = ReducedForm.H; % Covariance matrix of the measurement errors. +if isempty(H) + H = 0; +end + +% Initialization of the likelihood. +const_lik = log(2*pi)*number_of_observed_variables +log(det(H)) ; +lik = NaN(sample_size,1); + +% Get initial condition for the state vector. +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';%reduced_rank_cholesky(ReducedForm.StateVectorVariance)'; +if pruning + StateVectorMean_ = StateVectorMean; + StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot; +end + +% Get the rank of StateVectorVarianceSquareRoot +state_variance_rank = size(StateVectorVarianceSquareRoot,2); + +% Factorize the covariance matrix of the structural innovations +Q_lower_triangular_cholesky = chol(Q)'; + +% Set seed for randn(). +set_dynare_seed('default'); + +% Initialization of the weights across particles. +weights = ones(1,number_of_particles)/number_of_particles ; +StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); +if pruning + if order == 2 + StateVectors_ = StateVectors; + state_variables_steady_state_ = state_variables_steady_state; + mf0_ = mf0; + elseif order == 3 + StateVectors_ = repmat(StateVectors,3,1); + state_variables_steady_state_ = repmat(state_variables_steady_state,3,1); + mf0_ = repmat(mf0,1,3); + mask2 = number_of_state_variables+1:2*number_of_state_variables; + mask3 = 2*number_of_state_variables+1:3*number_of_state_variables; + mf0_(mask2) = mf0_(mask2)+size(ghx,1); + mf0_(mask3) = mf0_(mask3)+2*size(ghx,1); + else + error('Pruning is not available for orders > 3'); + end +end + +% Loop over observations +for t=1:sample_size + yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); + epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles); + if pruning + yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state_); + if order == 2 + [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning); + else + error('Pruning is not available for orders > 3'); + end + else + if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); + else + if order == 2 + tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning); + else + error('Order > 3: use_k_order_solver should be set to true'); + end + end + end + %PredictedObservedMean = tmp(mf1,:)*transpose(weights); + PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); + %dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean); + %PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H; + %PredictedObservedVariance = H; + if rcond(H) > 1e-16 + lnw = -.5*(const_lik+sum(PredictionError.*(H\PredictionError),1)); + else + LIK = NaN; + return + end + dfac = max(lnw); + wtilde = weights.*exp(lnw-dfac); + lik(t) = log(sum(wtilde))+dfac; + weights = wtilde/sum(wtilde); + if (ParticleOptions.resampling.status.generic && neff(weights). + +persistent init_flag restrict_variables_idx state_variables_idx mf0 mf1 number_of_state_variables + +info = 0; + +%---------------------------------------------------- +% 1. Get the structural parameters & define penalties +%---------------------------------------------------- + +% Test if some parameters are smaller than the lower bound of the prior domain. +if any(xparam1bounds.ub) + info = 42; + return +end + +% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H). +Q = Model.Sigma_e; +H = Model.H; +for i=1:EstimatedParameters.nvx + k =EstimatedParameters.var_exo(i,1); + Q(k,k) = xparam1(i)*xparam1(i); +end +offset = EstimatedParameters.nvx; +if EstimatedParameters.nvn + for i=1:EstimatedParameters.nvn + H(i,i) = xparam1(i+offset)*xparam1(i+offset); + end + offset = offset+EstimatedParameters.nvn; +else + H = zeros(size(DynareDataset.data, 2)); +end + +% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite. +if EstimatedParameters.ncx + for i=1:EstimatedParameters.ncx + k1 =EstimatedParameters.corrx(i,1); + k2 =EstimatedParameters.corrx(i,2); + Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2)); + Q(k2,k1) = Q(k1,k2); + end + % Try to compute the cholesky decomposition of Q (possible iff Q is positive definite) + [~, testQ] = chol(Q); + if testQ + % The variance-covariance matrix of the structural innovations is not definite positive. + info = 43; + return + end + offset = offset+EstimatedParameters.ncx; +end + +% Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite. +if EstimatedParameters.ncn + corrn_observable_correspondence = EstimatedParameters.corrn_observable_correspondence; + for i=1:EstimatedParameters.ncn + k1 = corrn_observable_correspondence(i,1); + k2 = corrn_observable_correspondence(i,2); + H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2)); + H(k2,k1) = H(k1,k2); + end + % Try to compute the cholesky decomposition of H (possible iff H is positive definite) + [~, testH] = chol(H); + if testH + % The variance-covariance matrix of the measurement errors is not definite positive. + info = 44; + return + end + offset = offset+EstimatedParameters.ncn; +end + +% Update estimated structural parameters in Mode.params. +if EstimatedParameters.np > 0 + Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end); +end + +% Update Model.Sigma_e and Model.H. +Model.Sigma_e = Q; +Model.H = H; + +%------------------------------------------------------------------------------ +% 2. call model setup & reduction program +%------------------------------------------------------------------------------ + +warning('off', 'MATLAB:nearlySingularMatrix') +[~, ~, ~, info, Model, DynareResults] = ... + dynare_resolve(Model, DynareOptions, DynareResults, 'restrict'); +warning('on', 'MATLAB:nearlySingularMatrix') + +if info(1)~=0 + if nargout==5 + ReducedForm = 0; + end + return +end + +% Get decision rules and transition equations. +dr = DynareResults.dr; + +% Set persistent variables (first call). +if isempty(init_flag) + mf0 = BayesInfo.mf0; + mf1 = BayesInfo.mf1; + restrict_variables_idx = dr.restrict_var_list; + state_variables_idx = restrict_variables_idx(mf0); + number_of_state_variables = length(mf0); + init_flag = true; +end + + +% Return reduced form model. +if nargout>4 + ReducedForm.ghx = dr.ghx(restrict_variables_idx,:); + ReducedForm.ghu = dr.ghu(restrict_variables_idx,:); + ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx)); + if DynareOptions.order==2 + ReducedForm.use_k_order_solver = false; + ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:); + ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:); + ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:); + ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx); + ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:); + elseif DynareOptions.order>=3 + ReducedForm.use_k_order_solver = true; + ReducedForm.dr = dr; + ReducedForm.udr = folded_to_unfolded_dr(dr, Model, DynareOptions); + else + n_states=size(dr.ghx,2); + n_shocks=size(dr.ghu,2); + ReducedForm.use_k_order_solver = false; + ReducedForm.ghxx = zeros(size(restrict_variables_idx,1),n_states^2); + ReducedForm.ghuu = zeros(size(restrict_variables_idx,1),n_shocks^2); + ReducedForm.ghxu = zeros(size(restrict_variables_idx,1),n_states*n_shocks); + ReducedForm.constant = ReducedForm.steadystate; + end + ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx)); + ReducedForm.Q = Q; + ReducedForm.H = H; + ReducedForm.mf0 = mf0; + ReducedForm.mf1 = mf1; +end + +% Set initial condition +if setinitialcondition + switch DynareOptions.particle.initialization + case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model. + StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0); + [A,B] = kalman_transition_matrix(dr,dr.restrict_var_list,dr.restrict_columns,Model.exo_nbr); + StateVectorVariance2 = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold); + StateVectorVariance = lyapunov_symm(A, B*ReducedForm.Q*B', DynareOptions.lyapunov_fixed_point_tol, ... + DynareOptions.qz_criterium, DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); + StateVectorVariance = StateVectorVariance(mf0,mf0); + case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model). + StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0); + old_DynareOptionsperiods = DynareOptions.periods; + DynareOptions.periods = 5000; + old_DynareOptionspruning = DynareOptions.pruning; + DynareOptions.pruning = DynareOptions.particle.pruning; + y_ = simult(dr.ys, dr, Model, DynareOptions, DynareResults); + y_ = y_(dr.order_var(state_variables_idx),2001:DynareOptions.periods); + StateVectorVariance = cov(y_'); + DynareOptions.periods = old_DynareOptionsperiods; + DynareOptions.pruning = old_DynareOptionspruning; + clear('old_DynareOptionsperiods','y_'); + case 3% Initial state vector covariance is a diagonal matrix. + StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0); + StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables); + otherwise + error('Unknown initialization option!') + end + ReducedForm.StateVectorMean = StateVectorMean; + ReducedForm.StateVectorVariance = StateVectorVariance; +end diff --git a/nonlinear-filters/spherical_radial_sigma_points.m b/nonlinear-filters/spherical_radial_sigma_points.m new file mode 100644 index 000000000..756cac28f --- /dev/null +++ b/nonlinear-filters/spherical_radial_sigma_points.m @@ -0,0 +1,36 @@ +function [nodes,weights] = spherical_radial_sigma_points(n) +% +% Computes nodes and weigths from a third-degree spherical-radial cubature +% rule. +% INPUTS +% n [integer] scalar, number of variables. +% +% OUTPUTS +% nodes [double] nodes of the cubature +% weigths [double] associated weigths +% +% REFERENCES +% +% Arasaratnam & Haykin 2008,2009. +% +% NOTES +% +% Copyright © 2009-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +nodes = (sqrt(n)*([eye(n) -eye(n)]))' ; +weights = (1/(2*n)) ; diff --git a/nonlinear-filters/traditional_resampling.m b/nonlinear-filters/traditional_resampling.m new file mode 100644 index 000000000..ab9089c8c --- /dev/null +++ b/nonlinear-filters/traditional_resampling.m @@ -0,0 +1,237 @@ +function return_resample = traditional_resampling(particles,weights,noise) +% Resamples particles. + +%@info: +%! @deftypefn {Function File} {@var{indx} =} traditional_resampling (@var{weights},@var{noise}) +%! @anchor{particle/traditional_resampling} +%! @sp 1 +%! Resamples particles (Resampling à la Kitagawa or stratified resampling). +%! @sp 2 +%! @strong{Inputs} +%! @sp 1 +%! @table @ @var +%! @item weights +%! n*1 vector of doubles, particles' weights. +%! @item noise +%! n*1 vector of doubles sampled from a [0,1] uniform distribution (stratified resampling) or scalar double +%! sampled from a [0,1] uniform distribution (Kitagawa resampling). +%! @end table +%! @sp 2 +%! @strong{Outputs} +%! @sp 1 +%! @table @ @var +%! @item indx +%! n*1 vector of intergers, indices. +%! @end table +%! @sp 2 +%! @strong{This function is called by:} +%! @sp 1 +%! @ref{particle/resample} +%! @sp 2 +%! @strong{This function calls:} +%! @sp 2 +%! @end deftypefn +%@eod: + +% Copyright © 2011-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr +% stephane DOT adjemian AT univ DASH lemans DOT fr + +% What is the number of particles? +number_of_particles = length(weights); + +% Initialize the returned argument. +indx = ones(number_of_particles,1); + +% Select method. +switch length(noise) + case 1 + kitagawa_resampling = 1; + case number_of_particles + kitagawa_resampling = 0; + otherwise + error(['particle::resampling: Unknown method! The size of the second argument (' inputname(3) ') is wrong.']) +end + +% Get the empirical CDF. +c = cumsum(weights); + +% Draw a starting point. +if kitagawa_resampling + randvec = (transpose(1:number_of_particles)-1+noise(:))/number_of_particles ; +else + randvec = fliplr(cumprod(noise.^(1./(number_of_particles:-1:1)))); +end + +% Start at the bottom of the CDF +if kitagawa_resampling + j = 1; + for i=1:number_of_particles + while (randvec(i)>c(j)) + j = j+1; + end + indx(i) = j; + end +else + for i=1:number_of_particles + indx(i) = sum(randvec(i)>c); + end + % Matlab's indices start at 1... + indx = indx+1; +end + +if particles==0 + return_resample = indx ; +else + return_resample = particles(indx,:) ; +end + +%@test:1 +%$ % Define the weights +%$ weights = randn(2000,1).^2; +%$ weights = weights/sum(weights); +%$ % Initialize t. +%$ t = ones(2,1); +%$ +%$ % First, try the stratified resampling. +%$ try +%$ indx1 = traditional_resampling(weights,rand(2000,1)); +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ % Second, try the Kitagawa resampling. +%$ try +%$ indx2 = traditional_resampling(weights,rand); +%$ catch +%$ t(2) = 0; +%$ end +%$ +%$ T = all(t); +%@eof:1 + +%@test:2 +%$ % Define the weights +%$ weights = exp(randn(2000,1)); +%$ weights = weights/sum(weights); +%$ % Initialize t. +%$ t = ones(2,1); +%$ +%$ % First, try the stratified resampling. +%$ try +%$ indx1 = traditional_resampling(weights,rand(2000,1)); +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ % Second, try the Kitagawa resampling. +%$ try +%$ indx2 = traditional_resampling(weights,rand); +%$ catch +%$ t(2) = 0; +%$ end +%$ +%$ T = all(t); +%@eof:2 + +%@test:3 +%$ % Set the number of particles. +%$ number_of_particles = 20000; +%$ +%$ show_plot = 0; +%$ show_time = 1; +%$ +%$ % Define the weights +%$ weights = randn(number_of_particles,1).^2; +%$ weights = weights/sum(weights); +%$ +%$ % Compute the empirical CDF +%$ c = cumsum(weights); +%$ +%$ % Stratified resampling. +%$ noise = rand(number_of_particles,1); +%$ +%$ if show_time +%$ disp('Stratified resampling timing:') +%$ tic +%$ end +%$ +%$ indx1 = traditional_resampling(weights,noise); +%$ +%$ if show_time +%$ toc +%$ tic +%$ end +%$ +%$ indx1_ = zeros(number_of_particles,1); +%$ randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles; +%$ for i=1:number_of_particles +%$ j = 1; +%$ while (randvec(i)>c(j)) +%$ j = j + 1; +%$ end +%$ indx1_(i) = j; +%$ end +%$ +%$ if show_time +%$ toc +%$ end +%$ +%$ % Kitagawa's resampling. +%$ noise = rand; +%$ +%$ if show_time +%$ disp('Kitagawa''s resampling timing:') +%$ tic +%$ end +%$ +%$ indx2 = traditional_resampling(weights,noise); +%$ +%$ if show_time +%$ toc +%$ tic +%$ end +%$ +%$ indx2_ = zeros(number_of_particles,1); +%$ randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles; +%$ j = 1; +%$ for i=1:number_of_particles +%$ while (randvec(i)>c(j)) +%$ j = j + 1; +%$ end +%$ indx2_(i) = j; +%$ end +%$ +%$ if show_time +%$ toc +%$ end +%$ +%$ % REMARK +%$ % Note that the alternative code used in this test is sensibly faster than the code proposed +%$ % in the routine for the resampling scheme à la Kitagawa... +%$ +%$ if show_plot +%$ plot(randvec,c,'-r'), hold on, plot([randvec(1),randvec(end)],[c(1),c(end)],'-k'), hold off, axis tight, box on +%$ end +%$ +%$ % Check results. +%$ t(1) = dassert(indx1,indx1_); +%$ t(2) = dassert(indx2,indx2_); +%$ T = all(t); +%@eof:3 \ No newline at end of file diff --git a/nonlinear-filters/univariate_smooth_resampling.m b/nonlinear-filters/univariate_smooth_resampling.m new file mode 100644 index 000000000..6661d8f0d --- /dev/null +++ b/nonlinear-filters/univariate_smooth_resampling.m @@ -0,0 +1,82 @@ +function new_particles = univariate_smooth_resampling(weights,particles,number_of_new_particles) +% Smooth Resampling of the particles. + +%@info: +%! @deftypefn {Function File} {@var{new_particles} =} univariate_smooth_resampling (@var{weights}, @var{number_of_new_particles}) +%! @anchor{particle/univariate_smooth_resampling} +%! @sp 1 +%! Smooth Resampling of the particles (univariate version). +%! @sp 2 +%! @strong{Inputs} +%! @sp 1 +%! @table @ @var +%! @item weights +%! n*1 vector of doubles, particles' weights. +%! @item particles +%! n*1 vector of doubles, particles. +%! @item number_of_new_particles +%! Integer scalar. +%! @end table +%! @sp 2 +%! @strong{Outputs} +%! @sp 1 +%! @table @ @var +%! @item indx +%! number_of_new_particles*1 vector of doubles, new particles. +%! @end table +%! @sp 2 +%! @strong{This function is called by:} +%! @sp 1 +%! @ref{particle/sequantial_importance_particle_filter} +%! @sp 2 +%! @strong{This function calls:} +%! @sp 2 +%! @end deftypefn +%@eod: + +% Copyright © 2012-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr +% stephane DOT adjemian AT univ DASH lemans DOT fr + +M = length(particles) ; +lambda_tilde = [ (.5*weights(1)) ; + (.5*(weights(1:M-1)+weights(2:M))) ; + (.5*weights(M)) ] ; +lambda_bar = cumsum(lambda_tilde) ; +u = rand(1,1) ; +new_particles = zeros(number_of_new_particles,1) ; +rj = 0 ; +i = 1 ; +j = 1 ; +while i<=number_of_new_particles + u_j = ( i-1 + u)/number_of_new_particles ; + while u_j>lambda_bar(j) + rj = j ; + j = j+1 ; + end + if rj==0 + new_particles(i) = particles(1) ; + elseif rj==M + new_particles(i) = particles(M) ; + else + u_star = (u_j - lambda_bar(rj))./lambda_tilde(rj+1) ; + new_particles(i) = (particles(rj+1) - particles(rj))*u_star + particles(rj) ; + end + i = i+1 ; +end \ No newline at end of file diff --git a/nonlinear-filters/unscented_sigma_points.m b/nonlinear-filters/unscented_sigma_points.m new file mode 100644 index 000000000..9ff7a27e8 --- /dev/null +++ b/nonlinear-filters/unscented_sigma_points.m @@ -0,0 +1,40 @@ +function [nodes,W_m,W_c] = unscented_sigma_points(n,ParticleOptions) +% +% Computes nodes and weigths for a scaled unscented transform cubature +% INPUTS +% n [integer] scalar, number of variables. +% +% OUTPUTS +% nodes [double] nodes of the cubature +% weigths [double] associated weigths +% +% REFERENCES +% +% +% +% NOTES +% +% Copyright © 2009-2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +lambda = (ParticleOptions.unscented.alpha^2)*(n+ParticleOptions.unscented.kappa) - n ; +nodes = [ zeros(n,1) ( sqrt(n+lambda).*([ eye(n) -eye(n)]) ) ]' ; +W_m = lambda/(n+lambda) ; +W_c = W_m + (1-ParticleOptions.unscented.alpha^2+ParticleOptions.unscented.beta) ; +temp = ones(2*n,1)/(2*(n+lambda)) ; +W_m = [W_m ; temp] ; +W_c = [W_c ; temp] ;