From b577de9b871c74e19d4509047f0d6f24dab06419 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Telemachus=29?= Date: Mon, 10 Nov 2014 19:00:16 +0100 Subject: [PATCH] Imported files from DynareTeam/dynare github repository. --- src/auxiliary_initialization.m | 116 ++++++ src/auxiliary_particle_filter.m | 159 ++++++++ src/conditional_filter_proposal.m | 123 ++++++ src/conditional_particle_filter.m | 125 ++++++ src/fit_gaussian_mixture.m | 52 +++ src/gaussian_densities.m | 52 +++ src/gaussian_filter.m | 162 ++++++++ src/gaussian_filter_bank.m | 122 ++++++ src/gaussian_mixture_densities.m | 58 +++ src/gaussian_mixture_filter.m | 225 +++++++++++ src/gaussian_mixture_filter_bank.m | 138 +++++++ src/importance_sampling.m | 29 ++ src/index_resample.m | 72 ++++ .../local_state_space_iteration_2.m | 275 +++++++++++++ src/measurement_equations.m | 33 ++ src/multivariate_smooth_resampling.m | 72 ++++ src/mykmeans.m | 53 +++ src/neff.m | 21 + src/online_auxiliary_filter.m | 373 ++++++++++++++++++ src/probability.m | 39 ++ src/probability2.m | 38 ++ src/resample.m | 74 ++++ src/residual_resampling.m | 143 +++++++ src/sequential_importance_particle_filter.m | 178 +++++++++ src/solve_model_for_online_filter.m | 344 ++++++++++++++++ src/spherical_radial_sigma_points.m | 36 ++ src/traditional_resampling.m | 237 +++++++++++ src/univariate_smooth_resampling.m | 82 ++++ src/unscented_sigma_points.m | 42 ++ 29 files changed, 3473 insertions(+) create mode 100644 src/auxiliary_initialization.m create mode 100644 src/auxiliary_particle_filter.m create mode 100644 src/conditional_filter_proposal.m create mode 100644 src/conditional_particle_filter.m create mode 100644 src/fit_gaussian_mixture.m create mode 100644 src/gaussian_densities.m create mode 100644 src/gaussian_filter.m create mode 100644 src/gaussian_filter_bank.m create mode 100644 src/gaussian_mixture_densities.m create mode 100644 src/gaussian_mixture_filter.m create mode 100644 src/gaussian_mixture_filter_bank.m create mode 100644 src/importance_sampling.m create mode 100644 src/index_resample.m create mode 100644 src/local_state_space_iteration/local_state_space_iteration_2.m create mode 100644 src/measurement_equations.m create mode 100644 src/multivariate_smooth_resampling.m create mode 100644 src/mykmeans.m create mode 100644 src/neff.m create mode 100644 src/online_auxiliary_filter.m create mode 100644 src/probability.m create mode 100644 src/probability2.m create mode 100644 src/resample.m create mode 100644 src/residual_resampling.m create mode 100644 src/sequential_importance_particle_filter.m create mode 100644 src/solve_model_for_online_filter.m create mode 100644 src/spherical_radial_sigma_points.m create mode 100644 src/traditional_resampling.m create mode 100644 src/univariate_smooth_resampling.m create mode 100644 src/unscented_sigma_points.m diff --git a/src/auxiliary_initialization.m b/src/auxiliary_initialization.m new file mode 100644 index 000000000..784550000 --- /dev/null +++ b/src/auxiliary_initialization.m @@ -0,0 +1,116 @@ +function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,DynareOptions) +% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling. +% +% INPUTS +% ReducedForm [structure] Matlab's structure describing the reduced form model. +% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors. +% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors. +% ReducedForm.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'. +% mf [integer] pp*1 vector of indices. +% number_of_particles [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 (C) 2013 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 . +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 = DynareOptions.particle.pruning; + +% Get steady state and mean. +%steadystate = ReducedForm.steadystate; +constant = ReducedForm.constant; +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 = DynareOptions.particle.number_of_particles; + init_flag = 1; +end + +% 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; + +% 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,DynareOptions.threads.local_state_space_iteration_2); +%else + tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); +%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',DynareOptions)' ; \ No newline at end of file diff --git a/src/auxiliary_particle_filter.m b/src/auxiliary_particle_filter.m new file mode 100644 index 000000000..e869484fe --- /dev/null +++ b/src/auxiliary_particle_filter.m @@ -0,0 +1,159 @@ +function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,DynareOptions) +% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling. +% +% INPUTS +% ReducedForm [structure] Matlab's structure describing the reduced form model. +% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors. +% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors. +% ReducedForm.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'. +% mf [integer] pp*1 vector of indices. +% number_of_particles [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 (C) 2011-2013 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 . +persistent init_flag mf0 mf1 number_of_particles +persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations + +% Set default +if isempty(start) + start = 1; +end + +% Set flag for prunning +pruning = DynareOptions.particle.pruning; + +% Get steady state and mean. +steadystate = ReducedForm.steadystate; +constant = ReducedForm.constant; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +% Set persistent variables. +if isempty(init_flag) + mf0 = ReducedForm.mf0; + mf1 = ReducedForm.mf1; + sample_size = size(Y,2); + number_of_state_variables = length(mf0); + number_of_observed_variables = length(mf1); + number_of_structural_innovations = length(ReducedForm.Q); + number_of_particles = DynareOptions.particle.number_of_particles; + init_flag = 1; +end + +% 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; + +% 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; +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); +if pruning + StateVectors_ = StateVectors; +end +for t=1:sample_size + 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,DynareOptions.threads.local_state_space_iteration_2); + else + tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); + 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 ; + sum_tau_tilde = sum(tau_tilde) ; + %var_wtilde = wtilde-sum_tau_tilde ; + %var_wtilde = var_wtilde'*var_wtilde/(number_of_particles-1) ; + lik(t) = log(sum_tau_tilde) ; %+ .5*var_wtilde/(number_of_particles*(sum_tau_tilde*sum_tau_tilde)) ; + tau_tilde = tau_tilde/sum_tau_tilde; + if pruning + temp = resample([yhat' yhat_'],tau_tilde',DynareOptions); + yhat = temp(:,1:number_of_state_variables)' ; + yhat_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)' ; + else + yhat = resample(yhat',tau_tilde',DynareOptions)' ; + end + if pruning + [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2); + else + tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); + 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))) ; + epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles); + if pruning + [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2); + StateVectors_ = tmp_(mf0,:); + else + tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); + end + StateVectors = tmp(mf0,:); + PredictedObservedMean = mean(tmp(mf1,:),2); + PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); + dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean); + PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles + H; + lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))); + wtilde = lnw./wtilde; + weights = wtilde/sum(wtilde); +end + +LIK = -sum(lik(start:end)); \ No newline at end of file diff --git a/src/conditional_filter_proposal.m b/src/conditional_filter_proposal.m new file mode 100644 index 000000000..aa82a449a --- /dev/null +++ b/src/conditional_filter_proposal.m @@ -0,0 +1,123 @@ +function [ProposalStateVector,Weights] = conditional_filter_proposal(ReducedForm,obs,StateVectors,SampleWeights,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions,normconst2) +% +% Computes the proposal for each past particle using Gaussian approximations +% for the state errors and the Kalman filter +% +% 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. +% +% 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 (C) 2012-2013 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 + +persistent init_flag2 mf0 mf1 +persistent number_of_state_variables number_of_observed_variables +persistent number_of_structural_innovations + +% Set local state space model (first-order approximation). +ghx = ReducedForm.ghx; +ghu = ReducedForm.ghu; +% Set local state space model (second-order approximation). +ghxx = ReducedForm.ghxx; +ghuu = ReducedForm.ghuu; +ghxu = ReducedForm.ghxu; + +if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ... + any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ... + any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4)) + ghx + ghu + ghxx + ghuu + ghxu +end + +constant = ReducedForm.constant; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +% Set persistent variables. +if isempty(init_flag2) + mf0 = ReducedForm.mf0; + mf1 = ReducedForm.mf1; + number_of_state_variables = length(mf0); + number_of_observed_variables = length(mf1); + number_of_structural_innovations = length(ReducedForm.Q); + init_flag2 = 1; +end + +if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.proposal_approximation.montecarlo + [nodes,weights] = spherical_radial_sigma_points(number_of_structural_innovations); + weights_c = weights ; +elseif DynareOptions.particle.proposal_approximation.unscented + [nodes,weights,weights_c] = unscented_sigma_points(number_of_structural_innovations,DynareOptions); +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)) ; + +tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); + +PredictedStateMean = tmp(mf0,:)*weights ; +PredictedObservedMean = tmp(mf1,:)*weights; + +if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.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)] ]; + [mat1,mat] = qr2(big_mat,0); + mat = mat'; + clear('mat1'); + PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables); + CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables); + StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables)); + StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*(obs - PredictedObservedMean); +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 ; + StateVectorMean = PredictedStateMean + KalmanFilterGain*(obs - PredictedObservedMean); + StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; + StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance'); + StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)'; +end + +ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ; +ypred = measurement_equations(ProposalStateVector,ReducedForm,DynareOptions) ; +foo = H_lower_triangular_cholesky \ (obs - ypred) ; +likelihood = exp(-0.5*(foo)'*foo)/normconst2 + 1e-99 ; +Weights = SampleWeights.*likelihood; diff --git a/src/conditional_particle_filter.m b/src/conditional_particle_filter.m new file mode 100644 index 000000000..3f365f85c --- /dev/null +++ b/src/conditional_particle_filter.m @@ -0,0 +1,125 @@ +function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,DynareOptions) +% +% Evaluates the likelihood of a non-linear model with a particle filter +% - 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. +% +% 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 (C) 2009-2010 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 + +persistent init_flag mf0 mf1 +persistent number_of_particles +persistent sample_size number_of_state_variables number_of_observed_variables + +% Set default +if isempty(start) + start = 1; +end + +% Set persistent variables. +if isempty(init_flag) + mf0 = ReducedForm.mf0; + mf1 = ReducedForm.mf1; + sample_size = size(Y,2); + number_of_state_variables = length(mf0); + number_of_observed_variables = length(mf1); + init_flag = 1; + number_of_particles = DynareOptions.particle.number_of_particles ; +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)'; + +% Set seed for randn(). +set_dynare_seed('default'); + +% Initialization of the likelihood. +normconst2 = log(2*pi)*number_of_observed_variables*prod(diag(H_lower_triangular_cholesky)) ; +lik = NaN(sample_size,1); +LIK = NaN; + +ks = 0 ; + +StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); +SampleWeights = ones(1,number_of_particles)/number_of_particles ; +for t=1:sample_size + for i=1:number_of_particles + [StateParticles(:,i),SampleWeights(i)] = ... + conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions,normconst2) ; + end + SumSampleWeights = sum(SampleWeights) ; + lik(t) = log(SumSampleWeights) ; + SampleWeights = SampleWeights./SumSampleWeights ; + if (DynareOptions.particle.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] = probability(StateMu,StateSqrtP,StateWeights,X); + 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/src/gaussian_densities.m b/src/gaussian_densities.m new file mode 100644 index 000000000..3a770118d --- /dev/null +++ b/src/gaussian_densities.m @@ -0,0 +1,52 @@ +function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,DynareOptions) +% +% 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 (C) 2009-2010 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,DynareOptions) ; +eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ; +yt_t_1 = sum(yt_t_1_i*weigths1,2) ; +tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ; +Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ; +sqr_det = sqrt(det(Pyy)) ; +foo = (eta_t_i/Pyy).*eta_t_i ; +likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ; +IncrementalWeights = likelihood.*prior./proposal ; diff --git a/src/gaussian_filter.m b/src/gaussian_filter.m new file mode 100644 index 000000000..da3840c74 --- /dev/null +++ b/src/gaussian_filter.m @@ -0,0 +1,162 @@ +function [LIK,lik] = gaussian_filter(ReducedForm,Y,start,DynareOptions) +% 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 Kronrod-Paterson gaussian quadrature with a limited number of nodes. +% Multidimensional quadrature is obtained by the Smolyak operator (ref: Winschel & Kratzig, 2010). +% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2008,2009). +% - a scaled unscented transform cubature (ref: ) +% - 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. +% => The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles +% filters since it treats a lesser number of particles and there is no need +% of resampling. +% However, estimations may reveal biaised if the model is truly non-gaussian +% since predictive and filtered densities are unimodal. +% +% 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 (C) 2009-2013 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 . + +persistent init_flag mf0 mf1 +persistent nodes2 weights2 weights_c2 number_of_particles +persistent sample_size number_of_state_variables number_of_observed_variables + +% Set default +if isempty(start) + start = 1; +end + +% Set persistent variables. +if isempty(init_flag) + mf0 = ReducedForm.mf0; + mf1 = ReducedForm.mf1; + sample_size = size(Y,2); + number_of_state_variables = length(mf0); + number_of_observed_variables = length(mf1); + number_of_particles = DynareOptions.particle.number_of_particles; + init_flag = 1; +end + +% compute gaussian quadrature nodes and weights on states and shocks +if isempty(nodes2) + if DynareOptions.particle.distribution_approximation.cubature + [nodes2,weights2] = spherical_radial_sigma_points(number_of_state_variables); + weights_c2 = weights2; + elseif DynareOptions.particle.distribution_approximation.unscented + [nodes2,weights2,weights_c2] = unscented_sigma_points(number_of_state_variables,DynareOptions); + else + if ~DynareOptions.particle.distribution_approximation.montecarlo + error('Estimation: This approximation for the proposal is not implemented or unknown!') + end + end +end + +if DynareOptions.particle.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; + +SampleWeights = 1/number_of_particles ; +ks = 0 ; +%Estimate = zeros(number_of_state_variables,sample_size) ; +%V_Estimate = zeros(number_of_state_variables,number_of_state_variables,sample_size) ; +for t=1:sample_size + % build the proposal + [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = ... + gaussian_filter_bank(ReducedForm,Y(:,t),StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions) ; + %Estimate(:,t) = PredictedStateMean ; + %V_Estimate(:,:,t) = PredictedStateVarianceSquareRoot ; + if DynareOptions.particle.distribution_approximation.cubature || DynareOptions.particle.distribution_approximation.unscented + StateParticles = bsxfun(@plus,StateVectorMean,StateVectorVarianceSquareRoot*nodes2') ; + IncrementalWeights = ... + gaussian_densities(Y(:,t),StateVectorMean,... + StateVectorVarianceSquareRoot,PredictedStateMean,... + PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,... + weights2,weights_c2,ReducedForm,DynareOptions) ; + SampleWeights = weights2.*IncrementalWeights ; + SumSampleWeights = sum(SampleWeights) ; + lik(t) = log(SumSampleWeights) ; + SampleWeights = SampleWeights./SumSampleWeights ; + else % Monte-Carlo draws + 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,DynareOptions) ; + SampleWeights = SampleWeights.*IncrementalWeights ; + SumSampleWeights = sum(SampleWeights) ; + %VarSampleWeights = IncrementalWeights-SumSampleWeights ; + %VarSampleWeights = VarSampleWeights*VarSampleWeights'/(number_of_particles-1) ; + lik(t) = log(SumSampleWeights) ; %+ .5*VarSampleWeights/(number_of_particles*(SumSampleWeights*SumSampleWeights)) ; + SampleWeights = SampleWeights./SumSampleWeights ; + Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ; + if (Neff<.5*sample_size && DynareOptions.particle.resampling.status.generic) || DynareOptions.particle.resampling.status.systematic + ks = ks + 1 ; + StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ; + StateVectorMean = mean(StateParticles,2) ; + StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )'; + SampleWeights = 1/number_of_particles ; + elseif DynareOptions.particle.resampling.status.none + StateVectorMean = (sampleWeights*StateParticles)' ; + temp = sqrt(SampleWeights').*StateParticles ; + StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp'*temp - StateVectorMean*(StateVectorMean') )'; + end + end +end + +LIK = -sum(lik(start:end)); diff --git a/src/gaussian_filter_bank.m b/src/gaussian_filter_bank.m new file mode 100644 index 000000000..6094dbcf1 --- /dev/null +++ b/src/gaussian_filter_bank.m @@ -0,0 +1,122 @@ +function [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = gaussian_filter_bank(ReducedForm,obs,StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions) +% +% Computes the proposal with a gaussian approximation for importance +% sampling +% This proposal is a gaussian distribution calculated à la Kalman +% +% 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. +% +% 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 (C) 2009-2012 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 . + +persistent init_flag2 mf0 mf1 +persistent number_of_state_variables number_of_observed_variables +persistent number_of_structural_innovations + +% Set local state space model (first-order approximation). +ghx = ReducedForm.ghx; +ghu = ReducedForm.ghu; +% Set local state space model (second-order approximation). +ghxx = ReducedForm.ghxx; +ghuu = ReducedForm.ghuu; +ghxu = ReducedForm.ghxu; + +if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ... + any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ... + any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4)) + ghx + ghu + ghxx + ghuu + ghxu +end + +constant = ReducedForm.constant; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +% Set persistent variables. +if isempty(init_flag2) + mf0 = ReducedForm.mf0; + mf1 = ReducedForm.mf1; + number_of_state_variables = length(mf0); + number_of_observed_variables = length(mf1); + number_of_structural_innovations = length(ReducedForm.Q); + init_flag2 = 1; +end + +if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.proposal_approximation.montecarlo + [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations) ; + weights_c = weights ; +elseif DynareOptions.particle.proposal_approximation.unscented + [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,DynareOptions); +else + error('Estimation: 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); +tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); +PredictedStateMean = tmp(mf0,:)*weights ; +PredictedObservedMean = tmp(mf1,:)*weights; + +if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.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)] ]; + [mat1,mat] = qr2(big_mat,0); + mat = mat'; + clear('mat1'); + PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables); + CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables); + StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables)); + PredictionError = obs - PredictedObservedMean; + StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*PredictionError; +else + dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean); + dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean); + 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 \ No newline at end of file diff --git a/src/gaussian_mixture_densities.m b/src/gaussian_mixture_densities.m new file mode 100644 index 000000000..9a5b28e05 --- /dev/null +++ b/src/gaussian_mixture_densities.m @@ -0,0 +1,58 @@ +function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... + StateMuPost,StateSqrtPPost,StateWeightsPost,... + StateParticles,H,normconst,weigths1,weigths2,ReducedForm,DynareOptions) + +% +% 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 (C) 2009-2012 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 +[ras,ras,prior] = probability(StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateParticles) ; +prior = prior' ; +% Compute the density of particles under the proposal distribution +[ras,ras,proposal] = probability(StateMuPost,StateSqrtPPost,StateWeightsPost,StateParticles) ; +proposal = proposal' ; +% Compute the density of the current observation conditionally to each particle +yt_t_1_i = measurement_equations(StateParticles,ReducedForm,DynareOptions) ; +eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ; +yt_t_1 = sum(yt_t_1_i*weigths1,2) ; +tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ; +Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ; +sqr_det = sqrt(det(Pyy)) ; +foo = (eta_t_i/Pyy).*eta_t_i ; +likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ; +IncrementalWeights = likelihood.*prior./proposal ; + diff --git a/src/gaussian_mixture_filter.m b/src/gaussian_mixture_filter.m new file mode 100644 index 000000000..5c60a81fb --- /dev/null +++ b/src/gaussian_mixture_filter.m @@ -0,0 +1,225 @@ +function [LIK,lik] = gaussian_mixture_filter(ReducedForm,Y,start,DynareOptions) +% 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 Smolyak quadrature à la Kronrod & Paterson (Heiss & Winschel 2010, Winschel & Kratzig 2010). +% - with a radial-spherical cubature +% - with scaled unscented sigma-points +% A Sparse grid Kalman Filter is implemented on each component of the mixture, +% 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 (C) 2009-2013 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 . + + +persistent init_flag mf0 mf1 Gprime Gsecond +persistent nodes weights weights_c I J G number_of_particles +persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations + +% Set default +if isempty(start) + start = 1; +end + +% Set persistent variables. +if isempty(init_flag) + mf0 = ReducedForm.mf0; + mf1 = ReducedForm.mf1; + sample_size = size(Y,2); + number_of_state_variables = length(mf0); + number_of_observed_variables = length(mf1); + number_of_structural_innovations = length(ReducedForm.Q); + G = DynareOptions.particle.mixture_state_variables; % number of GM components in state + I = DynareOptions.particle.mixture_structural_shocks ; % number of GM components in structural noise + J = DynareOptions.particle.mixture_measurement_shocks ; % number of GM components in observation noise + Gprime = G*I ; + Gsecond = G*I*J ; + number_of_particles = DynareOptions.particle.number_of_particles; + init_flag = 1; +end + +SampleWeights = ones(Gsecond,1)/Gsecond ; + +% compute gaussian quadrature nodes and weights on states and shocks +if isempty(nodes) + if DynareOptions.particle.distribution_approximation.cubature + [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables); + weights_c = weights; + elseif DynareOptions.particle.distribution_approximation.unscented + [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables,DynareOptions); + else + if ~DynareOptions.particle.distribution_approximation.montecarlo + error('Estimation: This approximation for the proposal is not implemented or unknown!') + end + end +end + +if DynareOptions.particle.distribution_approximation.montecarlo + set_dynare_seed('default'); + SampleWeights = 1/number_of_particles ; +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 all matrices +StateWeights = ones(1,G)/G ; +StateMu = ReducedForm.StateVectorMean*ones(1,G) ; +StateSqrtP = zeros(number_of_state_variables,number_of_state_variables,G) ; +for g=1:G + StateSqrtP(:,:,g) = reduced_rank_cholesky(ReducedForm.StateVectorVariance)' ; +end + +StructuralShocksWeights = ones(1,I)/I ; +StructuralShocksMu = zeros(number_of_structural_innovations,I) ; +StructuralShocksSqrtP = zeros(number_of_structural_innovations,number_of_structural_innovations,I) ; +for i=1:I + StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky ; +end + +ObservationShocksWeights = ones(1,J)/J ; +ObservationShocksMu = zeros(number_of_observed_variables,J) ; +ObservationShocksSqrtP = zeros(number_of_observed_variables,number_of_observed_variables,J) ; +for j=1:J + ObservationShocksSqrtP(:,:,j) = H_lower_triangular_cholesky ; +end + +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) ; + +%estimate = zeros(sample_size,number_of_state_variables,3) ; +const_lik = (2*pi)^(.5*number_of_observed_variables) ; + +ks = 0 ; +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 ; + a = g + (j-1)*G ; + b = a + (i-1)*Gprime ; + [StateMuPrior(:,a),StateSqrtPPrior(:,:,a),StateWeightsPrior(1,a),... + StateMuPost(:,b),StateSqrtPPost(:,:,b),StateWeightsPost(1,b)] =... + gaussian_mixture_filter_bank(ReducedForm,Y(:,t),StateMu(:,g),StateSqrtP(:,:,g),StateWeights(1,g),... + StructuralShocksMu(:,i),StructuralShocksSqrtP(:,:,i),StructuralShocksWeights(1,i),... + ObservationShocksMu(:,j),ObservationShocksSqrtP(:,:,j),ObservationShocksWeights(1,j),... + H,H_lower_triangular_cholesky,const_lik,DynareOptions) ; + end + end + end + + % Normalize weights + StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior,2) ; + StateWeightsPost = StateWeightsPost/sum(StateWeightsPost,2) ; + + if DynareOptions.particle.distribution_approximation.cubature || DynareOptions.particle.distribution_approximation.unscented + for i=1:Gsecond + StateParticles = bsxfun(@plus,StateMuPost(:,i),StateSqrtPPost(:,:,i)*nodes') ; + IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... + StateMuPost,StateSqrtPPost,StateWeightsPost,... + StateParticles,H,const_lik,weights,weights_c,ReducedForm,DynareOptions) ; + SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights) ; + end + SumSampleWeights = sum(SampleWeights) ; + lik(t) = log(SumSampleWeights) ; + SampleWeights = SampleWeights./SumSampleWeights ; + [ras,SortedRandomIndx] = sort(rand(1,Gsecond)); + SortedRandomIndx = SortedRandomIndx(1:G); + indx = index_resample(0,SampleWeights,DynareOptions) ; + 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) ; + % Compute prior, proposal and likelihood of particles + IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... + StateMuPost,StateSqrtPPost,StateWeightsPost,... + StateParticles,H,const_lik,1/number_of_particles,... + 1/number_of_particles,ReducedForm,DynareOptions) ; + % calculate importance weights of particles + SampleWeights = SampleWeights.*IncrementalWeights ; + SumSampleWeights = sum(SampleWeights,1) ; + SampleWeights = SampleWeights./SumSampleWeights ; + lik(t) = log(SumSampleWeights) ; + % First possible state point estimates + %estimate(t,:,1) = SampleWeights*StateParticles' ; + % Resampling if needed of required + Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ; + if (DynareOptions.particle.resampling.status.generic && Neff<.5*sample_size) || DynareOptions.particle.resampling.status.systematic + ks = ks + 1 ; + StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ; + StateVectorMean = mean(StateParticles,2) ; + StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/number_of_particles - StateVectorMean*(StateVectorMean') )'; + SampleWeights = 1/number_of_particles ; + elseif DynareOptions.particle.resampling.status.none + StateVectorMean = StateParticles*sampleWeights ; + temp = sqrt(SampleWeights').*StateParticles ; + StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )'; + end + % Use the information from particles to update the gaussian mixture on state variables + [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(StateParticles,StateMu,StateSqrtP,StateWeights,0.001,10,1) ; + %estimate(t,:,3) = StateWeights*StateMu' ; + end +end + +LIK = -sum(lik(start:end)) ; \ No newline at end of file diff --git a/src/gaussian_mixture_filter_bank.m b/src/gaussian_mixture_filter_bank.m new file mode 100644 index 000000000..3a3caa776 --- /dev/null +++ b/src/gaussian_mixture_filter_bank.m @@ -0,0 +1,138 @@ +function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPPost,StateWeightsPost] =... + gaussian_mixture_filter_bank(ReducedForm,obs,StateMu,StateSqrtP,StateWeights,... + StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,... + ObservationShocksMu,ObservationShocksSqrtP,ObservationShocksWeights,... + H,H_lower_triangular_cholesky,normfactO,DynareOptions) +% +% Computes the proposal with a gaussian approximation for importance +% sampling +% This proposal is a gaussian distribution calculated à la Kalman +% +% 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. +% +% 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 (C) 2009-2012 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 . + + + +persistent init_flag2 mf0 mf1 %nodes3 weights3 weights_c3 +persistent number_of_state_variables number_of_observed_variables +persistent number_of_structural_innovations + +% Set local state space model (first-order approximation). +ghx = ReducedForm.ghx; +ghu = ReducedForm.ghu; +% Set local state space model (second-order approximation). +ghxx = ReducedForm.ghxx; +ghuu = ReducedForm.ghuu; +ghxu = ReducedForm.ghxu; + +if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ... + any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ... + any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4)) + ghx + ghu + ghxx + ghuu + ghxu +end + +constant = ReducedForm.constant; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +% Set persistent variables. +if isempty(init_flag2) + mf0 = ReducedForm.mf0; + mf1 = ReducedForm.mf1; + number_of_state_variables = length(mf0); + number_of_observed_variables = length(mf1); + number_of_structural_innovations = length(ReducedForm.Q); + init_flag2 = 1; +end + +numb = number_of_state_variables+number_of_structural_innovations ; + +if DynareOptions.particle.proposal_approximation.cubature + [nodes3,weights3] = spherical_radial_sigma_points(numb); + weights_c3 = weights3; +elseif DynareOptions.particle.proposal_approximation.unscented + [nodes3,weights3,weights_c3] = unscented_sigma_points(numb,DynareOptions); +else + error('Estimation: This approximation for the proposal is not implemented or 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); +tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); +PredictedStateMean = tmp(mf0,:)*weights3; +PredictedObservedMean = tmp(mf1,:)*weights3; + +if DynareOptions.particle.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)] ]; + [mat1,mat] = qr2(big_mat,0); + mat = mat'; + clear('mat1'); + PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables); + CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables); + StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables)); + 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/src/importance_sampling.m b/src/importance_sampling.m new file mode 100644 index 000000000..403b6beb3 --- /dev/null +++ b/src/importance_sampling.m @@ -0,0 +1,29 @@ +function State_Particles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost,numP) + +% Copyright (C) 2013 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/src/index_resample.m b/src/index_resample.m new file mode 100644 index 000000000..e91645605 --- /dev/null +++ b/src/index_resample.m @@ -0,0 +1,72 @@ +function resampling_index = index_resample(particles,weights,DynareOptions) +% Resamples particles. + +%@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 (C) 2013 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 DynareOptions.particle.resampling.method.kitagawa + resampling_index = traditional_resampling(particles,weights,rand); + elseif DynareOptions.particle.resampling.method.stratified + resampling_index = traditional_resampling(particles,weights,rand(size(weights))); + else + error('Unknow sampling method!') + end +else + if DynareOptions.particle.resampling.method.kitagawa + resampled_particles = residual_resampling(particles,weights,rand); + elseif DynareOptions.particle.resampling.method.stratified + resampled_particles = residual_resampling(particles,weights,rand(size(weights))); + else + error('Unknown sampling method!') + end +end \ No newline at end of file diff --git a/src/local_state_space_iteration/local_state_space_iteration_2.m b/src/local_state_space_iteration/local_state_space_iteration_2.m new file mode 100644 index 000000000..0c6a52d07 --- /dev/null +++ b/src/local_state_space_iteration/local_state_space_iteration_2.m @@ -0,0 +1,275 @@ +function [y,y_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,a,b,c)%yhat_,ss) + +%@info: +%! @deftypefn {Function File} {@var{y}, @var{y_} =} local_state_equation_2 (@var{yhat},@var{epsilon}, @var{ghx}, @var{ghu}, @var{constant}, @var{ghxx}, @var{ghuu}, @var{ghxu}, @var{yhat_}, @var{ss}) +%! @anchor{particle/local_state_space_iteration_2} +%! @sp 1 +%! Given the states (y) and structural innovations (epsilon), this routine computes the level of selected endogenous variables when the +%! model is approximated by an order two taylor expansion around the deterministic steady state. Depending on the number of input/output +%! argument the pruning algorithm advocated by C. Sims is used or not (this version should not be used if the selected endogenous variables +%! are not the states of the model). +%! +%! @sp 2 +%! @strong{Inputs} +%! @sp 1 +%! @table @ @var +%! @item yhat +%! n*1 vector of doubles, initial condition, where n is the number of state variables. +%! @item epsilon +%! q*1 vector of doubles, structural innovations. +%! @item ghx +%! m*n matrix of doubles, restricted dr.ghx where we only consider the lines corresponding to a subset of endogenous variables. +%! @item ghu +%! m*q matrix of doubles, restricted dr.ghu where we only consider the lines corresponding to a subset of endogenous variables. +%! @item constant +%! m*1 vector of doubles, deterministic steady state plus second order correction for a subset of endogenous variables. +%! @item ghxx +%! m*n² matrix of doubles, restricted dr.ghxx where we only consider the lines corresponding to a subset of endogenous variables. +%! @item ghuu +%! m*q² matrix of doubles, restricted dr.ghuu where we only consider the lines corresponding to a subset of endogenous variables. +%! @item ghxu +%! m*(nq) matrix of doubles, subset of dr.ghxu where we only consider the lines corresponding to a subset of endogenous variables. +%! @item yhat_ +%! n*1 vector of doubles, spurious states for the pruning version. +%! @item ss +%! n*1 vector of doubles, steady state for the states. +%! @end table +%! @sp 2 +%! @strong{Outputs} +%! @sp 1 +%! @table @ @var +%! @item y +%! m*1 vector of doubles, selected endogenous variables. +%! @item y_ +%! m*1 vector of doubles, update of the latent variables needed for the pruning version (first order update). +%! @end table +%! @sp 2 +%! @strong{Remarks} +%! @sp 1 +%! [1] If the function has 10 input arguments then it must have 2 output arguments (pruning version). +%! @sp 1 +%! [2] If the function has 08 input arguments then it must have 1 output argument. +%! @sp 2 +%! @strong{This function is called by:} +%! @sp 2 +%! @strong{This function calls:} +%! +%! +%! @end deftypefn +%@eod: + +% Copyright (C) 2011-2012 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) stephane DOT adjemian AT univ DASH lemans DOT fr +% frederic DOT karame AT univ DASH evry DOT fr + +if nargin==9 + pruning = 0; numthreads = a; + if nargout>1 + error('local_state_space_iteration_2:: Numbers of input and output argument are inconsistent!') + end +elseif nargin==11 + pruning = 1; numthreads = c; yhat_ = a; ss = b; + if nargout~=2 + error('local_state_space_iteration_2:: Numbers of input and output argument are inconsistent!') + end +else + error('local_state_space_iteration_2:: Wrong number of input arguments!') +end + +number_of_threads = numthreads; + +switch pruning + case 0 + for i =1:size(yhat,2) + y(:,i) = constant + ghx*yhat(:,i) + ghu*epsilon(:,i) ... + + A_times_B_kronecker_C(.5*ghxx,yhat(:,i),number_of_threads) ... + + A_times_B_kronecker_C(.5*ghuu,epsilon(:,i),number_of_threads) ... + + A_times_B_kronecker_C(ghxu,yhat(:,i),epsilon(:,i),number_of_threads); + end + case 1 + for i =1:size(yhat,2) + y(:,i) = constant + ghx*yhat(:,i) + ghu*epsilon(:,i) ... + + A_times_B_kronecker_C(.5*ghxx,yhat_(:,i),number_of_threads) ... + + A_times_B_kronecker_C(.5*ghuu,epsilon(:,i),number_of_threads) ... + + A_times_B_kronecker_C(ghxu,yhat_(:,i),epsilon(:,i),number_of_threads); + end + y_ = ghx*yhat_+ghu*epsilon; + y_ = bsxfun(@plus,y_,ss); +end + +%@test:1 +%$ n = 2; +%$ q = 3; +%$ +%$ yhat = zeros(n,1); +%$ epsilon = zeros(q,1); +%$ ghx = rand(n,n); +%$ ghu = rand(n,q); +%$ constant = ones(n,1); +%$ ghxx = rand(n,n*n); +%$ ghuu = rand(n,q*q); +%$ ghxu = rand(n,n*q); +%$ yhat_ = zeros(n,1); +%$ ss = ones(n,1); +%$ +%$ % Call the tested routine. +%$ for i=1:10 +%$ y1 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1); +%$ [y2,y2_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,ss,1); +%$ end +%$ +%$ % Check the results. +%$ t(1) = dassert(y1,ones(n,1)); +%$ t(2) = dassert(y2,ones(n,1)); +%$ t(3) = dassert(y2_,ones(n,1)); +%$ T = all(t); +%@eof:1 + +%@test:2 +%$ old_path = pwd; +%$ cd([fileparts(which('dynare')) '/../tests/']); +%$ dynare('dsge_base2'); +%$ load dsge_base2; +%$ cd(old_path); +%$ dr = oo_.dr; +%$ clear('oo_','options_','M_'); +%$ delete([fileparts(which('dynare')) '/../tests/dsge_base2.mat']); +%$ istates = dr.nstatic+(1:dr.npred); +%$ n = dr.npred; +%$ q = size(dr.ghu,2); +%$ yhat = zeros(n,1); +%$ epsilon = zeros(q,1); +%$ ghx = dr.ghx(istates,:); +%$ ghu = dr.ghu(istates,:); +%$ constant = dr.ys(istates,:)+dr.ghs2(istates,:); +%$ ghxx = dr.ghxx(istates,:); +%$ ghuu = dr.ghuu(istates,:); +%$ ghxu = dr.ghxu(istates,:); +%$ yhat_ = zeros(n,1); +%$ ss = dr.ys(istates,:); +%$ +%$ t = ones(2,1); +%$ +%$ % Call the tested routine. +%$ try +%$ y1 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1); +%$ catch +%$ t(1) = 0; +%$ end +%$ try +%$ [y2,y2_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,ss,1); +%$ catch +%$ t(2) = 0; +%$ end +%$ +%$ % Check the results. +%$ T = all(t); +%@eof:2 + +%@test:3 +%$ Bohrbug = 1; % A bug that manifests reliably under a possibly unknown but well-defined set of conditions. +%$ if ~Bohrbug +%$ n = 2; +%$ q = 3; +%$ +%$ yhat = .01*randn(n,1); +%$ epsilon = .001*randn(q,1); +%$ ghx = rand(n,n); +%$ ghu = rand(n,q); +%$ constant = ones(n,1); +%$ ghxx = rand(n,n*n); +%$ ghuu = rand(n,q*q); +%$ ghxu = rand(n,n*q); +%$ yhat_ = zeros(n,1); +%$ ss = ones(n,1); +%$ +%$ % Call the tested routine (mex version). +%$ y1a = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1); +%$ [y2a,y2a_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,ss,1); +%$ +%$ % Call the tested routine (matlab version) +%$ path_to_mex = fileparts(which(['qmc_sequence.' mexext])); +%$ where_am_i_coming_from = pwd; +%$ cd(path_to_mex); +%$ tar('local_state_space_iteration_2.tar',['local_state_space_iteration_2.' mexext]); +%$ cd(where_am_i_coming_from); +%$ dynare_config([],0); +%$ y1b = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1); +%$ [y2b,y2b_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,ss,1); +%$ cd(path_to_mex); +%$ untar('local_state_space_iteration_2.tar'); +%$ delete('local_state_space_iteration_2.tar'); +%$ cd(where_am_i_coming_from); +%$ dynare_config([],0); +%$ % Check the results. +%$ t(1) = dassert(y1a,y1b); +%$ t(2) = dassert(y2a,y2b); +%$ t(3) = dassert(y2a_,y2b_); +%$ T = all(t); +%$ else +%$ t(1) = 1; +%$ T = all(t); +%$ end +%@eof:3 + + +%@test:4 +%$ % TIMING TEST (parallelization with openmp) +%$ old_path = pwd; +%$ cd([fileparts(which('dynare')) '/../tests/']); +%$ dynare('dsge_base2'); +%$ load dsge_base2; +%$ cd(old_path); +%$ dr = oo_.dr; +%$ clear('oo_','options_','M_'); +%$ delete([fileparts(which('dynare')) '/../tests/dsge_base2.mat']); +%$ istates = dr.nstatic+(1:dr.npred); +%$ n = dr.npred; +%$ q = size(dr.ghu,2); +%$ yhat = zeros(n,10000000); +%$ epsilon = zeros(q,10000000); +%$ ghx = dr.ghx(istates,:); +%$ ghu = dr.ghu(istates,:); +%$ constant = dr.ys(istates,:)+dr.ghs2(istates,:); +%$ ghxx = dr.ghxx(istates,:); +%$ ghuu = dr.ghuu(istates,:); +%$ ghxu = dr.ghxu(istates,:); +%$ yhat_ = zeros(n,10000000); +%$ ss = dr.ys(istates,:); +%$ +%$ t = NaN(4,1); +%$ tic, for i=1:10, y1 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1); end +%$ t1 = toc; +%$ tic, for i=1:10, y2 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,2); end +%$ t2 = toc; +%$ t(1) = dassert(y1,y2,1e-15); clear('y1'); +%$ tic, for i=1:10, y3 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,3); end +%$ t3 = toc; +%$ t(2) = dassert(y2,y3,1e-15); clear('y2'); +%$ tic, for i=1:10, y4 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,4); end +%$ t4 = toc; +%$ t(3) = dassert(y4,y3,1e-15); clear('y3','y4'); +%$ t(4) = (t1>t2) && (t2>t3) && (t3>t4); +%$ if ~t(4) +%$ disp('Timmings:') +%$ [t1, t2, t3, t4] +%$ end +%$ % Check the results. +%$ T = all(t); +%@eof:4 diff --git a/src/measurement_equations.m b/src/measurement_equations.m new file mode 100644 index 000000000..c00d8ebc4 --- /dev/null +++ b/src/measurement_equations.m @@ -0,0 +1,33 @@ +function measure = measurement_equations(StateVectors,ReducedForm,DynareOptions) + +% Copyright (C) 2013 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 . + +mf1 = ReducedForm.mf1; +ghx = ReducedForm.ghx(mf1,:); +ghu = ReducedForm.ghu(mf1,:); +ghxx = ReducedForm.ghxx(mf1,:); +ghuu = ReducedForm.ghuu(mf1,:); +ghxu = ReducedForm.ghxu(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) ; +measure = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,size(yhat,2)),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); + + + diff --git a/src/multivariate_smooth_resampling.m b/src/multivariate_smooth_resampling.m new file mode 100644 index 000000000..09ca287b9 --- /dev/null +++ b/src/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 (C) 2012-2013 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/src/mykmeans.m b/src/mykmeans.m new file mode 100644 index 000000000..56af5390b --- /dev/null +++ b/src/mykmeans.m @@ -0,0 +1,53 @@ +function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod) + +% Copyright (C) 2013 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/src/neff.m b/src/neff.m new file mode 100644 index 000000000..5251b9f17 --- /dev/null +++ b/src/neff.m @@ -0,0 +1,21 @@ +function n = neff(w) +% Evaluates the criterion for resampling + +% Copyright (C) 2013 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 = dot(w,w); diff --git a/src/online_auxiliary_filter.m b/src/online_auxiliary_filter.m new file mode 100644 index 000000000..4a8243a2d --- /dev/null +++ b/src/online_auxiliary_filter.m @@ -0,0 +1,373 @@ +function [xparam,std_param,lb_95,ub_95,median_param] = online_auxiliary_filter(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) + +% Carvalho & Lopes particle filter = auxiliary particle filter including Liu & West filter on parameters. +% +% INPUTS +% ReducedForm [structure] Matlab's structure describing the reduced form model. +% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors. +% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors. +% ReducedForm.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'. +% mf [integer] pp*1 vector of indices. +% number_of_particles [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 (C) 2013 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 . +persistent Y init_flag mf0 mf1 bounds number_of_particles number_of_parameters liu_west_delta liu_west_chol_sigma_bar +persistent start_param sample_size number_of_observed_variables number_of_structural_innovations + +% Set seed for randn(). +set_dynare_seed('default') ; +pruning = DynareOptions.particle.pruning; +second_resample = 1 ; +variance_update = 1 ; + +% initialization of state particles +[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... + solve_model_for_online_filter(1,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ; + +% Set persistent variables. +if isempty(init_flag) + 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,2); + number_of_observed_variables = length(mf1); + number_of_structural_innovations = length(ReducedForm.Q); + liu_west_delta = DynareOptions.particle.liu_west_delta ; + %liu_west_chol_sigma_bar = DynareOptions.particle.liu_west_chol_sigma_bar*eye(number_of_parameters) ; + start_param = xparam1 ; + %liu_west_chol_sigma_bar = sqrt(bsxfun(@times,eye(number_of_parameters),BayesInfo.p2)) ; + %start_param = BayesInfo.p1 ; + bounds = [BayesInfo.lb BayesInfo.ub] ; + init_flag = 1; +end + +% Get initial conditions for the state particles +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)'; +state_variance_rank = size(StateVectorVarianceSquareRoot,2); +StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); +if pruning + StateVectors_ = StateVectors; +end + +% parameters for the Liu & West filter +h_square = (3*liu_west_delta-1)/(2*liu_west_delta) ; +h_square = 1-h_square*h_square ; +small_a = sqrt(1-h_square) ; + +% Initialization of parameter particles +xparam = zeros(number_of_parameters,number_of_particles) ; +stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/100 ; +stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/50 ; +i = 1 ; +while i<=number_of_particles + %candidate = start_param + .001*liu_west_chol_sigma_bar*randn(number_of_parameters,1) ; + candidate = start_param + bsxfun(@times,stderr,randn(number_of_parameters,1)) ; + if all(candidate(:) >= bounds(:,1)) && all(candidate(:) <= bounds(:,2)) + xparam(:,i) = candidate(:) ; + i = i+1 ; + end +end + +%xparam = bsxfun(@plus,bounds(:,1),bsxfun(@times,(bounds(:,2)-bounds(:,1)),rand(number_of_parameters,number_of_particles))) ; + +% 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) ; +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 + disp(t) + % 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==1 + chol_sigma_bar = chol(h_square*sigma_bar)' ; + end + % Prediction (without shocks) + ObservedVariables = zeros(number_of_observed_variables,number_of_particles) ; + for i=1:number_of_particles + % model resolution + [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... + solve_model_for_online_filter(t,xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ; + steadystate = ReducedForm.steadystate; + state_variables_steady_state = ReducedForm.state_variables_steady_state; + % Set local state space model (second-order approximation). + constant = ReducedForm.constant; + ghx = ReducedForm.ghx; + ghu = ReducedForm.ghu; + ghxx = ReducedForm.ghxx; + ghuu = ReducedForm.ghuu; + ghxu = ReducedForm.ghxu; + % particle likelihood contribution + yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state); + if pruning + yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state); + [tmp, 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); + else + 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); + end + ObservedVariables(:,i) = tmp(mf1,:) ; + end + PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ; + PredictionError = bsxfun(@minus,Y(:,t),ObservedVariables); + dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean); + PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ; + wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ; + % unormalized weights and observation likelihood contribution + tau_tilde = weights.*wtilde ; + sum_tau_tilde = sum(tau_tilde) ; + % particles selection + tau_tilde = tau_tilde/sum_tau_tilde ; + indx = index_resample(0,tau_tilde',DynareOptions); + StateVectors = StateVectors(:,indx) ; + if pruning + StateVectors_ = StateVectors_(:,indx) ; + end + xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ; + xparam = xparam(:,indx) ; + wtilde = wtilde(indx) ; + % draw in the new distributions + i = 1 ; + while i<=number_of_particles + candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters,1) ; + if all(candidate >= bounds(:,1)) && all(candidate <= bounds(:,2)) + xparam(:,i) = candidate ; + % model resolution for new parameters particles + [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... + solve_model_for_online_filter(t,xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ; + steadystate = ReducedForm.steadystate; + state_variables_steady_state = ReducedForm.state_variables_steady_state; + % Set local state space model (second order approximation). + constant = ReducedForm.constant; + ghx = ReducedForm.ghx; + ghu = ReducedForm.ghu; + ghxx = ReducedForm.ghxx; + ghuu = ReducedForm.ghuu; + ghxu = ReducedForm.ghxu; + % 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 pruning + yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state); + [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2); + StateVectors_(:,i) = tmp_(mf0,:) ; + else + tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); + end + StateVectors(:,i) = tmp(mf0,:) ; + ObservedVariables(:,i) = tmp(mf1,:) ; + i = i+1 ; + end + end + PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ; + PredictionError = bsxfun(@minus,Y(:,t),ObservedVariables); + dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean); + PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ; + lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))); + % importance ratio + wtilde = lnw./wtilde ; + % normalization + weights = wtilde/sum(wtilde); + if (variance_update==1) && (neff(weights)0.025 && pass1==1 + lb95_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ; + pass1 = 2 ; + end + if cumulated_weights(j)>0.5 && pass2==1 + median_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ; + pass2 = 2 ; + end + if cumulated_weights(j)>0.975 && pass3==1 + ub95_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ; + pass3 = 2 ; + end + end + end + end + disp([lb95_xparam(:,t) mean_xparam(:,t) ub95_xparam(:,t)]) +end +distrib_param = xparam ; +xparam = mean_xparam(:,sample_size) ; +std_param = std_xparam(:,sample_size) ; +lb_95 = lb95_xparam(:,sample_size) ; +ub_95 = ub95_xparam(:,sample_size) ; +median_param = median_xparam(:,sample_size) ; + +%% Plot parameters trajectory +TeX = DynareOptions.TeX; + +[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters); + +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 + +z = 1:1:sample_size ; + +for plt = 1:nbplt, + if TeX + NAMES = []; + TeXNAMES = []; + end + hh = dyn_figure(DynareOptions,'Name','Parameters Trajectories'); + for k=1:min(nstar,length(xparam)-(plt-1)*nstar) + subplot(nr,nc,k) + kk = (plt-1)*nstar+k; + [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions); + if TeX + if isempty(NAMES) + NAMES = name; + TeXNAMES = texname; + else + NAMES = char(NAMES,name); + TeXNAMES = char(TeXNAMES,texname); + end + end + y = [mean_xparam(kk,:)' median_xparam(kk,:)' lb95_xparam(kk,:)' ub95_xparam(kk,:)' xparam(kk)*ones(sample_size,1)] ; + plot(z,y); + hold on + title(name,'interpreter','none') + hold off + axis tight + drawnow + end + dyn_saveas(hh,[ Model.fname '_param_traj' int2str(plt) ],DynareOptions); + if TeX + % TeX eps loader file + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for jj = 1:min(nstar,length(x)-(plt-1)*nstar) + fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:))); + end + 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, + if TeX + NAMES = []; + TeXNAMES = []; + end + hh = dyn_figure(DynareOptions,'Name','Parameters Densities'); + for k=1:min(nstar,length(xparam)-(plt-1)*nstar) + subplot(nr,nc,k) + kk = (plt-1)*nstar+k; + [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions); + if TeX + if isempty(NAMES) + NAMES = name; + TeXNAMES = texname; + else + NAMES = char(NAMES,name); + TeXNAMES = char(TeXNAMES,texname); + end + end + optimal_bandwidth = mh_optimal_bandwidth(distrib_param(kk,:)',number_of_particles,bandwidth,kernel_function); + [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(kk,:)',number_of_grid_points,... + number_of_particles,optimal_bandwidth,kernel_function); + plot(density(:,1),density(:,2)); + hold on + title(name,'interpreter','none') + hold off + axis tight + drawnow + end + dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions); + if TeX + % TeX eps loader file + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for jj = 1:min(nstar,length(x)-(plt-1)*nstar) + fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:))); + end + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParametersDensities%s}\n',Model.fname,int2str(plt)); + fprintf(fidTeX,'\\caption{ParametersDensities.}'); + fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt)); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,' \n'); + end +end + \ No newline at end of file diff --git a/src/probability.m b/src/probability.m new file mode 100644 index 000000000..5d76e49ee --- /dev/null +++ b/src/probability.m @@ -0,0 +1,39 @@ +function [prior,likelihood,C,posterior] = probability(mu,sqrtP,prior,X) + +% Copyright (C) 2013 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/src/probability2.m b/src/probability2.m new file mode 100644 index 000000000..57ab65f8f --- /dev/null +++ b/src/probability2.m @@ -0,0 +1,38 @@ +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 (C) 2009-2012 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 ; + diff --git a/src/resample.m b/src/resample.m new file mode 100644 index 000000000..a72894428 --- /dev/null +++ b/src/resample.m @@ -0,0 +1,74 @@ +function resampled_particles = resample(particles,weights,DynareOptions) +% Resamples particles. + +%@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 (C) 2011-2013 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 DynareOptions.particle.resampling.method.kitagawa + resampled_particles = traditional_resampling(particles,weights,rand); + elseif DynareOptions.particle.resampling.method.stratified + resampled_particles = traditional_resampling(particles,weights,rand(size(weights))); + elseif DynareOptions.particle.resampling.method.smooth + resampled_particles = multivariate_smooth_resampling(particles,weights); + else + error('Unknow sampling method!') + end +else + if DynareOptions.particle.resampling.method.kitagawa + resampled_particles = residual_resampling(particles,weights,rand); + elseif DynareOptions.particle.resampling.method.stratified + resampled_particles = residual_resampling(particles,weights,rand(size(weights))); + else + error('Unknown sampling method!') + end +end diff --git a/src/residual_resampling.m b/src/residual_resampling.m new file mode 100644 index 000000000..c2f09e3b8 --- /dev/null +++ b/src/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 (C) 2011-2013 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/src/sequential_importance_particle_filter.m b/src/sequential_importance_particle_filter.m new file mode 100644 index 000000000..eadaeba41 --- /dev/null +++ b/src/sequential_importance_particle_filter.m @@ -0,0 +1,178 @@ +function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,DynareOptions) +% Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling). +% Standard Sequential Monte Carlo approach with +% - the usual proposal (the state transition distribution) +% - options on resampling: none, adaptive or systematic +%@info: +%! @deftypefn {Function File} {@var{y}, @var{y_} =} sequential_importance_particle_filter (@var{ReducedForm},@var{Y}, @var{start}, @var{DynareOptions}) +%! @anchor{particle/sequential_importance_particle_filter} +%! @sp 1 +%! Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling). +%! +%! @sp 2 +%! @strong{Inputs} +%! @sp 1 +%! @table @ @var +%! @item ReducedForm +%! Structure describing the state space model (built in @ref{non_linear_dsge_likelihood}). +%! @item Y +%! p*smpl matrix of doubles (p is the number of observed variables), the (detrended) data. +%! @item start +%! Integer scalar, likelihood evaluation starts at observation 'start'. +%! @item DynareOptions +%! Structure specifying Dynare's options. +%! @end table +%! @sp 2 +%! @strong{Outputs} +%! @sp 1 +%! @table @ @var +%! @item LIK +%! double scalar, value of (minus) the logged likelihood. +%! @item lik +%! smpl*1 vector of doubles, density of the observations at each period. +%! @end table +%! @sp 2 +%! @strong{This function is called by:} +%! @ref{non_linear_dsge_likelihood} +%! @sp 2 +%! @strong{This function calls:} +%! +%! @end deftypefn +%@eod: + +% Copyright (C) 2011-2013 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 + +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 = DynareOptions.particle.pruning; + +% Get steady state and mean. +steadystate = ReducedForm.steadystate; +constant = ReducedForm.constant; +state_variables_steady_state = ReducedForm.state_variables_steady_state; + +% 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 = DynareOptions.particle.number_of_particles; + init_flag = 1; +end + +% 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; + +% 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; +lik = NaN(sample_size,1); + +% Get initial condition for the state vector. +StateVectorMean = ReducedForm.StateVectorMean; +StateVectorVarianceSquareRoot = 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 + StateVectors_ = StateVectors; +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); + [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2); + else + tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); + 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; + if rcond(PredictedObservedVariance) > 1e-16 + lnw = -.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\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 (DynareOptions.particle.resampling.status.generic && neff(weights). + +% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr +% frederic DOT karame AT univ DASH lemans DOT fr + +%global objective_function_penalty_base +% Declaration of the penalty as a persistent variable. +persistent init_flag +persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1 +persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations + +% Initialization of the returned arguments. +fval = []; +ys = []; +trend_coeff = []; +exit_flag = 1; + +% Set the number of observed variables +nvobs = DynareDataset.info.nvobs; + +%------------------------------------------------------------------------------ +% 1. Get the structural parameters & define penalties +%------------------------------------------------------------------------------ + +% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain. +%if (DynareOptions.mode_compute~=1) && any(xparam1BayesInfo.ub) +% k = find(xparam1(:)>BayesInfo.ub); +% fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2); +% exit_flag = 0; +% 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(nvobs); +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) +% [CholQ,testQ] = chol(Q); +% if testQ + % The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty. +% a = diag(eig(Q)); +% k = find(a < 0); +% if k > 0 +% fval = objective_function_penalty_base+sum(-a(k)); +% exit_flag = 0; +% info = 43; +% return +% end +% 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) +% [CholH,testH] = chol(H); +% if testH + % The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty. +% a = diag(eig(H)); +% k = find(a < 0); +% if k > 0 +% fval = objective_function_penalty_base+sum(-a(k)); +% exit_flag = 0; +% info = 44; +% return +% end +% 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 +%------------------------------------------------------------------------------ + +% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R). +[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict'); + +%if info(1) == 1 || info(1) == 2 || info(1) == 5 +% fval = objective_function_penalty_base+1; +% exit_flag = 0; +% return +%elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21 +% fval = objective_function_penalty_base+info(2); +% exit_flag = 0; +% return +%end + +% Define a vector of indices for the observed variables. Is this really usefull?... +BayesInfo.mf = BayesInfo.mf1; + +% Define the deterministic linear trend of the measurement equation. +if DynareOptions.noconstant + constant = zeros(nvobs,1); +else + if DynareOptions.loglinear + constant = log(SteadyState(BayesInfo.mfys)); + else + constant = SteadyState(BayesInfo.mfys); + end +end + +% Define the deterministic linear trend of the measurement equation. +if BayesInfo.with_trend + trend_coeff = zeros(DynareDataset.info.nvobs,1); + t = DynareOptions.trend_coeffs; + for i=1:length(t) + if ~isempty(t{i}) + trend_coeff(i) = evalin('base',t{i}); + end + end + trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs]; +else + trend = repmat(constant,1,DynareDataset.info.ntobs); +end + +% Get needed informations for kalman filter routines. +start = DynareOptions.presample+1; +np = size(T,1); +mf = BayesInfo.mf; +Y = transpose(DynareDataset.rawdata); + +%------------------------------------------------------------------------------ +% 3. Initial condition of the Kalman filter +%------------------------------------------------------------------------------ + +% 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 = BayesInfo.restrict_var_list; + observed_variables_idx = restrict_variables_idx(mf1); + state_variables_idx = restrict_variables_idx(mf0); + sample_size = size(Y,2); + number_of_state_variables = length(mf0); + number_of_observed_variables = length(mf1); + number_of_structural_innovations = length(Q); + init_flag = 1; +end + +ReducedForm.ghx = dr.ghx(restrict_variables_idx,:); +ReducedForm.ghu = dr.ghu(restrict_variables_idx,:); +ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:); +ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:); +ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:); +ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx)); +ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx); +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; + +% Set initial condition for t=1 +if observation_number==1 + 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.constant(mf0); + StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',1e-12,1e-12,[],[],DynareOptions.debug); + 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.constant(mf0); + old_DynareOptionsperiods = DynareOptions.periods; + DynareOptions.periods = 5000; + y_ = simult(oo_.steady_state, dr,Model,DynareOptions,DynareResults); + y_ = y_(state_variables_idx,2001:5000); + StateVectorVariance = cov(y_'); + DynareOptions.periods = old_DynareOptionsperiods; + clear('old_DynareOptionsperiods','y_'); + case 3% Initial state vector covariance is a diagonal matrix. + StateVectorMean = ReducedForm.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/src/spherical_radial_sigma_points.m b/src/spherical_radial_sigma_points.m new file mode 100644 index 000000000..754caacfa --- /dev/null +++ b/src/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 (C) 2009-2012 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/src/traditional_resampling.m b/src/traditional_resampling.m new file mode 100644 index 000000000..a9d86eb36 --- /dev/null +++ b/src/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 (C) 2011-2013 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/src/univariate_smooth_resampling.m b/src/univariate_smooth_resampling.m new file mode 100644 index 000000000..8238902a1 --- /dev/null +++ b/src/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 (C) 2012 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/src/unscented_sigma_points.m b/src/unscented_sigma_points.m new file mode 100644 index 000000000..26fcc4b5c --- /dev/null +++ b/src/unscented_sigma_points.m @@ -0,0 +1,42 @@ +function [nodes,W_m,W_c] = unscented_sigma_points(n,DynareOptions) +% +% 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 (C) 2009-2012 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 = (DynareOptions.particle.unscented.alpha^2)*(n+DynareOptions.particle.unscented.kappa) - n ; +nodes = [ zeros(n,1) ( sqrt(n+lambda).*([ eye(n) -eye(n)]) ) ]' ; +W_m = lambda/(n+lambda) ; +W_c = W_m + (1-DynareOptions.particle.unscented.alpha^2+DynareOptions.particle.unscented.beta) ; +temp = ones(2*n,1)/(2*(n+lambda)) ; +W_m = [W_m ; temp] ; +W_c = [W_c ; temp] ; + +