Commit 11bb96f1 authored by stepan's avatar stepan
Browse files

* Added subfolder for particle filters.

* Added a routine for gaussian particle filter. The likelihood is evaluated
by assuming that the state vector is normally distributed.

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3109 ac1d8469-bf42-47a9-8791-bf33cf982152
parent b9980474
function [LIK,lik] = gaussian_particle_filter(reduced_form_model,Y,start,mf,number_of_particles,grid_size)
% hparam,y,nbchocetat,nbchocmesure,smol_prec,nb_part,g,m,choix
% Evaluates the likelihood of a non linear model assuming that the particles are normally distributed.
%
% INPUTS
% reduced_form_model [structure] Matlab's structure desvcribing 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'.
% mf [integer] pp*1 vector of indices.
% number_of_particles [integer] scalar.
% grid_size [integer] scalar, size of the smoliak grid.
%
% 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 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/>.
global M_
dr = reduced_form_model.state.dr;
Q = reduced_form_model.state.Q;
H = reduced_form_model.measurement.H;
smpl = size(Y,2); % Sample size.
mm = size(dr.ghx,2); % Number of state variables.
pp = size(Y,1); % Maximum number of observed variables.
qq = length(Q);
lik = NaN(smpl,1);
LIK = 0 ;
if nargin==4
number_of_particles = 10 ;
method = 'mc' ;
end
if nargin==5
method = 'mc' ;
if isempty(number_of_particles)
number_of_particles = 5000 ;
end
end
if nargin==6
method = 'quadra' ;
if isempty(grid_size)
grid_size = 4 ;
end
end
if isempty(start)
start = 1;
end
k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*M_.endo_nbr;
moy_s_t_1_t_1 = dr.ys ;
tot = M_.endo_nbr ;
Pss_t_1_t_1 = lyapunov_symm(dr.ghx(k2,:),dr.ghu(k2,:)*Q*dr.ghu(k2,:)',1e-12,1e-12);
[V,D] = eig(Pss_t_1_t_1) ;
sqrtP = V*sqrt(D)*V';
dim_var = qq + mm ;
switch method
case 'mc'
seed = [ 362436069 ; 521288629 ];
randn('state',seed);
nodes = randn(number_of_particles,dim_var) ;
weights = 1/number_of_particles ;
nb_grid = number_of_particles ;
case 'quadra'
%{nodes,weights} = integration_smolyak(dim_var,grid_size) ;
nb_grid = size(nodes,1) ;
otherwise
error('undefined method') ;
end
var_idx = 1:mm;
inn_idx = mm+1:dim_var;
chol_Q = chol(Q);
iorder = 2;
s_t_1_t_1 = zeros(nb_grid,tot) ;
iorder=2;
for t=1:smpl
s_t_1_t_1(:,k2) = (repmat(moy_s_t_1_t_1(k2),1,nb_grid) + sqrtP*nodes(:,var_idx)')' ;
e_t = (chol_Q*nodes(:,inn_idx)')' ;
s_t_t_1 = 0 ;
y_t_t_1 = 0 ;
Pss_t_t_1 = 0 ;
Pyy_t_t_1 = 0 ;
Psy_t_t_1 = 0 ;
for i=1:nb_grid
tout = simult_(s_t_1_t_1(i,:)',dr,e_t(i,:),iorder) ;
s_t_t_1 = s_t_t_1 + tout(k2,2) ;
y_t_t_1 = y_t_t_1 + tout(mf,2) ;
Pyy_t_t_1 = Pyy_t_t_1 + tout(mf,2)*tout(mf,2)' ;
Psy_t_t_1 = Psy_t_t_1 + tout(k2,2)*tout(mf,2)' ;
Pss_t_t_1 = Pss_t_t_1 + tout(k2,2)*tout(k2,2)' ;
end
s_t_t_1 = s_t_t_1/nb_grid ;
y_t_t_1 = y_t_t_1/nb_grid ;
Pyy_t_t_1 = Pyy_t_t_1 + H - y_t_t_1*(y_t_t_1') ;
Psy_t_t_1 = Psy_t_t_1 - s_t_t_1*(y_t_t_1');
Pss_t_t_1 = Pss_t_t_1 - s_t_t_1*(s_t_t_1');
inv_Pyy_t_t_1 = inv(Pyy_t_t_1) ;
eta_t = Y(:,t) - y_t_t_1 ;
Kt = Psy_t_t_1*inv_Pyy_t_t_1 ;
moy_s_t_1_t_1(k2) = s_t_t_1 + Kt*eta_t ;
Pss_t_t_1 = Pss_t_t_1 - Kt*Pyy_t_t_1 *(Kt') ;
[V,D] = eig(Pss_t_1_t_1);
sqrtP = V*sqrt(D)*V';
lik(t) = -0.5*log(2*pi)*pp - 0.5*log(det(Pyy_t_t_1)) - 0.5*eta_t'*inv_Pyy_t_t_1*eta_t ;
if abs(imag(lik(t)))<1e-12
lik(t) = real(lik(t));
end
end
LIK = -sum(lik(start:end));
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment