Skip to content

GitLab

Commit 9f1ad556 authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

New nonlinear filters routines and bug fixes. 

Merge remote-tracking branch 'refs/remotes/fred/master'
parents 93e6b41d 6139af89
Hide whitespace changes
Inline Side-by-side
matlab/global_initialization.m
View file @ 9f1ad556
......@@ -200,7 +200,7 @@ options_.bnlms = bnlms;
particle.status = 0;
% How do we initialize the states?
particle.initialization = 1;
particle.initial_state_prior_std = .0001;
particle.initial_state_prior_std = .1;
% Set the default order of approximation of the model (perturbation).
particle.perturbation = 2;
% Set the default number of particles.
......@@ -211,12 +211,14 @@ particle.smolyak_accuracy = 3;
particle.pruning = 0;
% Set default algorithm
particle.algorithm = 'sequential_importance_particle_filter';
% Set the Gaussian approximation method
% Set the Gaussian approximation method for Importance sampling
particle.IS_approximation_method = 'unscented';
% Set the Gaussian approximation method for particles distributions
particle.approximation_method = 'unscented';
% Set unscented parameters alpha, beta and kappa for gaussian approximation
particle.unscented.alpha = 1;
particle.unscented.alpha = .01;
particle.unscented.beta = 2;
particle.unscented.kappa = 1;
particle.unscented.kappa = 0;
% Configuration of resampling in case of particles
particle.resampling.status = 'systematic'; % 'none', 'generic', 'smoothed'
particle.resampling.neff_threshold = .5;
......@@ -224,9 +226,9 @@ particle.resampling.neff_threshold = .5;
particle.resampling.method1 = 'traditional' ;
particle.resampling.method2 = 'kitagawa';
% Number of partitions for the smoothed resampling method
DynareOptions.particle.resampling.number_of_partitions = 200;
particle.resampling.number_of_partitions = 200;
% Configuration of the mixture filters
particle.mixture_method = 'particles' ;
%particle.mixture_method = 'particles' ;
% Size of the different mixtures
particle.mixture_state_variables = 5 ;
particle.mixture_structural_shocks = 1 ;
......
matlab/metropolis_hastings_initialization.m
View file @ 9f1ad556
......@@ -143,8 +143,8 @@ if ~options_.load_mh_file && ~options_.mh_recover
disp(' ')
else% Case 2: one chain (we start from the posterior mode)
fprintf(fidlog,[' Initial values of the parameters:\n']);
candidate = transpose(xparam1);
if all(candidate' > mh_bounds(:,1)) && all(candidate' < mh_bounds(:,2))
candidate = transpose(xparam1(:));%
if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2))
ix2 = candidate;
ilogpo2 = - feval(TargetFun,ix2',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
disp('MH: Initialization at the posterior mode.')
......
matlab/non_linear_dsge_likelihood.m
View file @ 9f1ad556
......@@ -142,7 +142,7 @@ nvobs = DynareDataset.info.nvobs;
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if (DynareOptions.mode_compute~=1) && any(xparam1<BayesInfo.lb)
k = find(xparam1 < BayesInfo.lb);
k = find(xparam1(:) < BayesInfo.lb);
fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
exit_flag = 0;
info = 41;
......@@ -151,7 +151,7 @@ end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if (DynareOptions.mode_compute~=1) && any(xparam1>BayesInfo.ub)
k = find(xparam1>BayesInfo.ub);
k = find(xparam1(:)>BayesInfo.ub);
fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
exit_flag = 0;
info = 42;
......@@ -345,7 +345,11 @@ ReducedForm.StateVectorVariance = StateVectorVariance;
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
DynareOptions.warning_for_steadystate = 0;
seed_norm_old = randn('state') ;
seed_uniform_old = rand('state') ;
LIK = feval(DynareOptions.particle.algorithm,ReducedForm,Y,[],DynareOptions);
randn('state',seed_norm_old);
rand('state',seed_uniform_old);
if imag(LIK)
likelihood = objective_function_penalty_base;
exit_flag = 0;
......@@ -359,5 +363,5 @@ DynareOptions.warning_for_steadystate = 1;
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
lnprior = priordens(xparam1(:),BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
fval = (likelihood-lnprior);
\ No newline at end of file
matlab/particle/auxiliary_particle_filter.m
View file @ 9f1ad556
......@@ -20,7 +20,7 @@ function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,DynareOptions
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2011-2012 Dynare Team
% Copyright (C) 2011, 2012 Dynare Team
%
% This file is part of Dynare.
%
......@@ -44,6 +44,9 @@ if isempty(start)
start = 1;
end
% Set flag for prunning
pruning = DynareOptions.particle.pruning;
% Get steady state and mean.
steadystate = ReducedForm.steadystate;
constant = ReducedForm.constant;
......@@ -80,6 +83,10 @@ 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');
......@@ -90,30 +97,46 @@ lik = NaN(sample_size,1);
LIK = NaN;
% Initialization of the weights across particles.
weights = ones(1,number_of_particles);
weights = ones(1,number_of_particles)/number_of_particles ;
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
StateVectors_ = StateVectors;
end
for t=1:sample_size
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
PredictedObservedMean = mean(tmp(mf1,:),2);
if pruning
yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state);
[tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
else
tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
end
PredictedObservedMean = weights*(tmp(mf1,:)');
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles+H;
dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean');
PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' +H;
wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
tau_tilde = weights.*wtilde ;
sum_tau_tilde = sum(tau_tilde) ;
lik(t) = log(sum_tau_tilde) ;
%var_wtilde = wtilde-sum_tau_tilde ;
%var_wtilde = var_wtilde'*var_wtilde/(number_of_particles-1) ;
lik(t) = log(sum_tau_tilde) ; %+ .5*var_wtilde/(number_of_particles*(sum_tau_tilde*sum_tau_tilde)) ;
tau_tilde = tau_tilde/sum_tau_tilde;
indx_resmpl = resample(tau_tilde,DynareOptions.particle.resampling.method1,DynareOptions.particle.resampling.method2);
yhat = yhat(:,indx_resmpl);
wtilde = wtilde(indx_resmpl);
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
if pruning
yhat_ = yhat_(:,indx_resmpl);
[tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
StateVectors_ = tmp_(mf0,:);
else
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
end
StateVectors = tmp(mf0,:);
PredictedObservedMean = mean(tmp(mf1,:),2);
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles+H;
PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles + H;
lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1)));
wtilde = lnw./wtilde;
weights = wtilde/sum(wtilde);
......
matlab/particle/fit_gaussian_mixture.m 0 → 100644
View file @ 9f1ad556
function [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(X,StateMu,StateSqrtP,StateWeights,crit,niters,check)
[dim,Ndata] = size(X);
M = size(StateMu,2) ;
if check % Ensure that covariances don't collapse
MIN_COVAR_SQRT = sqrt(eps);
init_covars = StateSqrtP;
end
eold = -Inf;
for n=1:niters
% Calculate posteriors based on old parameters
[prior,likelihood,marginal,posterior] = probability(StateMu,StateSqrtP,StateWeights,X);
e = sum(log(marginal));
if (n > 1 && abs((e - eold)/eold) < crit)
return;
else
eold = e;
end
new_pr = (sum(posterior,2))';
StateWeights = new_pr/Ndata;
StateMu = bsxfun(@rdivide,(posterior*X')',new_pr);
for j=1:M
diffs = bsxfun(@minus,X,StateMu(:,j));
tpost = (1/sqrt(new_pr(j)))*sqrt(posterior(j,:));
diffs = bsxfun(@times,diffs,tpost);
[foo,tcov] = qr2(diffs') ;
StateSqrtP(:,:,j) = tcov';
if check
if min(abs(diag(StateSqrtP(:,:,j)))) < MIN_COVAR_SQRT
StateSqrtP(:,:,j) = init_covars(:,:,j);
end
end
end
end
matlab/particle/gaussian_densities.m 0 → 100644
View file @ 9f1ad556
function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,DynareOptions)
%
% Elements to calculate the importance sampling ratio
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
% reduced_form_model.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% smolyak_accuracy [integer] scalar.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2010 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://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,DynareOptions) ;
eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ;
sqr_det = sqrt(det(Pyy)) ;
foo = (eta_t_i/Pyy).*eta_t_i ;
likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ;
IncrementalWeights = likelihood.*prior./proposal ;
matlab/particle/gaussian_filter.m 0 → 100644
View file @ 9f1ad556
function [LIK,lik] = gaussian_filter(ReducedForm,Y,start,DynareOptions)
% Evaluates the likelihood of a non-linear model approximating the
% predictive (prior) and filtered (posterior) densities for state variables
% by gaussian distributions.
% Gaussian approximation is done by:
% - a Kronrod-Paterson gaussian quadrature with a limited number of nodes.
% Multidimensional quadrature is obtained by the Smolyak operator (ref: Winschel & Kratzig, 2010).
% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2008,2009).
% - a scaled unscented transform cubature (ref: )
% - Monte-Carlo draws from a multivariate gaussian distribution.
% First and second moments of prior and posterior state densities are computed
% from the resulting nodes/particles and allows to generate new distributions at the
% following observation.
% => The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles
% filters since it treats a lesser number of particles and there is no need
% of resampling.
% However, estimations may reveal biaised if the model is truly non-gaussian
% since predictive and filtered densities are unimodal.
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
% reduced_form_model.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% smolyak_accuracy [integer] scalar.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2010 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent init_flag mf0 mf1
persistent nodes2 weights2 weights_c2 number_of_particles
persistent sample_size number_of_state_variables number_of_observed_variables
verbose = 1;
% Set default
if isempty(start)
start = 1;
end
% Set persistent variables.
if isempty(init_flag)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_particles = DynareOptions.particle.number_of_particles;
init_flag = 1;
end
% compute gaussian quadrature nodes and weights on states and shocks
if isempty(nodes2) && strcmpi(DynareOptions.particle.approximation_method,'quadrature')
[nodes2,weights2] = nwspgr('GQN',number_of_state_variables,DynareOptions.particle.smolyak_accuracy) ;
weights_c2 = weights2 ;
end
if isempty(nodes2) && strcmpi(DynareOptions.particle.approximation_method,'cubature')
[nodes2,weights2] = spherical_radial_sigma_points(number_of_state_variables) ;
weights_c2 = weights2 ;
end
if isempty(nodes2) && strcmpi(DynareOptions.particle.approximation_method,'unscented')
[nodes2,weights2,weights_c2] = unscented_sigma_points(number_of_state_variables,DynareOptions) ;
end
if isempty(nodes2) && strcmpi(DynareOptions.particle.approximation_method,'monte-carlo')
set_dynare_seed('default');
end
% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
if isempty(H)
H = 0;
H_lower_triangular_cholesky = 0;
else
H_lower_triangular_cholesky = reduced_rank_cholesky(H)';
end
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)';
% Initialization of the likelihood.
const_lik = (2*pi)^(number_of_observed_variables/2) ;
lik = NaN(sample_size,1);
LIK = NaN;
SampleWeights = 1/number_of_particles ;
ks = 0 ;
%Estimate = zeros(number_of_state_variables,sample_size) ;
%V_Estimate = zeros(number_of_state_variables,number_of_state_variables,sample_size) ;
for t=1:sample_size
% build the proposal
[PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = ...
gaussian_filter_bank(ReducedForm,Y(:,t),StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions) ;
%Estimate(:,t) = PredictedStateMean ;
%V_Estimate(:,:,t) = PredictedStateVarianceSquareRoot ;
if strcmpi(DynareOptions.particle.approximation_method,'quadrature') || ... % sparse grids approximations
strcmpi(DynareOptions.particle.approximation_method,'cubature') || ...
strcmpi(DynareOptions.particle.approximation_method,'unscented')
StateParticles = bsxfun(@plus,StateVectorMean,StateVectorVarianceSquareRoot*nodes2') ;
IncrementalWeights = ...
gaussian_densities(Y(:,t),StateVectorMean,...
StateVectorVarianceSquareRoot,PredictedStateMean,...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,...
weights2,weights_c2,ReducedForm,DynareOptions) ;
SampleWeights = weights2.*IncrementalWeights ;
SumSampleWeights = sum(SampleWeights) ;
lik(t) = log(SumSampleWeights) ;
SampleWeights = SampleWeights./SumSampleWeights ;
else % Monte-Carlo draws
StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean) ;
IncrementalWeights = ...
gaussian_densities(Y(:,t),StateVectorMean,...
StateVectorVarianceSquareRoot,PredictedStateMean,...
PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,...
1/number_of_particles,1/number_of_particles,ReducedForm,DynareOptions) ;
SampleWeights = SampleWeights.*IncrementalWeights ;
SumSampleWeights = sum(SampleWeights) ;
%VarSampleWeights = IncrementalWeights-SumSampleWeights ;
%VarSampleWeights = VarSampleWeights*VarSampleWeights'/(number_of_particles-1) ;
lik(t) = log(SumSampleWeights) ; %+ .5*VarSampleWeights/(number_of_particles*(SumSampleWeights*SumSampleWeights)) ;
SampleWeights = SampleWeights./SumSampleWeights ;
Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
if (Neff<.5*sample_size && strcmpi(DynareOptions.particle.resampling.status,'generic')) || ...
strcmpi(DynareOptions.particle.resampling.status,'systematic')
ks = ks + 1 ;
StateParticles = StateParticles(:,resample(SampleWeights',DynareOptions.particle.resampling.method1,DynareOptions.particle.resampling.method2)) ;
StateVectorMean = mean(StateParticles,2) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )';
SampleWeights = 1/number_of_particles ;
elseif strcmp(DynareOptions.particle.resampling.status,'smoothed')
StateParticles = multivariate_smooth_resampling(SampleWeights,StateParticles',number_of_particles,DynareOptions.particle.resampling.number_of_partitions)';
StateVectorMean = mean(StateParticles,2) ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )';
SampleWeights = 1/number_of_particles ;
elseif strcmpi(DynareOptions.particle.resampling.status,'none')
StateVectorMean = (sampleWeights*StateParticles)' ;
temp = sqrt(SampleWeights').*StateParticles ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp'*temp - StateVectorMean*(StateVectorMean') )';
end
end
end
LIK = -sum(lik(start:end));
matlab/particle/gaussian_filter_bank.m 0 → 100644
View file @ 9f1ad556
function [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = gaussian_filter_bank(ReducedForm,obs,StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions)
%
% Computes the proposal with a gaussian approximation for importance
% sampling
% This proposal is a gaussian distribution calculated à la Kalman
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
% reduced_form_model.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2010 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent init_flag2 mf0 mf1
persistent number_of_state_variables number_of_observed_variables
persistent number_of_structural_innovations
% Set local state space model (first-order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
% Set local state space model (second-order approximation).
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ...
any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ...
any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4))
ghx
ghu
ghxx
ghuu
ghxu
end
constant = ReducedForm.constant;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set persistent variables.
if isempty(init_flag2)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
init_flag2 = 1;
end
if strcmpi(DynareOptions.particle.IS_approximation_method,'cubature') || strcmpi(DynareOptions.particle.IS_approximation_method,'monte-carlo')
[nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations) ;
weights_c = weights ;
end
if strcmpi(DynareOptions.particle.IS_approximation_method,'quadrature')
[nodes,weights] = nwspgr('GQN',number_of_state_variables+number_of_structural_innovations,DynareOptions.particle.smolyak_accuracy) ;
weights_c = weights ;
end
if strcmpi(DynareOptions.particle.IS_approximation_method,'unscented')
[nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,DynareOptions) ;
end
xbar = [StateVectorMean ; zeros(number_of_structural_innovations,1) ] ;
sqr_Px = [ [ StateVectorVarianceSquareRoot zeros(number_of_state_variables,number_of_structural_innovations) ] ;
[ zeros(number_of_structural_innovations,number_of_state_variables) Q_lower_triangular_cholesky ] ] ;
sigma_points = bsxfun(@plus,xbar,sqr_Px*(nodes')) ;
StateVectors = sigma_points(1:number_of_state_variables,:) ;
epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:) ;
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
PredictedStateMean = tmp(mf0,:)*weights ;
PredictedObservedMean = tmp(mf1,:)*weights;
if strcmpi(DynareOptions.particle.IS_approximation_method,'cubature') || strcmpi(DynareOptions.particle.IS_approximation_method,'monte-carlo')
PredictedStateMean = sum(PredictedStateMean,2) ;
PredictedObservedMean = sum(PredictedObservedMean,2) ;
dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean)'.*sqrt(weights) ;
dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean)'.*sqrt(weights);
PredictedStateVarianceSquareRoot = chol(dState'*dState)';
big_mat = [dObserved dState ; [H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)] ] ;
[mat1,mat] = qr2(big_mat) ;
mat = mat' ;
clear('mat1');
PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables) ;
CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables) ;
StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables)) ;
PredictionError = obs - PredictedObservedMean ;
StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*PredictionError ;
end
if strcmpi(DynareOptions.particle.IS_approximation_method,'quadrature') || strcmpi(DynareOptions.particle.IS_approximation_method,'unscented')
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
matlab/particle/gaussian_mixture_densities.m 0 → 100644
View file @ 9f1ad556
function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
StateMuPost,StateSqrtPPost,StateWeightsPost,...
StateParticles,H,normconst,weigths1,weigths2,ReducedForm,DynareOptions)
%
% Elements to calculate the importance sampling ratio
%
% INPUTS
% reduced_form_model [structure] Matlab's structure describing the reduced form model.
% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
% reduced_form_model.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% smolyak_accuracy [integer] scalar.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES