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

Merge particles submodule into master.

merge-particles-into-master
Stéphane Adjemian (Ryûk) 2023-07-12 09:33:20 +02:00
parent 6265a6d2b1 21dcc911bc
commit 8fc9109108
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
32 changed files with 3377 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitmodules vendored
View File

@ -10,9 +10,6 @@
[submodule "matlab/utilities/tests"]
path = matlab/utilities/tests
url = ../../Dynare/m-unit-tests.git
[submodule "matlab/particles"]
path = matlab/particles
url = ../../Dynare/particles.git
[submodule "matlab/modules/dseries"]
path = matlab/modules/dseries
url = ../../Dynare/dseries.git

2
matlab/dynare_config.m
View File

@ -53,7 +53,7 @@ p = {'/distributions/' ; ...
'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ...
'/../contrib/jsonlab/' ; ...
'/parallel/' ; ...
'/particles/src' ; ...
'/nonlinear-filters/' ; ...
'/gsa/' ; ...
'/ep/' ; ...
'/backward/' ; ...

117
matlab/nonlinear-filters/auxiliary_initialization.m Normal file
View File

@ -0,0 +1,117 @@
function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,ParticleOptions,ThreadsOptions)
% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling.
% Copyright © 2011-2022 Dynare Team
%
% This file is part of Dynare (particles module).
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare particles module is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
persistent init_flag mf0 mf1 number_of_particles
persistent number_of_observed_variables number_of_structural_innovations
% Set default
if isempty(start)
start = 1;
end
% Set flag for prunning
%pruning = ParticleOptions.pruning;
% Get steady state and mean.
%steadystate = ReducedForm.steadystate;
constant = ReducedForm.constant;
ss = ReducedForm.ys;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set persistent variables.
if isempty(init_flag)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
number_of_particles = ParticleOptions.number_of_particles;
init_flag = 1;
end
order = DynareOptions.order;
% Set local state space model (first order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if (order == 3)
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
if isempty(H)
H = 0;
end
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
%Q_lower_triangular_cholesky = chol(Q)';
%if pruning
% StateVectorMean_ = StateVectorMean;
% StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
%end
% Set seed for randn().
set_dynare_seed('default');
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables;
% Initialization of the weights across particles.
weights = ones(1,number_of_particles)/number_of_particles ;
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
%if pruning
% StateVectors_ = StateVectors;
%end
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
%if pruning
% yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state);
% [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
%else
if (order == 2)
tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
elseif (order == 3)
tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ss, options_.threads.local_state_space_iteration_3, false);
else
error('Orders > 3 not allowed');
end
%end
PredictedObservedMean = weights*(tmp(mf1,:)');
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean');
PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + H;
wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
tau_tilde = weights.*wtilde ;
tau_tilde = tau_tilde/sum(tau_tilde);
initial_distribution = resample(StateVectors',tau_tilde',ParticleOptions)' ;

183
matlab/nonlinear-filters/auxiliary_particle_filter.m Normal file
View File

@ -0,0 +1,183 @@
function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model)
% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter
% allowing eventually resampling.
%
% Copyright © 2011-2022 Dynare Team
%
% This file is part of Dynare (particles module).
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare particles module is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Set default
if isempty(start)
start = 1;
end
% Get perturbation order
order = DynareOptions.order;
% Set flag for prunning
pruning = ParticleOptions.pruning;
% Get steady state and mean.
steadystate = ReducedForm.steadystate;
constant = ReducedForm.constant;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
number_of_particles = ParticleOptions.number_of_particles;
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
% Set local state space model (first order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if (order == 3)
% Set local state space model (third order approximation).
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
end
% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn().
set_dynare_seed('default');
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables+log(det(H));
lik = NaN(sample_size,1);
LIK = NaN;
% Initialization of the weights across particles.
weights = ones(1,number_of_particles)/number_of_particles ;
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
%StateVectors = bsxfun(@plus,zeros(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
if order == 2
StateVectors_ = StateVectors;
state_variables_steady_state_ = state_variables_steady_state;
mf0_ = mf0;
elseif order == 3
StateVectors_ = repmat(StateVectors,3,1);
state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
mf0_ = repmat(mf0,1,3);
mask2 = number_of_state_variables+1:2*number_of_state_variables;
mask3 = 2*number_of_state_variables+1:3*number_of_state_variables;
mf0_(mask2) = mf0_(mask2)+size(ghx,1);
mf0_(mask3) = mf0_(mask3)+2*size(ghx,1);
else
error('Pruning is not available for orders > 3');
end
end
for t=1:sample_size
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
if pruning
yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state_);
if order == 2
[tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
[tmp, tmp_] = local_state_space_iteration_3(yhat_, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
else
error('Pruning is not available for orders > 3');
end
else
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations,number_of_particles), dr, Model, DynareOptions, udr);
else
if order == 2
tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end
end
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
z = sum(PredictionError.*(H\PredictionError),1) ;
tau_tilde = weights.*(tpdf(z,3*ones(size(z)))+1e-99) ;
tau_tilde = tau_tilde/sum(tau_tilde) ;
indx = resample(0,tau_tilde',ParticleOptions);
if pruning
yhat_ = yhat_(:,indx) ;
end
yhat = yhat(:,indx) ;
weights_stage_1 = weights(indx)./tau_tilde(indx) ;
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
if pruning
if order == 2
[tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
[tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
else
error('Pruning is not available for orders > 3');
end
StateVectors_ = tmp_(mf0_,:);
else
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
else
if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end
end
StateVectors = tmp(mf0,:);
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
weights_stage_2 = weights_stage_1.*(exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) + 1e-99) ;
lik(t) = log(mean(weights_stage_2)) ;
weights = weights_stage_2/sum(weights_stage_2);
if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
if pruning
temp = resample([StateVectors' StateVectors_'],weights',ParticleOptions);
StateVectors = temp(:,1:number_of_state_variables)';
StateVectors_ = temp(:,number_of_state_variables+1:end)';
else
StateVectors = resample(StateVectors',weights',ParticleOptions)';
end
weights = ones(1,number_of_particles)/number_of_particles;
end
end
LIK = -sum(lik(start:end));

164
matlab/nonlinear-filters/conditional_filter_proposal.m Normal file
View File

@ -0,0 +1,164 @@
function [ProposalStateVector, Weights, flag] = conditional_filter_proposal(ReducedForm, y, StateVectors, SampleWeights, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ...
H, ParticleOptions, ThreadsOptions, DynareOptions, Model)
% Computes the proposal for each past particle using Gaussian approximations
% for the state errors and the Kalman filter
%
% INPUTS
% - ReducedForm [structure] Matlab's structure describing the reduced form model.
% - y [double] p×1 vector, current observation (p is the number of observed variables).
% - StateVectors
% - SampleWeights
% - Q_lower_triangular_cholesky
% - H_lower_triangular_cholesky
% - H
% - ParticleOptions
% - ThreadsOptions
% - DynareOptions
% - Model
%
% OUTPUTS
% - ProposalStateVector
% - Weights
% - flag
% Copyright © 2012-2022 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
flag = false;
order = DynareOptions.order;
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if order == 3
% Set local state space model (third order approximation).
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
end
constant = ReducedForm.constant;
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
if ParticleOptions.proposal_approximation.montecarlo
nodes = randn(ParticleOptions.number_of_particles/10, number_of_structural_innovations);
weights = 1.0/ParticleOptions.number_of_particles;
weights_c = weights;
elseif ParticleOptions.proposal_approximation.cubature
[nodes, weights] = spherical_radial_sigma_points(number_of_structural_innovations);
weights_c = weights;
elseif ParticleOptions.proposal_approximation.unscented
[nodes, weights, weights_c] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions);
else
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
epsilon = Q_lower_triangular_cholesky*nodes';
yhat = repmat(StateVectors-state_variables_steady_state, 1, size(epsilon, 2));
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
else
if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end
PredictedStateMean = tmp(mf0,:)*weights;
PredictedObservedMean = tmp(mf1,:)*weights;
if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_approximation.montecarlo
PredictedStateMean = sum(PredictedStateMean, 2);
PredictedObservedMean = sum(PredictedObservedMean, 2);
dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean)'.*sqrt(weights);
dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean)'.*sqrt(weights);
PredictedStateVariance = dState*dState';
big_mat = [dObserved dState; H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)];
[~, mat] = qr2(big_mat,0);
mat = mat';
PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables, 1:number_of_observed_variables);
CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables);
StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables));
Error = y-PredictedObservedMean;
StateVectorMean = PredictedStateMean+(CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*Error;
if ParticleOptions.cpf_weights_method.amisanotristani
Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1), PredictedObservedVarianceSquareRoot, Error);
end
else
dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean);
dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean);
PredictedStateVariance = dState*diag(weights_c)*dState';
PredictedObservedVariance = dObserved*diag(weights_c)*dObserved'+H;
PredictedStateAndObservedCovariance = dState*diag(weights_c)*dObserved';
KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance;
Error = y-PredictedObservedMean;
StateVectorMean = PredictedStateMean+KalmanFilterGain*Error;
StateVectorVariance = PredictedStateVariance-KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
StateVectorVariance = 0.5*(StateVectorVariance+StateVectorVariance');
[StateVectorVarianceSquareRoot, p] = chol(StateVectorVariance, 'lower') ;
if p
flag = true;
ProposalStateVector = zeros(number_of_state_variables, 1);
Weights = 0.0;
return
end
if ParticleOptions.cpf_weights_method.amisanotristani
Weights = SampleWeights.*probability2(zeros(number_of_observed_variables, 1), chol(PredictedObservedVariance)', Error);
end
end
ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot, 2), 1)+StateVectorMean;
if ParticleOptions.cpf_weights_method.murrayjonesparslow
PredictedStateVariance = 0.5*(PredictedStateVariance+PredictedStateVariance');
[PredictedStateVarianceSquareRoot, p] = chol(PredictedStateVariance, 'lower');
if p
flag = true;
ProposalStateVector = zeros(number_of_state_variables,1);
Weights = 0.0;
return
end
Prior = probability2(PredictedStateMean, PredictedStateVarianceSquareRoot, ProposalStateVector);
Posterior = probability2(StateVectorMean, StateVectorVarianceSquareRoot, ProposalStateVector);
Likelihood = probability2(y, H_lower_triangular_cholesky, measurement_equations(ProposalStateVector, ReducedForm, ThreadsOptions, DynareOptions, Model));
Weights = SampleWeights.*Likelihood.*(Prior./Posterior);
end

110
matlab/nonlinear-filters/conditional_particle_filter.m Normal file
View File

@ -0,0 +1,110 @@
function [LIK,lik] = conditional_particle_filter(ReducedForm, Y, s, ParticleOptions, ThreadsOptions, DynareOptions, Model)
% Evaluates the likelihood of a non-linear model with a particle filter
%
% INPUTS
% - ReducedForm [structure] Matlab's structure describing the reduced form model.
% - Y [double] p×T matrix of (detrended) data, where p is the number of observed variables.
% - s [integer] scalar, likelihood evaluation starts at s (has to be smaller than T, the sample length provided in Y).
% - ParticlesOptions [struct]
% - ThreadsOptions [struct]
% - DynareOptions [struct]
% - Model [struct]
%
% OUTPUTS
% - LIK [double] scalar, likelihood
% - lik [double] (T-s+1)×1 vector, density of observations in each period.
%
% REMARKS
% - The proposal is built using the Kalman updating step for each particle.
% - we need draws in the errors distributions
% Whether we use Monte-Carlo draws from a multivariate gaussian distribution
% as in Amisano & Tristani (JEDC 2010).
% Whether we use multidimensional Gaussian sparse grids approximations:
% - a univariate Kronrod-Paterson Gaussian quadrature combined by the Smolyak
% operator (ref: Winschel & Kratzig, 2010).
% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009a,2009b).
% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1997, van der
% Merwe & Wan 2003).
%
% Pros:
% - Allows using current observable information in the proposal
% - The use of sparse grids Gaussian approximation is much faster than the Monte-Carlo approach
% Cons:
% - The use of the Kalman updating step may biais the proposal distribution since
% it has been derived in a linear context and is implemented in a nonlinear
% context. That is why particle resampling is performed.
% Copyright © 2009-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Set default for third input argument.
if isempty(s)
s = 1;
end
T = size(Y,2);
p = length(ReducedForm.mf1);
n = ParticleOptions.number_of_particles;
% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
if isempty(H)
H = 0;
H_lower_triangular_cholesky = 0;
else
H_lower_triangular_cholesky = chol(H)';
end
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot, 2);
Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn().
set_dynare_seed('default');
% Initialization of the likelihood.
lik = NaN(T, 1);
ks = 0;
StateParticles = bsxfun(@plus, StateVectorVarianceSquareRoot*randn(state_variance_rank, n), StateVectorMean);
SampleWeights = ones(1, n)/n;
for t=1:T
flags = false(n, 1);
for i=1:n
[StateParticles(:,i), SampleWeights(i), flags(i)] = ...
conditional_filter_proposal(ReducedForm, Y(:,t), StateParticles(:,i), SampleWeights(i), Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ParticleOptions, ThreadsOptions, DynareOptions, Model);
end
if any(flags)
LIK = -Inf;
lik(t) = -Inf;
return
end
SumSampleWeights = sum(SampleWeights);
lik(t) = log(SumSampleWeights);
SampleWeights = SampleWeights./SumSampleWeights;
if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*T) || ParticleOptions.resampling.status.systematic
ks = ks + 1;
StateParticles = resample(StateParticles', SampleWeights', ParticleOptions)';
SampleWeights = ones(1, n)/n;
end
end
LIK = -sum(lik(s:end));

51
matlab/nonlinear-filters/fit_gaussian_mixture.m Normal file
View File

@ -0,0 +1,51 @@
function [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(X,X_weights,StateMu,StateSqrtP,StateWeights,crit,niters,check)
% Copyright © 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[dim,Ndata] = size(X);
M = size(StateMu,2) ;
if check % Ensure that covariances don't collapse
MIN_COVAR_SQRT = sqrt(eps);
init_covars = StateSqrtP;
end
eold = -Inf;
for n=1:niters
% Calculate posteriors based on old parameters
[prior,likelihood,marginal,posterior] = probability3(StateMu,StateSqrtP,StateWeights,X,X_weights);
e = sum(log(marginal));
if (n > 1 && abs((e - eold)/eold) < crit)
return;
else
eold = e;
end
new_pr = (sum(posterior,2))';
StateWeights = new_pr/Ndata;
StateMu = bsxfun(@rdivide,(posterior*X')',new_pr);
for j=1:M
diffs = bsxfun(@minus,X,StateMu(:,j));
tpost = (1/sqrt(new_pr(j)))*sqrt(posterior(j,:));
diffs = bsxfun(@times,diffs,tpost);
[foo,tcov] = qr2(diffs',0);
StateSqrtP(:,:,j) = tcov';
if check
if min(abs(diag(StateSqrtP(:,:,j)))) < MIN_COVAR_SQRT
StateSqrtP(:,:,j) = init_covars(:,:,j);
end
end
end
end

50
matlab/nonlinear-filters/gaussian_densities.m Normal file
View File

@ -0,0 +1,50 @@
function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions,DynareOptions, Model)
%
% Elements to calculate the importance sampling ratio
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
% reduced_form_model.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% smolyak_accuracy [integer] scalar.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright © 2009-2019 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% proposal density
proposal = probability2(mut_t, sqr_Pss_t_t, particles);
% prior density
prior = probability2(st_t_1, sqr_Pss_t_t_1, particles);
% likelihood
yt_t_1_i = measurement_equations(particles, ReducedForm, ThreadsOptions, DynareOptions, Model);
likelihood = probability2(obs, sqrt(H), yt_t_1_i);
IncrementalWeights = likelihood.*prior./proposal;

135
matlab/nonlinear-filters/gaussian_filter.m Normal file
View File

@ -0,0 +1,135 @@
function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model)
% Evaluates the likelihood of a non-linear model approximating the
% predictive (prior) and filtered (posterior) densities for state variables
% by gaussian distributions.
% Gaussian approximation is done by:
% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009).
% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1995)
% - Monte-Carlo draws from a multivariate gaussian distribution.
% First and second moments of prior and posterior state densities are computed
% from the resulting nodes/particles and allows to generate new distributions at the
% following observation.
% Pros: The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles
% filters since it treats a lesser number of particles. Furthermore, in all cases, there is no need
% of resampling.
% Cons: estimations may be biaised if the model is truly non-gaussian
% since predictive and filtered densities are unimodal.
%
% INPUTS
% Reduced_Form [structure] Matlab's structure describing the reduced form model.
% Y [double] matrix of original observed variables.
% start [double] structural parameters.
% ParticleOptions [structure] Matlab's structure describing options concerning particle filtering.
% ThreadsOptions [structure] Matlab's structure.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright © 2009-2019 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Set default
if isempty(start)
start = 1;
end
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_particles = ParticleOptions.number_of_particles;
% compute gaussian quadrature nodes and weights on states and shocks
if ParticleOptions.distribution_approximation.cubature
[nodes2, weights2] = spherical_radial_sigma_points(number_of_state_variables);
weights_c2 = weights2;
elseif ParticleOptions.distribution_approximation.unscented
[nodes2, weights2, weights_c2] = unscented_sigma_points(number_of_state_variables,ParticleOptions);
else
if ~ParticleOptions.distribution_approximation.montecarlo
error('This approximation for the proposal is unknown!')
end
end
if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default');
end
% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
if isempty(H)
H = 0;
H_lower_triangular_cholesky = 0;
else
H_lower_triangular_cholesky = reduced_rank_cholesky(H)';
end
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)';
% Initialization of the likelihood.
const_lik = (2*pi)^(number_of_observed_variables/2) ;
lik = NaN(sample_size,1);
LIK = NaN;
for t=1:sample_size
[PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm, Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ...
H, ParticleOptions, ThreadsOptions, DynareOptions, Model);
if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented
StateParticles = bsxfun(@plus, StateVectorMean, StateVectorVarianceSquareRoot*nodes2');
IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ...
PredictedStateVarianceSquareRoot, StateParticles, H, const_lik, ...
weights2, weights_c2, ReducedForm, ThreadsOptions, ...
DynareOptions, Model);
SampleWeights = weights2.*IncrementalWeights;
else
StateParticles = bsxfun(@plus, StateVectorVarianceSquareRoot*randn(state_variance_rank, number_of_particles), StateVectorMean) ;
IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik, ...
1/number_of_particles,1/number_of_particles,ReducedForm,ThreadsOptions, ...
DynareOptions, Model);
SampleWeights = IncrementalWeights/number_of_particles;
end
SampleWeights = SampleWeights + 1e-6*ones(size(SampleWeights, 1), 1);
SumSampleWeights = sum(SampleWeights);
lik(t) = log(SumSampleWeights);
SampleWeights = SampleWeights./SumSampleWeights;
if not(ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented)
if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
StateParticles = resample(StateParticles', SampleWeights, ParticleOptions)';
SampleWeights = ones(number_of_particles, 1)/number_of_particles;
end
end
StateVectorMean = StateParticles*SampleWeights;
temp = bsxfun(@minus, StateParticles, StateVectorMean);
StateVectorVarianceSquareRoot = reduced_rank_cholesky(bsxfun(@times,SampleWeights',temp)*temp')';
end
LIK = -sum(lik(start:end));

139
matlab/nonlinear-filters/gaussian_filter_bank.m Normal file
View File

@ -0,0 +1,139 @@
function [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm, obs, StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ...
ParticleOptions, ThreadsOptions, DynareOptions, Model)
%
% Computes the proposal with a gaussian approximation for importance
% sampling
% 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 © 2009-2022 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
order = DynareOptions.order;
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if order == 3
% Set local state space model (third order approximation).
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
end
constant = ReducedForm.constant;
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
if ParticleOptions.proposal_approximation.montecarlo
nodes = randn(ParticleOptions.number_of_particles, number_of_state_variables+number_of_structural_innovations) ;
weights = 1/ParticleOptions.number_of_particles ;
weights_c = weights ;
elseif ParticleOptions.proposal_approximation.cubature
[nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations) ;
weights_c = weights ;
elseif ParticleOptions.proposal_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations, ParticleOptions);
else
error('This approximation for the proposal is not implemented or unknown!')
end
xbar = [StateVectorMean ; zeros(number_of_structural_innovations,1)] ;
sqr_Px = [ StateVectorVarianceSquareRoot, zeros(number_of_state_variables, number_of_structural_innovations);
zeros(number_of_structural_innovations, number_of_state_variables) Q_lower_triangular_cholesky];
sigma_points = bsxfun(@plus, xbar, sqr_Px*(nodes'));
StateVectors = sigma_points(1:number_of_state_variables,:);
epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
else
if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end
PredictedStateMean = tmp(mf0,:)*weights;
PredictedObservedMean = tmp(mf1,:)*weights;
if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_approximation.montecarlo
PredictedStateMean = sum(PredictedStateMean, 2);
PredictedObservedMean = sum(PredictedObservedMean, 2);
dState = bsxfun(@minus,tmp(mf0,:), PredictedStateMean)'.*sqrt(weights);
dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean)'.*sqrt(weights);
PredictedStateVarianceSquareRoot = chol(dState'*dState)';
big_mat = [dObserved, dState ; H_lower_triangular_cholesky, zeros(number_of_observed_variables,number_of_state_variables)];
[~, mat] = qr2(big_mat, 0);
mat = mat';
PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables);
CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables);
StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables));
PredictionError = obs - PredictedObservedMean;
StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*PredictionError;
else
dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean);
dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean);
PredictedStateVariance = dState*diag(weights_c)*dState';
PredictedObservedVariance = dObserved*diag(weights_c)*dObserved' + H;
PredictedStateAndObservedCovariance = dState*diag(weights_c)*dObserved';
PredictedStateVarianceSquareRoot = chol(PredictedStateVariance)';
PredictionError = obs - PredictedObservedMean;
KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance;
StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError;
StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';
end

55
matlab/nonlinear-filters/gaussian_mixture_densities.m Normal file
View File

@ -0,0 +1,55 @@
function IncrementalWeights = gaussian_mixture_densities(obs, StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
ReducedForm, ThreadsOptions, DynareOptions, Model)
% Elements to calculate the importance sampling ratio
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
% reduced_form_model.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% smolyak_accuracy [integer] scalar.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright © 2009-2019 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Compute the density of particles under the prior distribution
[~, ~, prior] = probability(StateMuPrior, StateSqrtPPrior, StateWeightsPrior, StateParticles);
prior = prior';
% Compute the density of particles under the proposal distribution
[~, ~, proposal] = probability(StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles);
proposal = proposal';
% Compute the density of the current observation conditionally to each particle
yt_t_1_i = measurement_equations(StateParticles, ReducedForm, ThreadsOptions, DynareOptions, Model);
% likelihood
likelihood = probability2(obs, sqrt(H), yt_t_1_i);
IncrementalWeights = likelihood.*prior./proposal;

226
matlab/nonlinear-filters/gaussian_mixture_filter.m Normal file
View File

@ -0,0 +1,226 @@
function [LIK, lik] = gaussian_mixture_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model)
% Evaluates the likelihood of a non-linear model approximating the state
% variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing
% a wide variety of generalized distributions (when multimodal for instance).
% Each gaussian distribution is obtained whether
% - with a radial-spherical cubature
% - with scaled unscented sigma-points
% A Sparse grid Kalman Filter is implemented on each component of the mixture,
% which confers it a weight about current information.
% Information on the current observables is then embodied in the proposal
% distribution in which we draw particles, which allows
% - reaching a greater precision relatively to a standard particle filter,
% - reducing the number of particles needed,
% - still being faster.
%
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
% reduced_form_model.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% Van der Meerwe & Wan, Gaussian Mixture Sigma-Point Particle Filters for Sequential
% Probabilistic Inference in Dynamic State-Space Models.
% Heiss & Winschel, 2010, Journal of Applied Economics.
% Winschel & Kratzig, 2010, Econometrica.
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright © 2009-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Set default
if isempty(start)
start = 1;
end
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
G = ParticleOptions.mixture_state_variables; % number of GM components in state
number_of_particles = ParticleOptions.number_of_particles;
% compute gaussian quadrature nodes and weights on states and shocks
if ParticleOptions.distribution_approximation.cubature
[nodes, weights] = spherical_radial_sigma_points(number_of_state_variables);
elseif ParticleOptions.distribution_approximation.unscented
[nodes, weights] = unscented_sigma_points(number_of_state_variables, ParticleOptions);
else
if ~ParticleOptions.distribution_approximation.montecarlo
error('This approximation for the proposal is unknown!')
end
end
if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default');
end
% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
if isempty(H)
H = 0;
H_lower_triangular_cholesky = 0;
else
H_lower_triangular_cholesky = reduced_rank_cholesky(H)';
end
Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)';
% Initialize mixtures
StateWeights = ones(1, G)/G;
StateMu = ReducedForm.StateVectorMean;
StateSqrtP = zeros(number_of_state_variables, number_of_state_variables, G);
temp = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
StateMu = bsxfun(@plus, StateMu, bsxfun(@times,diag(temp), (-(G-1)/2:1:(G-1)/2))/10);
for g=1:G
StateSqrtP(:,:,g) = temp/sqrt(G) ;
end
if ~ParticleOptions.mixture_structural_shocks
StructuralShocksMu = zeros(1, number_of_structural_innovations);
StructuralShocksWeights = 1;
I = 1;
StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu';
StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I);
StructuralShocksSqrtP(:,:,1) = Q_lower_triangular_cholesky;
elseif ParticleOptions.mixture_structural_shocks==1
if ParticleOptions.proposal_approximation.cubature
[StructuralShocksMu, StructuralShocksWeights] = spherical_radial_sigma_points(number_of_structural_innovations);
StructuralShocksWeights = ones(size(StructuralShocksMu, 1), 1)*StructuralShocksWeights;
elseif ParticleOptions.proposal_approximation.unscented
[StructuralShocksMu, StructuralShocksWeights] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions);
else
if ~ParticleOptions.distribution_approximation.montecarlo
error('This approximation for the proposal is unknown!')
end
end
I = size(StructuralShocksWeights, 1);
StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu';
StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I);
for i=1:I
StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky;
end
else
if ParticleOptions.proposal_approximation.cubature
[StructuralShocksMu, StructuralShocksWeights] = spherical_radial_sigma_points(number_of_structural_innovations);
StructuralShocksWeights = ones(size(StructuralShocksMu, 1), 1)*StructuralShocksWeights ;
elseif ParticleOptions.proposal_approximation.unscented
[StructuralShocksMu, StructuralShocksWeights] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions);
else
if ~ParticleOptions.distribution_approximation.montecarlo
error('This approximation for the proposal is unknown!')
end
end
I = size(StructuralShocksWeights, 1);
StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu';
StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I);
for i=1:I
StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky/sqrt(StructuralShocksWeights(i));
end
end
ObservationShocksWeights = 1;
J = 1 ;
Gprime = G*I;
Gsecond = G*I*J;
SampleWeights = ones(Gsecond, 1)/Gsecond;
StateWeightsPrior = zeros(1,Gprime);
StateMuPrior = zeros(number_of_state_variables,Gprime);
StateSqrtPPrior = zeros(number_of_state_variables, number_of_state_variables, Gprime);
StateWeightsPost = zeros(1, Gsecond);
StateMuPost = zeros(number_of_state_variables, Gsecond);
StateSqrtPPost = zeros(number_of_state_variables, number_of_state_variables, Gsecond);
const_lik = (2*pi)^(.5*number_of_observed_variables);
lik = NaN(sample_size, 1);
LIK = NaN;
for t=1:sample_size
% Build the proposal joint quadratures of Gaussian on states, structural
% shocks and observation shocks based on each combination of mixtures
for i=1:I
for j=1:J
for g=1:G
gprime = g + (i-1)*G;
gsecond = gprime + (j-1)*Gprime;
[StateMuPrior(:,gprime), StateSqrtPPrior(:,:,gprime), StateWeightsPrior(1,gprime), ...
StateMuPost(:,gsecond), StateSqrtPPost(:,:,gsecond), StateWeightsPost(1,gsecond)] = ...
gaussian_mixture_filter_bank(ReducedForm,Y(:,t), StateMu(:,g), StateSqrtP(:,:,g), StateWeights(g),...
StructuralShocksMu(:,i), StructuralShocksSqrtP(:,:,i), StructuralShocksWeights(i),...
ObservationShocksWeights(j), H, H_lower_triangular_cholesky, const_lik, ...
ParticleOptions, ThreadsOptions, DynareOptions, Model);
end
end
end
% Normalize weights
StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior, 2);
StateWeightsPost = StateWeightsPost/sum(StateWeightsPost, 2);
if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented
for i=1:Gsecond
StateParticles = bsxfun(@plus, StateMuPost(:,i), StateSqrtPPost(:,:,i)*nodes');
IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
ReducedForm, ThreadsOptions, DynareOptions, Model);
SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights);
end
SumSampleWeights = sum(SampleWeights);
lik(t) = log(SumSampleWeights);
SampleWeights = SampleWeights./SumSampleWeights;
[~, SortedRandomIndx] = sort(rand(1,Gsecond));
SortedRandomIndx = SortedRandomIndx(1:G);
indx = resample(0,SampleWeights,ParticleOptions);
indx = indx(SortedRandomIndx);
StateMu = StateMuPost(:,indx);
StateSqrtP = StateSqrtPPost(:,:,indx);
StateWeights = ones(1,G)/G;
else
% Sample particle in the proposal distribution, ie the posterior state GM
StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles);
IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
ReducedForm, ThreadsOptions, DynareOptions, Model);
SampleWeights = IncrementalWeights/number_of_particles;
SumSampleWeights = sum(SampleWeights,1);
SampleWeights = SampleWeights./SumSampleWeights;
lik(t) = log(SumSampleWeights);
if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
StateParticles = resample(StateParticles',SampleWeights',ParticleOptions)';
SampleWeights = ones(number_of_particles,1)/number_of_particles;
end
[StateMu, StateSqrtP, StateWeights] = fit_gaussian_mixture(StateParticles, SampleWeights', StateMu, StateSqrtP, StateWeights, 0.001, 10, 1);
end
end
LIK = -sum(lik(start:end));

146
matlab/nonlinear-filters/gaussian_mixture_filter_bank.m Normal file
View File

@ -0,0 +1,146 @@
function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPPost,StateWeightsPost] =...
gaussian_mixture_filter_bank(ReducedForm, obs, StateMu, StateSqrtP, StateWeights, ...
StructuralShocksMu, StructuralShocksSqrtP, StructuralShocksWeights, ...
ObservationShocksWeights, H, H_lower_triangular_cholesky, normfactO, ...
ParticleOptions, ThreadsOptions, DynareOptions, Model)
% 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 © 2009-2022 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
order = DynareOptions.order;
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if order == 3
% Set local state space model (third order approximation).
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
end
constant = ReducedForm.constant;
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
numb = number_of_state_variables+number_of_structural_innovations;
if ParticleOptions.proposal_approximation.cubature
[nodes3, weights3] = spherical_radial_sigma_points(numb);
weights_c3 = weights3;
elseif ParticleOptions.proposal_approximation.unscented
[nodes3, weights3, weights_c3] = unscented_sigma_points(numb, ParticleOptions);
else
error('This approximation for the proposal is unknown!')
end
epsilon = bsxfun(@plus, StructuralShocksSqrtP*nodes3(:,number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations)', StructuralShocksMu);
StateVectors = bsxfun(@plus, StateSqrtP*nodes3(:,1:number_of_state_variables)', StateMu);
yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
else
if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end
PredictedStateMean = tmp(mf0,:)*weights3;
PredictedObservedMean = tmp(mf1,:)*weights3;
if ParticleOptions.proposal_approximation.cubature
PredictedStateMean = sum(PredictedStateMean, 2);
PredictedObservedMean = sum(PredictedObservedMean, 2);
dState = (bsxfun(@minus, tmp(mf0,:), PredictedStateMean)').*sqrt(weights3);
dObserved = (bsxfun(@minus, tmp(mf1,:), PredictedObservedMean)').*sqrt(weights3);
PredictedStateVariance = dState'*dState;
big_mat = [dObserved, dState ; H_lower_triangular_cholesky, zeros(number_of_observed_variables, number_of_state_variables)];
[~, mat] = qr2(big_mat, 0);
mat = mat';
PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables, 1:number_of_observed_variables);
CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables), 1:number_of_observed_variables);
StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables), number_of_observed_variables+(1:number_of_state_variables));
iPredictedObservedVarianceSquareRoot = inv(PredictedObservedVarianceSquareRoot);
iPredictedObservedVariance = iPredictedObservedVarianceSquareRoot'*iPredictedObservedVarianceSquareRoot;
sqrdet = 1/sqrt(det(iPredictedObservedVariance));
PredictionError = obs - PredictedObservedMean;
StateVectorMean = PredictedStateMean + CovarianceObservedStateSquareRoot*iPredictedObservedVarianceSquareRoot*PredictionError;
else
dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean);
dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean);
PredictedStateVariance = dState*diag(weights_c3)*dState';
PredictedObservedVariance = dObserved*diag(weights_c3)*dObserved' + H;
PredictedStateAndObservedCovariance = dState*diag(weights_c3)*dObserved';
sqrdet = sqrt(det(PredictedObservedVariance));
iPredictedObservedVariance = inv(PredictedObservedVariance);
PredictionError = obs - PredictedObservedMean;
KalmanFilterGain = PredictedStateAndObservedCovariance*iPredictedObservedVariance;
StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError;
StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';
end
data_lik_GM_g = exp(-0.5*PredictionError'*iPredictedObservedVariance*PredictionError)/abs(normfactO*sqrdet) + 1e-99;
StateMuPrior = PredictedStateMean;
StateSqrtPPrior = reduced_rank_cholesky(PredictedStateVariance)';
StateWeightsPrior = StateWeights*StructuralShocksWeights;
StateMuPost = StateVectorMean;
StateSqrtPPost = StateVectorVarianceSquareRoot;
StateWeightsPost = StateWeightsPrior*ObservationShocksWeights*data_lik_GM_g;

28
matlab/nonlinear-filters/importance_sampling.m Normal file
View File

@ -0,0 +1,28 @@
function State_Particles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost,numP)
% Copyright © 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[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);

