Skip to content
GitLab
Commit 2063c094 authored by stepan's avatar stepan
Browse files

* Removed simult_ and added local_state_space_iteration_2 (a specialized version of simult_). 

* Renamed gaussian_particle_filter.m. Changed the names of the variables in this file. The 
gaussian nodes are not stored anymore in a big matrix.
* Added a function to compute the cholesky decomposition of a symetric semidefinite positive
matrix.  



git-svn-id: https://www.dynare.org/svn/dynare/trunk@3138 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 26185780
Hide whitespace changes
Inline Side-by-side
matlab/particle/DsgeLikelihood.m
View file @ 2063c094
......@@ -255,7 +255,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
rfm.measurement.H = H;
number_of_particles = 10;
LIK = gaussian_particle_filter(rfm,Y,start,mf,number_of_particles);
LIK = gaussian_particle_filter(rfm,Y);
% ------------------------------------------------------------------------------
......
matlab/particle/gaussian_particle_filter.m deleted 100755 → 0
View file @ 26185780
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));
matlab/particle/local_state_space_iteration_2.m 0 → 100644
View file @ 2063c094
function y = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,half_ghxx,half_ghuu,ghxu)
% Given an initial condition (y) and an innovation (epsilon), this
% routines computes the next value of the endogenous variables if the
% model is approximated by an order two taylor expansion around the
% deterministic steady state.
%
% INPUTS
% yhat [double] n*1 vector, initial condition, where n is the number of state variables.
% epsilon [double] q*1 vector, structural innovations.
% ys [double] m*1 vector, steady state (the variables are ordered as in ghx) where m
% is the number of elements in the union of the states and observed
% variables.
% ghx [double] m*n matrix, is a subset of dr.ghx we only consider the lines corresponding
% to the states and the observed variables.
% ghu [double] m*q matrix, is a subset of dr.ghu
% constant [double] m*1 vector (steady state + second order correction).
% half_ghxx [double] m*n² matrix, subset of .5*dr.ghxx.
% half_ghuu [double] m*q² matrix, subset of .5*dr.ghuu.
% ghxu [double] m*nq matrix, subset of dr.ghxu.
% index [integer]
%
% OUTPUTS
% y [double] stochastic simulations results
%
% SPECIAL REQUIREMENTS
% none
% 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/>.
y = constant + ghx*yhat + ghu*epsilon ...
+ A_times_B_kronecker_C(half_ghxx,yhat) ...
+ A_times_B_kronecker_C(half_ghuu,epsilon) ...
+ A_times_B_kronecker_C(ghxu,yhat,epsilon);
\ No newline at end of file
matlab/particle/monte_carlo_gaussian_particle_filter.m 0 → 100755
View file @ 2063c094
function [LIK,lik] = monte_carlo_gaussian_particle_filter(reduced_form_model,Y,start,number_of_particles)
% 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_ bayestopt_
persistent init_flag
persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1
persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
% Set defaults.
if (nargin<4) || (nargin==4 && isempty(number_of_particles))
number_of_particles = 1000 ;
end
if nargin==2 || isempty(start)
start = 1;
end
dr = reduced_form_model.state.dr;% Decision rules and transition equations.
Q = reduced_form_model.state.Q;% Covariance matrix of the structural innovations.
H = reduced_form_model.measurement.H;% Covariance matrix of the measurement errors.
% Set persistent variables.
if isempty(init_flag)
mf0 = bayestopt_.mf0;
mf1 = bayestopt_.mf1;
restrict_variables_idx = bayestopt_.restrict_var_list;
observed_variables_idx = restrict_variables_idx(mf1);
state_variables_idx = restrict_variables_idx(mf0);
sample_size = size(Y,2);
number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(Q);
init_flag = 1;
end
% Set local state space model (second order approximation).
ghx = dr.ghx(restrict_variables_idx,:);
ghu = dr.ghu(restrict_variables_idx,:);
half_ghxx = .5*dr.ghxx(restrict_variables_idx,:);
half_ghuu = .5*dr.ghuu(restrict_variables_idx,:);
ghxu = dr.ghxu(restrict_variables_idx,:);
steadystate = dr.ys(dr.order_var(restrict_variables_idx));
constant = steadystate + .5*dr.ghs2(restrict_variables_idx);
state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
StateVectorMean = state_variables_steady_state;
StateVectorVariance = lyapunov_symm(ghx(mf0,:),ghu(mf0,:)*Q*ghu(mf0,:)',1e-12,1e-12);
StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
%state_idx = 1:state_variance_rank;
%innovation_idx = 1+state_variance_rank:state_variance_rank+number_of_structural_innovations;
Q_lower_triangular_cholesky = chol(Q)';
% Set seed for randn().
seed = [ 362436069 ; 521288629 ];
randn('state',seed);
const_lik = log(2*pi)*number_of_observed_variables;
lik = NaN(sample_size,1);
for t=1:sample_size
PredictedStateMean = zeros(number_of_state_variables,1);
PredictedObservedMean = zeros(number_of_observed_variables,1);
PredictedStateVariance = zeros(number_of_state_variables);
PredictedObservedVariance = zeros(number_of_observed_variables);
PredictedStateAndObservedCovariance = zeros(number_of_state_variables,number_of_observed_variables);
for i=1:number_of_particles
StateVector = StateVectorMean + StateVectorVarianceSquareRoot*randn(state_variance_rank,1);
yhat = StateVector-state_variables_steady_state;
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,1);
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,half_ghxx,half_ghuu,ghxu);
PredictedStateMean_old = PredictedStateMean;
PredictedObservedMean_old = PredictedObservedMean;
PredictedStateMean = PredictedStateMean + (tmp(mf0)-PredictedStateMean)/i;
PredictedObservedMean = PredictedObservedMean + (tmp(mf1)-PredictedObservedMean)/i;
psm = PredictedStateMean*PredictedStateMean';
pom = PredictedObservedMean*PredictedObservedMean';
pcm = PredictedStateMean*PredictedObservedMean';
PredictedStateVariance = PredictedStateVariance ...
+ ( (tmp(mf0)*tmp(mf0)'-psm-PredictedStateVariance)+(i-1)*(PredictedStateMean_old*PredictedStateMean_old'-psm) )/i;
PredictedObservedVariance = PredictedObservedVariance ...
+ ( (tmp(mf1)*tmp(mf1)'-pom-PredictedObservedVariance)+(i-1)*(PredictedObservedMean_old*PredictedObservedMean_old'-pom) )/i;
PredictedStateAndObservedCovariance = PredictedStateAndObservedCovariance ...
+ ( (tmp(mf0)*tmp(mf1)'-pcm-PredictedStateAndObservedCovariance)+(i-1)*(PredictedStateMean_old*PredictedObservedMean_old'-pcm) )/i;
end
PredictedObservedVariance = PredictedObservedVariance + H;
iPredictedObservedVariance = inv(PredictedObservedVariance);
prediction_error = Y(:,t) - PredictedObservedMean;
filter_gain = PredictedStateAndObservedCovariance*iPredictedObservedVariance;
StateVectorMean = PredictedStateMean + filter_gain*prediction_error;
StateVectorVariance = PredictedStateVariance - filter_gain*PredictedObservedVariance*filter_gain';
StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
lik(t) = -.5*(const_lik + log(det(PredictedObservedVariance)) + prediction_error'*iPredictedObservedVariance*prediction_error);
end
LIK = -sum(lik(start:end));
\ No newline at end of file
matlab/particle/reduced_rank_cholesky.m 0 → 100644
View file @ 2063c094
function T = reduced_rank_cholesky(X)
% Computes the cholesky decomposition of a symetric semidefinite matrix or of a definite positive matrix.
%
% INPUTS:
% X [double] n*n matrix to be factorized.
%
% OUTPUTS
% T [double] q*n matrix such that T'*T = X, where q is the number of positive eigenvalues in X.
%
% NOTES:
% If X is not positive definite, then X has to be a symetric semidefinite matrix.
% The matrix T is upper triangular iff X is positive definite.
% 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/>.
[T,X_is_not_positive_definite] = chol(X);
if X_is_not_positive_definite
n = length(X);
[U,D] = eig(.5*(X+X'));
[tmp,max_elements_indices] = max(abs(U),[],1);
negloc = (U(max_elements_indices+(0:n:(n-1)*n))<0);
U(:,negloc) = -U(:,negloc);
D = diag(D);
tol = eps(max(D)) * length(D)*100;
t = (abs(D) > tol);
D = D(t);
if ~(sum(D<0))
T = diag(sqrt(D)) * U(:,t)';
else
disp('reduced_rank_cholesky:: Input matrix is not semidefinite positive!')
T = NaN;
end
end
\ No newline at end of file
matlab/particle/simult_.m deleted 100644 → 0
View file @ 26185780
function y_=simult_(y0,dr,ex_,iorder)
% function y_=simult_(y0,dr,ex_,iorder)
%
% Simulates the model, given the path of exogenous variables and the
% decision rules.
%
% INPUTS
% y0: starting values
% dr: structure of decisions rules for stochastic simulations
% ex_: matrix of shocks
% iorder=0: first-order approximation
% iorder=1: second-order approximation
%
% OUTPUTS
% y_: stochastic simulations results
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2001-2007 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_ options_ it_
iter = size(ex_,1);
if ~isempty(dr.ghu)
nx = size(dr.ghu,2);
end
y_ = zeros(size(y0,1),iter+M_.maximum_lag);
y_(:,1:M_.maximum_lag) = y0;
k1 = [M_.maximum_lag:-1:1];
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;
k3 = M_.lead_lag_incidence(1:M_.maximum_lag,:)';
k3 = find(k3(:));
k4 = dr.kstate(find(dr.kstate(:,2) < M_.maximum_lag+1),[1 2]);
k4 = k4(:,1)+(M_.maximum_lag+1-k4(:,2))*M_.endo_nbr;
if iorder == 1
if ~isempty(dr.ghu)
for i = M_.maximum_lag+1: iter+M_.maximum_lag
tempx1 = y_(dr.order_var,k1);
tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag);
tempx = tempx2(k2);
y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx+dr.ghu* ...
ex_(i-M_.maximum_lag,:)';
k1 = k1+1;
end
else
for i = M_.maximum_lag+1: iter+M_.maximum_lag
tempx1 = y_(dr.order_var,k1);
tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag);
tempx = tempx2(k2);
y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx;
k1 = k1+1;
end
end
elseif iorder == 2
for i = M_.maximum_lag+1: iter+M_.maximum_lag
tempx1 = y_(dr.order_var,k1);
tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag);
tempx = tempx2(k2);
tempu = ex_(i-M_.maximum_lag,:)';
%tempuu = kron(tempu,tempu);
% tempxx = kron(tempx,tempx);
% tempxu = kron(tempx,tempu);
%y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghs2/2+dr.ghx*tempx+ ...
% dr.ghu*tempu+0.5*(dr.ghxx*tempxx+dr.ghuu*tempuu)+dr.ghxu* ...
% tempxu;
y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghs2/2+dr.ghx*tempx+ ...
dr.ghu*tempu+0.5*(A_times_B_kronecker_C(dr.ghxx,tempx)+A_times_B_kronecker_C(dr.ghuu,tempu))+A_times_B_kronecker_C(dr.ghxu,tempx,tempu);
k1 = k1+1;
end
end
% MJ 08/30/02 corrected bug at order 2
\ No newline at end of file
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