diff --git a/matlab/particle/conditional_particle_filter.m b/matlab/particle/conditional_particle_filter.m index 05b725534eb4b1b2b891c6d08a6f1cf7e30aa0a9..941e01331ecd2c3f625607a509a8b2ed6d1af9b6 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 0c89f0862f07bd6dabd7f0e82a1fe8baa4fbae62..3a770118dd7e1b11255e24226e56d723dce779d3 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 8bd67edf705fe903698ba4d8a8cbcf95761e252e..09ca287b977011bfff201ecce95e2e531b1611eb 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 68a2cd2e9858b8a9207356e3432d521594be8b57..7089a9038977a3b7d1ccde357b8ca68cde94dbe1 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 7ab91083e6aedcdc21b2700075296435e1ae7650..05e5f033fc906c531c9c35f35369da69000222e3 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 b8cc6f92613e6cfb93e4ac4eb00463e4e09cec07..9db9632d55f4bf2e8115ec2dcce94125d1916c53 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 04f7be230fc2150ed3d0105c348c01c7b539bbe7..575d0f16b01d2b0c5a78cc529352a43bf2bdbe90 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 990e310cdb12de896ce463728cebc7b2d1318963..50aacc2c19bdf15da133892d8a3ddeceb79d5f3d 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);