57
matlab/nonlinear-filters/measurement_equations.m Normal file
View File

@ -0,0 +1,57 @@
function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions, DynareOptions, Model)
% Copyright © 2013-2022 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
order = DynareOptions.order;
mf1 = ReducedForm.mf1;
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
ghx = ReducedForm.ghx(mf1,:);
ghu = ReducedForm.ghu(mf1,:);
ghxx = ReducedForm.ghxx(mf1,:);
ghuu = ReducedForm.ghuu(mf1,:);
ghxu = ReducedForm.ghxu(mf1,:);
ghs2 = ReducedForm.ghs2(mf1,:);
if order == 3
ghxxx = ReducedForm.ghxxx(mf1,:);
ghuuu = ReducedForm.ghuuu(mf1,:);
ghxxu = ReducedForm.ghxxu(mf1,:);
ghxuu = ReducedForm.ghxuu(mf1,:);
ghxss = ReducedForm.ghxss(mf1,:);
ghuss = ReducedForm.ghuss(mf1,:);
end
end
steadystate = ReducedForm.steadystate(mf1,:);
constant = ReducedForm.constant(mf1,:);
state_variables_steady_state = ReducedForm.state_variables_steady_state;
number_of_structural_innovations = length(ReducedForm.Q);
yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions, udr);
measure = tmp(mf1,:);
else
if order == 2
measure = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, size(yhat,2)), ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
measure = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, size(yhat,2)), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end

