From 8820c63f5e590982f19ba311d40205a5ca7b6e99 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Karam=C3=A9?=
<frederic.karame@univ-lemans.fr>
Date: Tue, 18 Jun 2013 16:15:16 +0200
Subject: [PATCH] Bug fixes.
---
matlab/particle/conditional_particle_filter.m | 4 +-
matlab/particle/gaussian_densities.m | 2 +-
.../particle/multivariate_smooth_resampling.m | 1 -
matlab/particle/online_auxiliary_filter.m | 29 +++++++-------
matlab/particle/resample.m | 1 +
.../sequential_importance_particle_filter.m | 2 +-
.../particle/solve_model_for_online_filter.m | 3 +-
tests/particle/dsge_base2.mod | 40 +++++++++++++------
8 files changed, 48 insertions(+), 34 deletions(-)
diff --git a/matlab/particle/conditional_particle_filter.m b/matlab/particle/conditional_particle_filter.m
index 05b725534e..941e01331e 100644
--- a/matlab/particle/conditional_particle_filter.m
+++ b/matlab/particle/conditional_particle_filter.m
@@ -37,7 +37,7 @@ function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,DynareOptio
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
-% Copyright (C) 2009-2013 Dynare Team
+% Copyright (C) 2009-2010 Dynare Team
%
% This file is part of Dynare.
%
@@ -116,7 +116,7 @@ for t=1:sample_size
if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(SampleWeights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || ...
strcmp(DynareOptions.particle.resampling.status,'systematic')
ks = ks + 1 ;
- StateParticles = resample(StateParticles',SampleWeights,DynareOptions)';
+ StateParticles = resample(StateParticles',SampleWeights',DynareOptions)';
SampleWeights = ones(1,number_of_particles)/number_of_particles ;
end
end
diff --git a/matlab/particle/gaussian_densities.m b/matlab/particle/gaussian_densities.m
index 0c89f0862f..3a770118dd 100644
--- a/matlab/particle/gaussian_densities.m
+++ b/matlab/particle/gaussian_densities.m
@@ -19,7 +19,7 @@ function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sq
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
-% Copyright (C) 2009-2012 Dynare Team
+% Copyright (C) 2009-2010 Dynare Team
%
% This file is part of Dynare.
%
diff --git a/matlab/particle/multivariate_smooth_resampling.m b/matlab/particle/multivariate_smooth_resampling.m
index 8bd67edf70..09ca287b97 100644
--- a/matlab/particle/multivariate_smooth_resampling.m
+++ b/matlab/particle/multivariate_smooth_resampling.m
@@ -59,7 +59,6 @@ function new_particles = multivariate_smooth_resampling(particles,weights)
% stephane DOT adjemian AT univ DASH lemans DOT fr
number_of_particles = length(weights);
-weights = weights' ;
number_of_states = size(particles,2);
[P,D] = eig(particles'*(bsxfun(@times,1/number_of_particles,particles))) ;
D = diag(D) ;
diff --git a/matlab/particle/online_auxiliary_filter.m b/matlab/particle/online_auxiliary_filter.m
index 68a2cd2e98..7089a90389 100644
--- a/matlab/particle/online_auxiliary_filter.m
+++ b/matlab/particle/online_auxiliary_filter.m
@@ -41,7 +41,6 @@ persistent Y init_flag mf0 mf1 bounds number_of_particles number_of_parameters l
persistent start_param sample_size number_of_observed_variables number_of_structural_innovations
% Set seed for randn().
-%set_dynare_seed('mt19937ar',1234) ;
set_dynare_seed('default') ;
pruning = DynareOptions.particle.pruning;
second_resample = 1 ;
@@ -62,11 +61,10 @@ if isempty(init_flag)
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
liu_west_delta = DynareOptions.particle.liu_west_delta ;
- liu_west_chol_sigma_bar = DynareOptions.particle.liu_west_chol_sigma_bar*eye(number_of_parameters) ;
- %start_param = xparam1 ;
- % Conditions initiales
- %liu_west_chol_sigma_bar = bsxfun(@times,eye(number_of_parameters),BayesInfo.p2) ;
- start_param = BayesInfo.p1 ;
+ %liu_west_chol_sigma_bar = DynareOptions.particle.liu_west_chol_sigma_bar*eye(number_of_parameters) ;
+ start_param = xparam1 ;
+ %liu_west_chol_sigma_bar = sqrt(bsxfun(@times,eye(number_of_parameters),BayesInfo.p2)) ;
+ %start_param = BayesInfo.p1 ;
bounds = [BayesInfo.lb BayesInfo.ub] ;
init_flag = 1;
end
@@ -87,11 +85,12 @@ small_a = sqrt(1-h_square) ;
% Initialization of parameter particles
xparam = zeros(number_of_parameters,number_of_particles) ;
-stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/1000 ;
+stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/100 ;
+stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/50 ;
i = 1 ;
while i<=number_of_particles
- candidate = start_param + 10*liu_west_chol_sigma_bar*randn(number_of_parameters,1) ;
- %candidate = start_param + bsxfun(@times,stderr,randn(number_of_parameters,1)) ;
+ %candidate = start_param + .001*liu_west_chol_sigma_bar*randn(number_of_parameters,1) ;
+ candidate = start_param + bsxfun(@times,stderr,randn(number_of_parameters,1)) ;
if all(candidate(:) >= bounds(:,1)) && all(candidate(:) <= bounds(:,2))
xparam(:,i) = candidate(:) ;
i = i+1 ;
@@ -113,7 +112,7 @@ ub95_xparam = zeros(number_of_parameters,sample_size) ;
%% The Online filter
for t=1:sample_size
- disp(t)
+ disp(t)
% Moments of parameters particles distribution
m_bar = xparam*(weights') ;
temp = bsxfun(@minus,xparam,m_bar) ;
@@ -210,7 +209,7 @@ for t=1:sample_size
if (variance_update==1) && (neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size)
variance_update = 0 ;
end
- % final resampling (advised)
+ % final resampling (advised)
if second_resample==1
indx = index_resample(0,weights,DynareOptions);
StateVectors = StateVectors(:,indx) ;
@@ -259,7 +258,8 @@ for t=1:sample_size
end
disp([lb95_xparam(:,t) mean_xparam(:,t) ub95_xparam(:,t)])
end
-xparam = mean_xparam(:,sample_size) ;
+distrib_param = xparam ;
+xparam = mean_xparam(:,sample_size) ;
std_param = std_xparam(:,sample_size) ;
lb_95 = lb95_xparam(:,sample_size) ;
ub_95 = ub95_xparam(:,sample_size) ;
@@ -345,8 +345,8 @@ for plt = 1:nbplt,
TeXNAMES = char(TeXNAMES,texname);
end
end
- optimal_bandwidth = mh_optimal_bandwidth(xparam(kk,:)',number_of_particles,bandwidth,kernel_function);
- [density(:,1),density(:,2)] = kernel_density_estimate(xparam(kk,:)',number_of_grid_points,...
+ optimal_bandwidth = mh_optimal_bandwidth(distrib_param(kk,:)',number_of_particles,bandwidth,kernel_function);
+ [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(kk,:)',number_of_grid_points,...
number_of_particles,optimal_bandwidth,kernel_function);
plot(density(:,1),density(:,2));
hold on
@@ -370,3 +370,4 @@ for plt = 1:nbplt,
fprintf(fidTeX,' \n');
end
end
+
\ No newline at end of file
diff --git a/matlab/particle/resample.m b/matlab/particle/resample.m
index 7ab91083e6..05e5f033fc 100644
--- a/matlab/particle/resample.m
+++ b/matlab/particle/resample.m
@@ -54,6 +54,7 @@ function resampled_particles = resample(particles,weights,DynareOptions)
% AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr
% stephane DOT adjemian AT univ DASH lemans DOT fr
+
switch DynareOptions.particle.resampling.method1
case 'residual'
if strcmpi(DynareOptions.particle.resampling.method2,'kitagawa')
diff --git a/matlab/particle/sequential_importance_particle_filter.m b/matlab/particle/sequential_importance_particle_filter.m
index b8cc6f9261..9db9632d55 100644
--- a/matlab/particle/sequential_importance_particle_filter.m
+++ b/matlab/particle/sequential_importance_particle_filter.m
@@ -161,7 +161,7 @@ for t=1:sample_size
StateVectors = temp(:,1:number_of_state_variables)' ;
StateVectors_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)';
else
- StateVectors = resample(tmp(mf0,:)',weights,DynareOptions)';
+ StateVectors = resample(tmp(mf0,:)',weights',DynareOptions)';
end
weights = ones(1,number_of_particles)/number_of_particles;
elseif strcmp(DynareOptions.particle.resampling.status,'none')
diff --git a/matlab/particle/solve_model_for_online_filter.m b/matlab/particle/solve_model_for_online_filter.m
index 04f7be230f..575d0f16b0 100644
--- a/matlab/particle/solve_model_for_online_filter.m
+++ b/matlab/particle/solve_model_for_online_filter.m
@@ -168,8 +168,7 @@ end
offset = EstimatedParameters.nvx;
if EstimatedParameters.nvn
for i=1:EstimatedParameters.nvn
- k = EstimatedParameters.var_endo(i,1);
- H(k,k) = xparam1(i+offset)*xparam1(i+offset);
+ H(i,i) = xparam1(i+offset)*xparam1(i+offset);
end
offset = offset+EstimatedParameters.nvn;
else
diff --git a/tests/particle/dsge_base2.mod b/tests/particle/dsge_base2.mod
index 990e310cdb..50aacc2c19 100644
--- a/tests/particle/dsge_base2.mod
+++ b/tests/particle/dsge_base2.mod
@@ -35,12 +35,12 @@ steady;
//disp(oo_.mean) ;
estimated_params;
-alp, uniform_pdf,,, 0.0001, 1;
-bet, uniform_pdf,,, 0.75, 0.999;
+alp, uniform_pdf,,, 0.0001, 0.5;
+bet, uniform_pdf,,, 0.0001, 0.99;
tet, uniform_pdf,,, 0.0001, 1;
tau, uniform_pdf,,, 0.0001, 100;
delt, uniform_pdf,,, 0.0001, 0.05;
-rho, uniform_pdf,,, 0.0001, 0.999;
+rho, uniform_pdf,,, 0.8, 0.99;
stderr e_a, uniform_pdf,,, 0.00001, 0.1;
stderr y, uniform_pdf,,, 0.00001, 0.1;
stderr l, uniform_pdf,,, 0.00001, 0.1;
@@ -49,11 +49,11 @@ end;
estimated_params_init;
alp, 0.4;
-bet, 0.99;
+bet, 0.97;
tet, 0.357 ;
tau, 50;
delt, 0.02;
-rho, 0.95;
+rho, 0.9 ;
stderr e_a, .035;
stderr y, .0175;//.00158;
stderr l, .00312;//.0011;
@@ -66,28 +66,42 @@ varobs y l i ;
//options_.gstep(2) = .1;
options_.particle.status = 1;
-options_.particle.algorithm = 'sequential_importance_particle_filter';
options_.particle.initialization = 1;
-options_.particle.pruning = 1;
-options_.particle.number_of_particles = 2000;
+options_.particle.pruning = 0;
+options_.particle.number_of_particles = 1000 ;
options_.particle.resampling.status = 'systematic';
-options_.particle.resampling.neff_threshold = .1;
+
+//options_.particle.resampling.method1 = 'traditional' ;
+//options_.particle.resampling.method1 = 'residual' ;
+options_.particle.resampling.method1 = 'smooth' ;
+
+options_.particle.reampling.method2 = 'kitagawa' ;
+//options_.particle.resampling.method2 = 'stratified' ;
+
+options_.particle.resampling.number_of_partitions = 1;
+options_.particle.resampling.neff_threshold = .5;
set_dynare_threads('local_state_space_iteration_2',3);
options_.particle.algorithm = 'sequential_importance_particle_filter';
//options_.particle.algorithm = 'auxiliary_particle_filter';
//options_.particle.algorithm = 'gaussian_mixture_filter';
-//options_.particle.algorithm = 'each_gaussian_filter';
+//options_.particle.algorithm = 'conditional_particle_filter';
//options_.particle.algorithm = 'gaussian_filter';
//options_.particle.IS_approximation_method = 'quadrature' ;
-//options_.particle.IS_approximation_method = 'cubature' ;
+options_.particle.IS_approximation_method = 'cubature' ;
//options_.particle.IS_approximation_method = 'unscented' ;
//options_.particle.approximation_method = 'quadrature' ;
-//options_.particle.approximation_method = 'cubature' ;
+options_.particle.approximation_method = 'cubature' ;
//options_.particle.approximation_method = 'unscented' ;
//options_.particle.approximation_method = 'MonteCarlo' ;
-estimation(datafile=data_risky_perturb2,nograph,order=2,nobs=100,mh_replic=0,mode_compute=7,mode_check);
+//options_.mh_posterior_mode_estimation=1 ;
+
+// online
+options_.particle.liu_west_delta = 0.9 ;
+options_.mode_check_node_number = 250 ;
+
+estimation(datafile=data_risky_perturb3,nograph,order=2,nobs=100,mh_replic=0,mode_compute=7,mode_check);
--
GitLab