Compare commits
113 Commits
master
...
rm-particl
Author | SHA1 | Date |
---|---|---|
Stéphane Adjemian (Ryûk) | 8af6fbd829 | |
Stéphane Adjemian (Ryûk) | 683fa6e028 | |
Stéphane Adjemian (Ryûk) | bb8724afa9 | |
Sébastien Villemot | 4a679b26ee | |
Normann Rion | 558990e447 | |
Stéphane Adjemian (Ryûk) | bb16df7f2e | |
Sébastien Villemot | 7926e75aa5 | |
Normann Rion | 80d48b7745 | |
Normann Rion | 25401abeb8 | |
Normann Rion | 363056a3d1 | |
Sébastien Villemot | d71f8be6e7 | |
Frédéric Karamé | 4942488287 | |
Sébastien Villemot | d20be2c271 | |
Sébastien Villemot | 2e18ca8ff5 | |
Johannes Pfeifer | e5a7688433 | |
Sébastien Villemot | 383df9c85c | |
Johannes Pfeifer | 13fe810453 | |
Sébastien Villemot | 1451bf64a2 | |
Johannes Pfeifer | 648c825040 | |
Stéphane Adjemian (Ryuk) | a5c2c4af4b | |
Stéphane Adjemian (Ryuk) | e31e4a0626 | |
Stéphane Adjemian (Ryuk) | ec50899f7a | |
Stéphane Adjemian (Ryuk) | 6f791870b4 | |
Sébastien Villemot | 99aa1fc04b | |
Stéphane Adjemian | 41320041c3 | |
Johannes Pfeifer | e685fa396f | |
Johannes Pfeifer | 60f71f2720 | |
Johannes Pfeifer | f4003536db | |
Sébastien Villemot | 1c3153741f | |
Johannes Pfeifer | 3dca797ba8 | |
Sébastien Villemot | 6258cca240 | |
Johannes Pfeifer | 178eaceb9c | |
Johannes Pfeifer | c440779a11 | |
Stéphane Adjemian (Charybdis) | cef2015162 | |
Stéphane Adjemian (Charybdis) | 2ed534e52e | |
Stéphane Adjemian (Charybdis) | 313003b145 | |
Stéphane Adjemian (Charybdis) | 472d755d98 | |
Stéphane Adjemian (Charybdis) | 89a941605f | |
Stéphane Adjemian (Charybdis) | 8eaccada38 | |
Stéphane Adjemian (Charybdis) | 13d79cb438 | |
Stéphane Adjemian (Charybdis) | a0fb5c7348 | |
Frédéric Karamé | 22be3797a8 | |
Sébastien Villemot | 7c85187319 | |
Frédéric Karamé | 66867386ef | |
Frédéric Karamé | fb7ec278d6 | |
Frédéric Karamé | 35c3a9fcd9 | |
Frédéric Karamé | ef22c716d3 | |
Frédéric Karamé | 43615ce4f6 | |
Frédéric Karamé | 1f08164b05 | |
Frédéric Karamé | 28f7c7621b | |
Frédéric Karamé | ee6eaa8449 | |
Frédéric Karamé | 427e88e6b4 | |
Frédéric Karamé | 857168ddf8 | |
Stéphane Adjemian (Charybdis) | 84d213ea3e | |
Stéphane Adjemian (Charybdis) | 4921000f0a | |
Houtan Bastani | 823d947484 | |
Houtan Bastani | 20978ce117 | |
Stéphane Adjemian | 79b42ec947 | |
Johannes Pfeifer | ddc5b5b665 | |
Frédéric Karamé | 285d5c1711 | |
Frédéric Karamé | 6d1036b738 | |
Houtan Bastani | 0fbab3d401 | |
Frédéric Karamé | 81c606abb2 | |
Frédéric Karamé | 30540e3940 | |
Frédéric Karamé | a9800bef4d | |
Stéphane Adjemian (Charybdis) | af66ea9c28 | |
Stéphane Adjemian (Charybdis) | 8c65631645 | |
Frédéric Karamé | a73e4d94cb | |
Frédéric Karamé | a6a3e1b109 | |
Frédéric Karamé | c8c14f426b | |
Frédéric Karamé | c50392de4b | |
Frédéric Karamé | 7dbc2accae | |
Frédéric Karamé | 8d401a22f5 | |
Frédéric Karamé | 5922f88182 | |
Frédéric Karamé | f37094b20a | |
Frédéric Karamé | 18059c01a6 | |
Frédéric Karamé | 5b5e88b50d | |
Frédéric Karamé | 93039cd92e | |
Frédéric Karamé | 04ad104bfb | |
Frédéric Karamé | cdc7f6ddf9 | |
Stéphane Adjemian | 7235de438c | |
Johannes Pfeifer | 5220ec2a1a | |
Stéphane Adjemian (Hermes) | 12a5e664f4 | |
Frédéric Karamé | 3e22ed1bc2 | |
Frédéric Karamé | c89801b2df | |
Frédéric Karamé | ef04269486 | |
Frédéric Karamé | 5cd887b048 | |
Frédéric Karamé | ec4f2da2fd | |
Frédéric Karamé | 6e7090f1f5 | |
Frédéric Karamé | 2995562e68 | |
Frédéric Karamé | 5880db7537 | |
Frédéric Karamé | d974c77df9 | |
Frédéric Karamé | aad6c297fe | |
Frédéric Karamé | 473b2b8bb1 | |
Frédéric Karamé | 3b04e52667 | |
Frédéric Karamé | 25abaf2f14 | |
Stéphane Adjemian | 7a418a60b8 | |
Stéphane Adjemian (Charybdis) | f512e39453 | |
Stéphane Adjemian (Charybdis) | f66b6b6e8c | |
Stéphane Adjemian (Charybdis) | 7ad84f3f2e | |
Stéphane Adjemian (Charybdis) | 80af764efc | |
Stéphane Adjemian (Charybdis) | f4a1f065ef | |
Stéphane Adjemian (Charybdis) | 21f97c8a98 | |
Stéphane Adjemian (Charybdis) | f82afdf8aa | |
Stéphane Adjemian (Charybdis) | f187c3d414 | |
Frédéric Karamé | e55f3fd272 | |
Frédéric Karamé | bf9e0490ae | |
Stéphane Adjemian (Telemachus) | ddd006ffc3 | |
Stéphane Adjemian (Telemachus) | 0d11f56e06 | |
Stéphane Adjemian (Telemachus) | f7b5d53f5b | |
Stéphane Adjemian (Telemachus) | b577de9b87 | |
Stéphane Adjemian (Telemachus) | bd951bdcce | |
Stéphane Adjemian (Telemachus) | 01c35cc1dc |
|
@ -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
|
||||
|
|
|
@ -53,7 +53,7 @@ p = {'/distributions/' ; ...
|
|||
'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ...
|
||||
'/../contrib/jsonlab/' ; ...
|
||||
'/parallel/' ; ...
|
||||
'/particles/src' ; ...
|
||||
'/nonlinear-filters' ; ...
|
||||
'/gsa/' ; ...
|
||||
'/ep/' ; ...
|
||||
'/backward/' ; ...
|
||||
|
|
|
@ -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)' ;
|
|
@ -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));
|
|
@ -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
|
|
@ -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));
|
|
@ -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
|
|
@ -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;
|
|
@ -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));
|
|
@ -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
|
|
@ -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;
|
|
@ -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));
|
|
@ -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;
|
|
@ -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);
|
|
@ -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
|
|
@ -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') ;
|
|
@ -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
|
|
@ -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);
|
|
@ -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));
|
|
@ -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
|
|
@ -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
|
|
@ -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 ;
|
|
@ -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
|
|
@ -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
|
|
@ -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
|
|
@ -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));
|
|
@ -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
|
|
@ -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)) ;
|
|
@ -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
|
|
@ -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
|
|
@ -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 +0,0 @@
|
|||
Subproject commit 4a679b26eef0463169705cde460079feddc14819
|
Loading…
Reference in New Issue