72
matlab/nonlinear-filters/multivariate_smooth_resampling.m Normal file
View File

@ -0,0 +1,72 @@
function new_particles = multivariate_smooth_resampling(particles,weights)
% Smooth Resampling of the particles.
%@info:
%! @deftypefn {Function File} {@var{new_particles} =} multivariate_smooth_resampling (@var{weights}, @var{particles}, @var{number_of_new_particles}, @var{number_of_partitions})
%! @anchor{particle/multivariate_smooth_resampling}
%! @sp 1
%! Smooth Resampling of the particles (multivariate version).
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item weights
%! n*1 vector of doubles, particles' weights.
%! @item particles
%! n*1 vector of doubles, particles.
%! @item number_of_new_particles
%! Integer scalar.
%! @item number_of_partitions
%! Integer scalar.
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item indx
%! number_of_new_particles*1 vector of doubles, new particles.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{particle/sequantial_importance_particle_filter}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{particle/univariate_smooth_resampling}
%! @sp 2
%! @end deftypefn
%@eod:
% Copyright © 2012-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% 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') ;

53
matlab/nonlinear-filters/mykmeans.m Normal file
View File

@ -0,0 +1,53 @@
function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod)
% Copyright © 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[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

21
matlab/nonlinear-filters/neff.m Normal file
View File

@ -0,0 +1,21 @@
function n = neff(w)
% Evaluates the criterion for resampling
% Copyright © 2013-2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
n = 1/dot(w,w);

182
matlab/nonlinear-filters/nonlinear_kalman_filter.m Normal file
View File

@ -0,0 +1,182 @@
function [LIK,lik] = nonlinear_kalman_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model)
% Evaluates the likelihood of a non-linear model approximating the
% predictive (prior) and filtered (posterior) densities for state variables
% by a Kalman filter.
% Gaussian distribution approximation is done by:
% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009).
% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1995)
% - Monte-Carlo draws from a multivariate gaussian distribution.
% First and second moments of prior and posterior state densities are computed
% from the resulting nodes/particles and allows to generate new distributions at the
% following observation.
% Pros: The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles
% filters since it treats a lesser number of particles.
% Cons: 1. Application a linear projection formulae in a nonlinear context.
% 2. Parameter estimations may be biaised if the model is truly non-gaussian since predictive and
% filtered densities are unimodal.
%
% INPUTS
% Reduced_Form [structure] Matlab's structure describing the reduced form model.
% Y [double] matrix of original observed variables.
% start [double] structural parameters.
% ParticleOptions [structure] Matlab's structure describing options concerning particle filtering.
% ThreadsOptions [structure] Matlab's structure.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright © 2009-2022 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Set default
if isempty(start)
start = 1;
end
order = DynareOptions.order;
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if (order == 3)
% Set local state space model (third order approximation).
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
end
constant = ReducedForm.constant;
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
% compute gaussian quadrature nodes and weights on states and shocks
if ParticleOptions.proposal_approximation.montecarlo
nodes = randn(ParticleOptions.number_of_particles,number_of_state_variables+number_of_structural_innovations);
weights = 1/ParticleOptions.number_of_particles;
weights_c = weights;
elseif ParticleOptions.proposal_approximation.cubature
[nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations);
weights_c = weights;
elseif ParticleOptions.proposal_approximation.unscented
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,ParticleOptions);
else
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
if ParticleOptions.distribution_approximation.montecarlo
set_dynare_seed('default');
end
% Get covariance matrices
H = ReducedForm.H;
H_lower_triangular_cholesky = chol(H)' ;
Q_lower_triangular_cholesky = chol(ReducedForm.Q)';
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
% Initialization of the likelihood.
lik = NaN(sample_size,1);
LIK = NaN;
for t=1:sample_size
xbar = [StateVectorMean ; zeros(number_of_structural_innovations,1) ] ;
sqr_Px = [StateVectorVarianceSquareRoot zeros(number_of_state_variables,number_of_structural_innovations);
zeros(number_of_structural_innovations,number_of_state_variables) Q_lower_triangular_cholesky];
sigma_points = bsxfun(@plus,xbar,sqr_Px*(nodes'));
StateVectors = sigma_points(1:number_of_state_variables,:);
epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
else
if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
end
end
PredictedStateMean = tmp(mf0,:)*weights ;
PredictedObservedMean = tmp(mf1,:)*weights;
if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_approximation.montecarlo
PredictedStateMean = sum(PredictedStateMean,2);
PredictedObservedMean = sum(PredictedObservedMean,2);
dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean)'.*sqrt(weights);
dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean)'.*sqrt(weights);
big_mat = [dObserved dState ; [H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)] ];
[~, mat] = qr2(big_mat,0);
mat = mat';
PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables);
CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables);
StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables));
PredictionError = Y(:,t) - PredictedObservedMean;
StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*PredictionError;
else
dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean);
dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
PredictedStateVariance = dState*diag(weights_c)*dState';
PredictedObservedVariance = dObserved*diag(weights_c)*dObserved' + H;
PredictedStateAndObservedCovariance = dState*diag(weights_c)*dObserved';
PredictionError = Y(:,t) - PredictedObservedMean;
KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance;
StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError;
StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
[StateVectorVarianceSquareRoot, p]= chol(StateVectorVariance,'lower');
if p
LIK=-Inf;
lik(t)=-Inf;
return
end
[~, p]= chol(PredictedObservedVariance,'lower');
if p
LIK=-Inf;
lik(t)=-Inf;
return
end
end
lik(t) = log( sum(probability2(Y(:,t),H_lower_triangular_cholesky,tmp(mf1,:)).*weights,1) ) ;
end
LIK = -sum(lik(start:end));

455
matlab/nonlinear-filters/online_auxiliary_filter.m Normal file
View File

@ -0,0 +1,455 @@
function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, DynareResults)
% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
%
% INPUTS
% - xparam1 [double] n×1 vector, Initial condition for the estimated parameters.
% - DynareDataset [dseries] Sample used for estimation.
% - dataset_info [struct] Description of the sample.
% - DynareOptions [struct] Option values (options_).
% - Model [struct] Description of the model (M_).
% - EstimatedParameters [struct] Description of the estimated parameters (estim_params_).
% - BayesInfo [struct] Prior definition (bayestopt_).
% - DynareResults [struct] Results (oo_).
%
% OUTPUTS
% - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters).
% - pmode [double] n×1 vector, mode of the particles at the end of the sample (for the parameters).
% - pmedian [double] n×1 vector, median of the particles at the end of the sample (for the parameters).
% - pstdev [double] n×1 vector, st. dev. of the particles at the end of the sample (for the parameters).
% - p025 [double] n×1 vector, 2.5 percent of the particles are below p025(i) for i=1,…,n.
% - p975 [double] n×1 vector, 97.5 percent of the particles are below p975(i) for i=1,…,n.
% - covariance [double] n×n matrix, covariance of the particles at the end of the sample.
% Copyright © 2013-2022 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Set seed for randn().
set_dynare_seed('default');
pruning = DynareOptions.particle.pruning;
second_resample = DynareOptions.particle.resampling.status.systematic;
variance_update = true;
bounds = prior_bounds(BayesInfo, DynareOptions.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
% initialization of state particles
[~, Model, DynareOptions, DynareResults, ReducedForm] = solve_model_for_online_filter(true, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
order = DynareOptions.order;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_particles = DynareOptions.particle.number_of_particles;
number_of_parameters = size(xparam1,1);
Y = DynareDataset.data;
sample_size = size(Y,1);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
liu_west_delta = DynareOptions.particle.liu_west_delta;
% Get initial conditions for the state particles
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
if order == 2
StateVectors_ = StateVectors;
elseif order == 3
StateVectors_ = repmat(StateVectors,3,1);
else
error('Pruning is not available for orders > 3');
end
end
% parameters for the Liu & West filter
small_a = (3*liu_west_delta-1)/(2*liu_west_delta);
b_square = 1-small_a*small_a;
% Initialization of parameter particles
xparam = zeros(number_of_parameters,number_of_particles);
Prior = dprior(BayesInfo, DynareOptions.prior_trunc);
for i=1:number_of_particles
info = 12042009;
while info
candidate = Prior.draw();
[info, Model, DynareOptions, DynareResults] = solve_model_for_online_filter(false, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
if ~info
xparam(:,i) = candidate(:);
end
end
end
% Initialization of the weights of particles.
weights = ones(1,number_of_particles)/number_of_particles;
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables;
mean_xparam = zeros(number_of_parameters,sample_size);
mode_xparam = zeros(number_of_parameters,sample_size);
median_xparam = zeros(number_of_parameters,sample_size);
std_xparam = zeros(number_of_parameters,sample_size);
lb95_xparam = zeros(number_of_parameters,sample_size);
ub95_xparam = zeros(number_of_parameters,sample_size);
%% The Online filter
for t=1:sample_size
if t>1
fprintf('\nSubsample with %s first observations.\n\n', int2str(t))
else
fprintf('\nSubsample with only the first observation.\n\n')
end
% Moments of parameters particles distribution
m_bar = xparam*(weights');
temp = bsxfun(@minus,xparam,m_bar);
sigma_bar = (bsxfun(@times,weights,temp))*(temp');
if variance_update
chol_sigma_bar = chol(b_square*sigma_bar)';
end
% Prediction (without shocks)
fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam);
tau_tilde = zeros(1,number_of_particles);
for i=1:number_of_particles
% model resolution
[info, Model, DynareOptions, DynareResults, ReducedForm] = ...
solve_model_for_online_filter(false, fore_xparam(:,i), DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
if ~info(1)
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set local state space model (second-order approximation).
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
steadystate = ReducedForm.steadystate;
constant = ReducedForm.constant;
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if (order == 3)
% Set local state space model (third order approximation).
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
if pruning
if order == 2
state_variables_steady_state_ = state_variables_steady_state;
elseif order == 3
state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
else
error('Pruning is not available for orders > 3');
end
end
end
% particle likelihood contribution
yhat = bsxfun(@minus, StateVectors(:,i), state_variables_steady_state);
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions, udr);
else
if pruning
yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state_);
if order == 2
[tmp, ~] = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2);
elseif order == 3
[tmp, tmp_] = local_state_space_iteration_3(yhat_, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning);
else
error('Pruning is not available for orders > 3');
end
else
if order == 2
tmp = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end
end
PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
% Replace Gaussian density with a Student density with 3 degrees of freedom for fat tails.
z = sum(PredictionError.*(ReducedForm.H\PredictionError), 1) ;
tau_tilde(i) = weights(i).*(tpdf(z, 3*ones(size(z)))+1e-99) ;
end
end
% particles selection
tau_tilde = tau_tilde/sum(tau_tilde);
indx = resample(0, tau_tilde', DynareOptions.particle);
StateVectors = StateVectors(:,indx);
xparam = fore_xparam(:,indx);
if pruning
StateVectors_ = StateVectors_(:,indx);
end
w_stage1 = weights(indx)./tau_tilde(indx);
% draw in the new distributions
wtilde = zeros(1, number_of_particles);
for i=1:number_of_particles
info = 12042009;
counter=0;
while info(1) && counter <DynareOptions.particle.liu_west_max_resampling_tries
counter=counter+1;
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
% model resolution for new parameters particles
[info, Model, DynareOptions, DynareResults, ReducedForm] = ...
solve_model_for_online_filter(false, candidate, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults) ;
if ~info(1)
xparam(:,i) = candidate ;
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set local state space model (second order approximation).
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
constant = ReducedForm.constant;
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if (order == 3)
% Set local state space model (third order approximation).
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
if pruning
if order == 2
state_variables_steady_state_ = state_variables_steady_state;
mf0_ = mf0;
elseif order == 3
state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
mf0_ = repmat(mf0,1,3);
mask2 = number_of_state_variables+1:2*number_of_state_variables;
mask3 = 2*number_of_state_variables+1:number_of_state_variables;
mf0_(mask2) = mf0_(mask2)+size(ghx,1);
mf0_(mask3) = mf0_(mask3)+2*size(ghx,1);
else
error('Pruning is not available for orders > 3');
end
end
end
% Get covariance matrices and structural shocks
epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations, 1);
% compute particles likelihood contribution
yhat = bsxfun(@minus,StateVectors(:,i), state_variables_steady_state);
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
else
if pruning
yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state_);
if order == 2
[tmp, tmp_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2);
elseif order == 3
[tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning);
else
error('Pruning is not available for orders > 3');
end
StateVectors_(:,i) = tmp_(mf0_,:);
else
if order == 2
tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end
end
StateVectors(:,i) = tmp(mf0,:);
PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1)));
end
end
if counter==DynareOptions.particle.liu_west_max_resampling_tries
fprintf('\nLiu & West particle filter: I haven''t been able to solve the model in %u tries.\n',DynareOptions.particle.liu_west_max_resampling_tries)
fprintf('Liu & West particle filter: The last error message was: %s\n',get_error_message(info))
fprintf('Liu & West particle filter: You can try to increase liu_west_max_resampling_tries, but most\n')
fprintf('Liu & West particle filter: likely there is an issue with the model.\n')
error('Liu & West particle filter: unable to solve the model.')
end
end
end
% normalization
weights = wtilde/sum(wtilde);
if variance_update && (neff(weights)<DynareOptions.particle.resampling.threshold*sample_size)
variance_update = false;
end
% final resampling (not advised)
if second_resample
[~, idmode] = max(weights);
mode_xparam(:,t) = xparam(:,idmode);
indx = resample(0, weights,DynareOptions.particle);
StateVectors = StateVectors(:,indx) ;
if pruning
StateVectors_ = StateVectors_(:,indx);
end
xparam = xparam(:,indx);
weights = ones(1, number_of_particles)/number_of_particles;
mean_xparam(:,t) = mean(xparam, 2);
mat_var_cov = bsxfun(@minus, xparam, mean_xparam(:,t));
mat_var_cov = (mat_var_cov*mat_var_cov')/(number_of_particles-1);
std_xparam(:,t) = sqrt(diag(mat_var_cov));
for i=1:number_of_parameters
temp = sortrows(xparam(i,:)');
lb95_xparam(i,t) = temp(0.025*number_of_particles);
median_xparam(i,t) = temp(0.5*number_of_particles);
ub95_xparam(i,t) = temp(0.975*number_of_particles);
end
end
if second_resample
[~, idmode] = max(weights);
mode_xparam(:,t) = xparam(:,idmode);
mean_xparam(:,t) = xparam*(weights');
mat_var_cov = bsxfun(@minus, xparam,mean_xparam(:,t));
mat_var_cov = mat_var_cov*(bsxfun(@times, mat_var_cov, weights)');
std_xparam(:,t) = sqrt(diag(mat_var_cov));
for i=1:number_of_parameters
temp = sortrows([xparam(i,:)' weights'], 1);
cumulated_weights = cumsum(temp(:,2));
pass1 = false;
pass2 = false;
pass3 = false;
for j=1:number_of_particles
if ~pass1 && cumulated_weights(j)>=0.025
lb95_xparam(i,t) = temp(j,1);
pass1 = true;
end
if ~pass2 && cumulated_weights(j)>=0.5
median_xparam(i,t) = temp(j,1);
pass2 = true;
end
if ~pass3 && cumulated_weights(j)>=0.975
ub95_xparam(i,t) = temp(j,1);
pass3 = true;
end
end
end
end
str = sprintf(' Lower Bound (95%%) \t Mean \t\t\t Upper Bound (95%%)');
for l=1:size(xparam,1)
str = sprintf('%s\n %5.4f \t\t %7.5f \t\t %5.4f', str, lb95_xparam(l,t), mean_xparam(l,t), ub95_xparam(l,t));
end
disp(str)
disp('')
end
pmean = xparam(:,sample_size);
pmode = mode_xparam(:,sample_size);
pstdev = std_xparam(:,sample_size) ;
p025 = lb95_xparam(:,sample_size) ;
p975 = ub95_xparam(:,sample_size) ;
pmedian = median_xparam(:,sample_size) ;
covariance = mat_var_cov;
%% Plot parameters trajectory
TeX = DynareOptions.TeX;
nr = ceil(sqrt(number_of_parameters)) ;
nc = floor(sqrt(number_of_parameters));
nbplt = 1 ;
if TeX
fidTeX = fopen([Model.fname '_param_traj.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
for plt = 1:nbplt
hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories');
for k=1:length(pmean)
subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
% Draw the surface for an interval containing 95% of the particles.
area(1:sample_size, ub95_xparam(k,:), 'FaceColor', [.9 .9 .9], 'BaseValue', min(lb95_xparam(k,:)));
hold on
area(1:sample_size, lb95_xparam(k,:), 'FaceColor', [1 1 1], 'BaseValue', min(lb95_xparam(k,:)));
% Draw the mean of particles.
plot(1:sample_size, mean_xparam(k,:), '-k', 'linewidth', 2)
if TeX
title(texname,'interpreter','latex')
else
title(name,'interpreter','none')
end
hold off
axis tight
drawnow
end
dyn_saveas(hh, [Model.fname '_param_traj' int2str(plt)], DynareOptions.nodisplay, DynareOptions.graph_format);
if TeX
% TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParamTraj%s}\n',Model.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Parameters trajectories.}');
fprintf(fidTeX,'\\label{Fig:ParametersPlots:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
end
% Plot Parameter Densities
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation.
for plt = 1:nbplt
hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities');
for k=1:length(pmean)
subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ...
number_of_particles, optimal_bandwidth, kernel_function);
plot(density(:,1), density(:,2));
hold on
if TeX
title(texname,'interpreter','latex')
else
title(name,'interpreter','none')
end
hold off
axis tight
drawnow
end
dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
if TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format)))
% TeX eps loader file
fprintf(fidTeX, '\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/nc,1),M_.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Parameter densities based on the Liu/West particle filter.}');
fprintf(fidTeX,'\\label{Fig:ParameterDensities:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
end

39
matlab/nonlinear-filters/probability.m Normal file
View File

@ -0,0 +1,39 @@
function [prior,likelihood,C,posterior] = probability(mu,sqrtP,prior,X)
% Copyright © 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[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

37
matlab/nonlinear-filters/probability2.m Normal file
View File

@ -0,0 +1,37 @@
function [density] = probability2(mu,S,X)
%
% Multivariate gaussian density
%
% INPUTS
% n [integer] scalar, number of variables.
%
% OUTPUTS
% nodes [double] nodes of the cubature
% weigths [double] associated weigths
%
% REFERENCES
%
%
% NOTES
%
% Copyright © 2009-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
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 ;

39
matlab/nonlinear-filters/probability3.m Normal file
View File

@ -0,0 +1,39 @@
function [prior,likelihood,C,posterior] = probability3(mu,sqrtP,prior,X,X_weights)
% Copyright © 2013-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[dim,nov] = size(X);
M = size(mu,2) ;
if nargout>1
likelihood = zeros(M,nov);
normfact = (2*pi)^(dim/2);
for k=1:M
XX = bsxfun(@minus,X,mu(:,k));
S = sqrtP(:,:,k);
foo = S \ XX;
likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S))));
end
end
wlikelihood = bsxfun(@times,X_weights,likelihood) + 1e-99;
if nargout>2
C = prior*wlikelihood + 1e-99;
end
if nargout>3
posterior = bsxfun(@rdivide,bsxfun(@times,prior',wlikelihood),C) + 1e-99 ;
posterior = bsxfun(@rdivide,posterior,sum(posterior,1));
end

79
matlab/nonlinear-filters/resample.m Normal file
View File

@ -0,0 +1,79 @@
function resampled_output = resample(particles,weights,ParticleOptions)
% Resamples particles.
% if particles = 0, returns the resampling index (except for smooth resampling)
% Otherwise, returns the resampled particles set.
%@info:
%! @deftypefn {Function File} {@var{indx} =} resample (@var{weights}, @var{method})
%! @anchor{particle/resample}
%! @sp 1
%! Resamples particles.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item weights
%! n*1 vector of doubles, particles' weights.
%! @item method
%! string equal to 'residual' or 'traditional'.
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item indx
%! n*1 vector of intergers, indices.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{particle/sequantial_importance_particle_filter}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{residual_resampling}, @ref{traditional_resampling}
%! @sp 2
%! @end deftypefn
%@eod:
% Copyright © 2011-2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
defaultmethod = 1; % For residual based method set this variable equal to 0.
if defaultmethod
if ParticleOptions.resampling.method.kitagawa
resampled_output = traditional_resampling(particles,weights,rand);
elseif ParticleOptions.resampling.method.stratified
resampled_output = traditional_resampling(particles,weights,rand(size(weights)));
elseif ParticleOptions.resampling.method.smooth
if particles==0
error('Particle = 0 is incompatible with this resampling method!')
end
resampled_output = multivariate_smooth_resampling(particles,weights);
else
error('Unknown sampling method!')
end
else
if ParticleOptions.resampling.method.kitagawa
resampled_output = residual_resampling(particles,weights,rand);
elseif ParticleOptions.resampling.method.stratified
resampled_output = residual_resampling(particles,weights,rand(size(weights)));
else
error('Unknown sampling method!')
end
end

143
matlab/nonlinear-filters/residual_resampling.m Normal file
View File

@ -0,0 +1,143 @@
function return_resample = residual_resampling(particles,weights,noise)
% Resamples particles.
%@info:
%! @deftypefn {Function File} {@var{indx} =} residual_resampling (@var{weights})
%! @anchor{particle/residual_resampling}
%! @sp 1
%! Resamples particles.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item weights
%! n*1 vector of doubles, particles' weights.
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item indx
%! n*1 vector of intergers, indices.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{particle/resample}
%! @sp 2
%! @strong{This function calls:}
%! @sp 2
%! @end deftypefn
%@eod:
% Copyright © 2011-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% 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

185
matlab/nonlinear-filters/sequential_importance_particle_filter.m Normal file
View File

@ -0,0 +1,185 @@
function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model)
% Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
% Copyright © 2011-2022 Dynare Team
%
% This file is part of Dynare (particles module).
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare particles module is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
persistent init_flag
persistent mf0 mf1
persistent number_of_particles number_of_state_variables
persistent sample_size number_of_observed_variables number_of_structural_innovations
% Set default value for start
if isempty(start)
start = 1;
end
% Set flag for prunning
pruning = ParticleOptions.pruning;
% Get steady state and mean.
steadystate = ReducedForm.steadystate;
constant = ReducedForm.constant;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
order = DynareOptions.order;
% Set persistent variables (if needed).
if isempty(init_flag)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
number_of_particles = ParticleOptions.number_of_particles;
init_flag = 1;
end
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
udr = ReducedForm.udr;
else
% Set local state space model (first order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
ghs2 = ReducedForm.ghs2;
if order == 3
% Set local state space model (third order approximation).
ghxxx = ReducedForm.ghxxx;
ghuuu = ReducedForm.ghuuu;
ghxxu = ReducedForm.ghxxu;
ghxuu = ReducedForm.ghxuu;
ghxss = ReducedForm.ghxss;
ghuss = ReducedForm.ghuss;
end
end
% Get covariance matrices.
Q = ReducedForm.Q; % Covariance matrix of the structural innovations.
H = ReducedForm.H; % Covariance matrix of the measurement errors.
if isempty(H)
H = 0;
end
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables +log(det(H)) ;
lik = NaN(sample_size,1);
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';%reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
if pruning
StateVectorMean_ = StateVectorMean;
StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
end
% Get the rank of StateVectorVarianceSquareRoot
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
% Factorize the covariance matrix of the structural innovations
Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn().
set_dynare_seed('default');
% Initialization of the weights across particles.
weights = ones(1,number_of_particles)/number_of_particles ;
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
if order == 2
StateVectors_ = StateVectors;
state_variables_steady_state_ = state_variables_steady_state;
mf0_ = mf0;
elseif order == 3
StateVectors_ = repmat(StateVectors,3,1);
state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
mf0_ = repmat(mf0,1,3);
mask2 = number_of_state_variables+1:2*number_of_state_variables;
mask3 = 2*number_of_state_variables+1:3*number_of_state_variables;
mf0_(mask2) = mf0_(mask2)+size(ghx,1);
mf0_(mask3) = mf0_(mask3)+2*size(ghx,1);
else
error('Pruning is not available for orders > 3');
end
end
% Loop over observations
for t=1:sample_size
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
if pruning
yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state_);
if order == 2
[tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
[tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
else
error('Pruning is not available for orders > 3');
end
else
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
else
if order == 2
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
elseif order == 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
else
error('Order > 3: use_k_order_solver should be set to true');
end
end
end
%PredictedObservedMean = tmp(mf1,:)*transpose(weights);
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
%dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
%PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H;
%PredictedObservedVariance = H;
if rcond(H) > 1e-16
lnw = -.5*(const_lik+sum(PredictionError.*(H\PredictionError),1));
else
LIK = NaN;
return
end
dfac = max(lnw);
wtilde = weights.*exp(lnw-dfac);
lik(t) = log(sum(wtilde))+dfac;
weights = wtilde/sum(wtilde);
if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
if pruning
temp = resample([tmp(mf0,:)' tmp_(mf0_,:)'],weights',ParticleOptions);
StateVectors = temp(:,1:number_of_state_variables)';
StateVectors_ = temp(:,number_of_state_variables+1:end)';
else
StateVectors = resample(tmp(mf0,:)',weights',ParticleOptions)';
end
weights = ones(1,number_of_particles)/number_of_particles;
elseif ParticleOptions.resampling.status.none
StateVectors = tmp(mf0,:);
if pruning
StateVectors_ = tmp_(mf0_,:);
end
end
end
LIK = -sum(lik(start:end));

215
matlab/nonlinear-filters/solve_model_for_online_filter.m Normal file
View File

@ -0,0 +1,215 @@
function [info, Model, DynareOptions, DynareResults, ReducedForm] = ...
solve_model_for_online_filter(setinitialcondition, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults)
% Solves the dsge model for an particular parameters set.
%
% INPUTS
% - setinitialcondition [logical] return initial condition if true.
% - xparam1 [double] n×1 vector, parameter values.
% - DynareDataset [struct] Dataset for estimation (dataset_).
% - DynareOptions [struct] Dynare options (options_).
% - Model [struct] Model description (M_).
% - EstimatedParameters [struct] Estimated parameters (estim_params_).
% - BayesInfo [struct] Prior definition (bayestopt_).
% - DynareResults [struct] Dynare results (oo_).
%
% OUTPUTS
% - info [integer] scalar, nonzero if any problem occur when computing the reduced form.
% - Model [struct] Model description (M_).
% - DynareOptions [struct] Dynare options (options_).
% - DynareResults [struct] Dynare results (oo_).
% - ReducedForm [struct] Reduced form model.
% Copyright © 2013-2022 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
persistent init_flag restrict_variables_idx state_variables_idx mf0 mf1 number_of_state_variables
info = 0;
%----------------------------------------------------
% 1. Get the structural parameters & define penalties
%----------------------------------------------------
% Test if some parameters are smaller than the lower bound of the prior domain.
if any(xparam1<bounds.lb)
info = 41;
return
end
% Test if some parameters are greater than the upper bound of the prior domain.
if any(xparam1>bounds.ub)
info = 42;
return
end
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
Q = Model.Sigma_e;
H = Model.H;
for i=1:EstimatedParameters.nvx
k =EstimatedParameters.var_exo(i,1);
Q(k,k) = xparam1(i)*xparam1(i);
end
offset = EstimatedParameters.nvx;
if EstimatedParameters.nvn
for i=1:EstimatedParameters.nvn
H(i,i) = xparam1(i+offset)*xparam1(i+offset);
end
offset = offset+EstimatedParameters.nvn;
else
H = zeros(size(DynareDataset.data, 2));
end
% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
if EstimatedParameters.ncx
for i=1:EstimatedParameters.ncx
k1 =EstimatedParameters.corrx(i,1);
k2 =EstimatedParameters.corrx(i,2);
Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
Q(k2,k1) = Q(k1,k2);
end
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
[~, testQ] = chol(Q);
if testQ
% The variance-covariance matrix of the structural innovations is not definite positive.
info = 43;
return
end
offset = offset+EstimatedParameters.ncx;
end
% Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite.
if EstimatedParameters.ncn
corrn_observable_correspondence = EstimatedParameters.corrn_observable_correspondence;
for i=1:EstimatedParameters.ncn
k1 = corrn_observable_correspondence(i,1);
k2 = corrn_observable_correspondence(i,2);
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
H(k2,k1) = H(k1,k2);
end
% Try to compute the cholesky decomposition of H (possible iff H is positive definite)
[~, testH] = chol(H);
if testH
% The variance-covariance matrix of the measurement errors is not definite positive.
info = 44;
return
end
offset = offset+EstimatedParameters.ncn;
end
% Update estimated structural parameters in Mode.params.
if EstimatedParameters.np > 0
Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end);
end
% Update Model.Sigma_e and Model.H.
Model.Sigma_e = Q;
Model.H = H;
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
warning('off', 'MATLAB:nearlySingularMatrix')
[~, ~, ~, info, Model, DynareResults] = ...
dynare_resolve(Model, DynareOptions, DynareResults, 'restrict');
warning('on', 'MATLAB:nearlySingularMatrix')
if info(1)~=0
if nargout==5
ReducedForm = 0;
end
return
end
% Get decision rules and transition equations.
dr = DynareResults.dr;
% Set persistent variables (first call).
if isempty(init_flag)
mf0 = BayesInfo.mf0;
mf1 = BayesInfo.mf1;
restrict_variables_idx = dr.restrict_var_list;
state_variables_idx = restrict_variables_idx(mf0);
number_of_state_variables = length(mf0);
init_flag = true;
end
% Return reduced form model.
if nargout>4
ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
if DynareOptions.order==2
ReducedForm.use_k_order_solver = false;
ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:);
elseif DynareOptions.order>=3
ReducedForm.use_k_order_solver = true;
ReducedForm.dr = dr;
ReducedForm.udr = folded_to_unfolded_dr(dr, Model, DynareOptions);
else
n_states=size(dr.ghx,2);
n_shocks=size(dr.ghu,2);
ReducedForm.use_k_order_solver = false;
ReducedForm.ghxx = zeros(size(restrict_variables_idx,1),n_states^2);
ReducedForm.ghuu = zeros(size(restrict_variables_idx,1),n_shocks^2);
ReducedForm.ghxu = zeros(size(restrict_variables_idx,1),n_states*n_shocks);
ReducedForm.constant = ReducedForm.steadystate;
end
ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
ReducedForm.Q = Q;
ReducedForm.H = H;
ReducedForm.mf0 = mf0;
ReducedForm.mf1 = mf1;
end
% Set initial condition
if setinitialcondition
switch DynareOptions.particle.initialization
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
[A,B] = kalman_transition_matrix(dr,dr.restrict_var_list,dr.restrict_columns,Model.exo_nbr);
StateVectorVariance2 = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
StateVectorVariance = lyapunov_symm(A, B*ReducedForm.Q*B', DynareOptions.lyapunov_fixed_point_tol, ...
DynareOptions.qz_criterium, DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
StateVectorVariance = StateVectorVariance(mf0,mf0);
case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
old_DynareOptionsperiods = DynareOptions.periods;
DynareOptions.periods = 5000;
old_DynareOptionspruning = DynareOptions.pruning;
DynareOptions.pruning = DynareOptions.particle.pruning;
y_ = simult(dr.ys, dr, Model, DynareOptions, DynareResults);
y_ = y_(dr.order_var(state_variables_idx),2001:DynareOptions.periods);
StateVectorVariance = cov(y_');
DynareOptions.periods = old_DynareOptionsperiods;
DynareOptions.pruning = old_DynareOptionspruning;
clear('old_DynareOptionsperiods','y_');
case 3% Initial state vector covariance is a diagonal matrix.
StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables);
otherwise
error('Unknown initialization option!')
end
ReducedForm.StateVectorMean = StateVectorMean;
ReducedForm.StateVectorVariance = StateVectorVariance;
end

36
matlab/nonlinear-filters/spherical_radial_sigma_points.m Normal file
View File

@ -0,0 +1,36 @@
function [nodes,weights] = spherical_radial_sigma_points(n)
%
% Computes nodes and weigths from a third-degree spherical-radial cubature
% rule.
% INPUTS
% n [integer] scalar, number of variables.
%
% OUTPUTS
% nodes [double] nodes of the cubature
% weigths [double] associated weigths
%
% REFERENCES
%
% Arasaratnam & Haykin 2008,2009.
%
% NOTES
%
% Copyright © 2009-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
nodes = (sqrt(n)*([eye(n) -eye(n)]))' ;
weights = (1/(2*n)) ;

237
matlab/nonlinear-filters/traditional_resampling.m Normal file
View File

@ -0,0 +1,237 @@
function return_resample = traditional_resampling(particles,weights,noise)
% Resamples particles.
%@info:
%! @deftypefn {Function File} {@var{indx} =} traditional_resampling (@var{weights},@var{noise})
%! @anchor{particle/traditional_resampling}
%! @sp 1
%! Resamples particles (Resampling à la Kitagawa or stratified resampling).
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item weights
%! n*1 vector of doubles, particles' weights.
%! @item noise
%! n*1 vector of doubles sampled from a [0,1] uniform distribution (stratified resampling) or scalar double
%! sampled from a [0,1] uniform distribution (Kitagawa resampling).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item indx
%! n*1 vector of intergers, indices.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{particle/resample}
%! @sp 2
%! @strong{This function calls:}
%! @sp 2
%! @end deftypefn
%@eod:
% Copyright © 2011-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% 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

82
matlab/nonlinear-filters/univariate_smooth_resampling.m Normal file
View File

@ -0,0 +1,82 @@
function new_particles = univariate_smooth_resampling(weights,particles,number_of_new_particles)
% Smooth Resampling of the particles.
%@info:
%! @deftypefn {Function File} {@var{new_particles} =} univariate_smooth_resampling (@var{weights}, @var{number_of_new_particles})
%! @anchor{particle/univariate_smooth_resampling}
%! @sp 1
%! Smooth Resampling of the particles (univariate version).
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item weights
%! n*1 vector of doubles, particles' weights.
%! @item particles
%! n*1 vector of doubles, particles.
%! @item number_of_new_particles
%! Integer scalar.
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item indx
%! number_of_new_particles*1 vector of doubles, new particles.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{particle/sequantial_importance_particle_filter}
%! @sp 2
%! @strong{This function calls:}
%! @sp 2
%! @end deftypefn
%@eod:
% Copyright © 2012-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% 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

40
matlab/nonlinear-filters/unscented_sigma_points.m Normal file
View File

@ -0,0 +1,40 @@
function [nodes,W_m,W_c] = unscented_sigma_points(n,ParticleOptions)
%
% Computes nodes and weigths for a scaled unscented transform cubature
% INPUTS
% n [integer] scalar, number of variables.
%
% OUTPUTS
% nodes [double] nodes of the cubature
% weigths [double] associated weigths
%
% REFERENCES
%
%
%
% NOTES
%
% Copyright © 2009-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
lambda = (ParticleOptions.unscented.alpha^2)*(n+ParticleOptions.unscented.kappa) - n ;
nodes = [ zeros(n,1) ( sqrt(n+lambda).*([ eye(n) -eye(n)]) ) ]' ;
W_m = lambda/(n+lambda) ;
W_c = W_m + (1-ParticleOptions.unscented.alpha^2+ParticleOptions.unscented.beta) ;
temp = ones(2*n,1)/(2*(n+lambda)) ;
W_m = [W_m ; temp] ;
W_c = [W_c ; temp] ;

1
matlab/particles

@ -1 +0,0 @@
Subproject commit 4a679b26eef0463169705cde460079feddc14819