diff --git a/nonlinear-filters/src/auxiliary_initialization.m b/nonlinear-filters/src/auxiliary_initialization.m index 7034882b7..98f79aecb 100644 --- a/nonlinear-filters/src/auxiliary_initialization.m +++ b/nonlinear-filters/src/auxiliary_initialization.m @@ -87,7 +87,7 @@ yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); % 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,ThreadsOptions.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,ThreadsOptions.local_state_space_iteration_2); +tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); %end PredictedObservedMean = weights*(tmp(mf1,:)'); PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); diff --git a/nonlinear-filters/src/auxiliary_particle_filter.m b/nonlinear-filters/src/auxiliary_particle_filter.m index a7934115c..42d57ca45 100644 --- a/nonlinear-filters/src/auxiliary_particle_filter.m +++ b/nonlinear-filters/src/auxiliary_particle_filter.m @@ -1,6 +1,6 @@ function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions) -% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter +% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter % allowing eventually resampling. % % Copyright (C) 2011-2015 Dynare Team @@ -91,11 +91,11 @@ end % [nodes,nodes_weights,nodes_weights_c] = unscented_sigma_points(number_of_structural_innovations,ParticleOptions); % else % error('Estimation: This approximation for the proposal is not implemented or unknown!') -% end +% end % nodes = (Q_lower_triangular_cholesky*(nodes'))' ; nodes = zeros(1,number_of_structural_innovations) ; -nodes_weights = ones(number_of_structural_innovations,1) ; +nodes_weights = ones(number_of_structural_innovations,1) ; for t=1:sample_size yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); @@ -126,7 +126,7 @@ for t=1:sample_size yhat_ = yhat_(:,indx) ; end yhat = yhat(:,indx) ; - weights_stage_1 = weights(indx)./tau_tilde(indx) ; + weights_stage_1 = weights(indx)./tau_tilde(indx) ; epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles); if pruning [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); @@ -137,7 +137,7 @@ for t=1:sample_size StateVectors = tmp(mf0,:); PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); weights_stage_2 = weights_stage_1.*(exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) + 1e-99) ; - lik(t) = log(mean(weights_stage_2)) ; + lik(t) = log(mean(weights_stage_2)) ; weights = weights_stage_2/sum(weights_stage_2); if (ParticleOptions.resampling.status.generic && neff(weights)1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4)) + 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 @@ -106,7 +106,7 @@ if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_a Error = obs - PredictedObservedMean ; StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*Error ; if ParticleOptions.cpf_weights_method.amisanotristani - Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1),PredictedObservedVarianceSquareRoot,Error) ; + Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1),PredictedObservedVarianceSquareRoot,Error) ; end else dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean); @@ -120,15 +120,15 @@ else StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; StateVectorVarianceSquareRoot = chol(StateVectorVariance + eye(number_of_state_variables)*1e-6)' ; if ParticleOptions.cpf_weights_method.amisanotristani - Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1),chol(PredictedObservedVariance)',Error) ; + Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1),chol(PredictedObservedVariance)',Error) ; end end PredictedStateVarianceSquareRoot = chol(PredictedStateVariance + eye(number_of_state_variables)*1e-6)' ; ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ; if ParticleOptions.cpf_weights_method.murrayjonesparslow - Prior = probability2(PredictedStateMean,PredictedStateVarianceSquareRoot,ProposalStateVector) ; - Posterior = probability2(StateVectorMean,StateVectorVarianceSquareRoot,ProposalStateVector) ; - Likelihood = probability2(obs,H_lower_triangular_cholesky,measurement_equations(ProposalStateVector,ReducedForm,ThreadsOptions)) ; + Prior = probability2(PredictedStateMean,PredictedStateVarianceSquareRoot,ProposalStateVector) ; + Posterior = probability2(StateVectorMean,StateVectorVarianceSquareRoot,ProposalStateVector) ; + Likelihood = probability2(obs,H_lower_triangular_cholesky,measurement_equations(ProposalStateVector,ReducedForm,ThreadsOptions)) ; Weights = SampleWeights.*Likelihood.*(Prior./Posterior) ; end diff --git a/nonlinear-filters/src/conditional_particle_filter.m b/nonlinear-filters/src/conditional_particle_filter.m index 89c37b954..01437773e 100644 --- a/nonlinear-filters/src/conditional_particle_filter.m +++ b/nonlinear-filters/src/conditional_particle_filter.m @@ -1,24 +1,24 @@ function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions) -% +% % Evaluates the likelihood of a non-linear model with a particle filter -% - the proposal is built using the Kalman updating step for each particle. -% - we need draws in the errors distributions -% Whether we use Monte-Carlo draws from a multivariate gaussian distribution -% as in Amisano & Tristani (JEDC 2010). -% Whether we use multidimensional Gaussian sparse grids approximations: -% - a univariate Kronrod-Paterson Gaussian quadrature combined by the Smolyak -% operator (ref: Winschel & Kratzig, 2010). +% - the proposal is built using the Kalman updating step for each particle. +% - we need draws in the errors distributions +% Whether we use Monte-Carlo draws from a multivariate gaussian distribution +% as in Amisano & Tristani (JEDC 2010). +% Whether we use multidimensional Gaussian sparse grids approximations: +% - a univariate Kronrod-Paterson Gaussian quadrature combined by the Smolyak +% operator (ref: Winschel & Kratzig, 2010). % - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009a,2009b). -% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1997, van der +% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1997, van der % Merwe & Wan 2003). -% -% Pros: -% - Allows using current observable information in the proposal -% - The use of sparse grids Gaussian approximation is much faster than the Monte-Carlo approach -% Cons: -% - The use of the Kalman updating step may biais the proposal distribution since +% +% Pros: +% - Allows using current observable information in the proposal +% - The use of sparse grids Gaussian approximation is much faster than the Monte-Carlo approach +% Cons: +% - The use of the Kalman updating step may biais the proposal distribution since % it has been derived in a linear context and is implemented in a nonlinear -% context. That is why particle resampling is performed. +% context. That is why particle resampling is performed. % % INPUTS % reduced_form_model [structure] Matlab's structure describing the reduced form model. @@ -58,8 +58,8 @@ function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,ParticleOpt % stephane DOT adjemian AT univ DASH lemans DOT fr persistent init_flag mf1 -persistent number_of_particles -persistent sample_size number_of_observed_variables +persistent number_of_particles +persistent sample_size number_of_observed_variables % Set default if isempty(start) @@ -82,14 +82,14 @@ if isempty(H) H = 0; H_lower_triangular_cholesky = 0; else - H_lower_triangular_cholesky = chol(H)'; + H_lower_triangular_cholesky = chol(H)'; end % Get initial condition for the state vector. StateVectorMean = ReducedForm.StateVectorMean; StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)'; state_variance_rank = size(StateVectorVarianceSquareRoot,2); -Q_lower_triangular_cholesky = chol(Q)'; +Q_lower_triangular_cholesky = chol(Q)'; % Set seed for randn(). set_dynare_seed('default'); @@ -102,13 +102,13 @@ ks = 0 ; StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); SampleWeights = ones(1,number_of_particles)/number_of_particles ; for t=1:sample_size - for i=1:number_of_particles - [StateParticles(:,i),SampleWeights(i)] = ... - conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2) ; + for i=1:number_of_particles + [StateParticles(:,i),SampleWeights(i)] = ... + conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2) ; end SumSampleWeights = sum(SampleWeights) ; - lik(t) = log(SumSampleWeights) ; - SampleWeights = SampleWeights./SumSampleWeights ; + lik(t) = log(SumSampleWeights) ; + SampleWeights = SampleWeights./SumSampleWeights ; if (ParticleOptions.resampling.status.generic && neff(SampleWeights). -[dim,Ndata] = size(X); +[dim,Ndata] = size(X); M = size(StateMu,2) ; if check % Ensure that covariances don't collapse - MIN_COVAR_SQRT = sqrt(eps); - init_covars = StateSqrtP; + 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] = probability3(StateMu,StateSqrtP,StateWeights,X,X_weights); - 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',0); - StateSqrtP(:,:,j) = tcov'; - if check - if min(abs(diag(StateSqrtP(:,:,j)))) < MIN_COVAR_SQRT - StateSqrtP(:,:,j) = init_covars(:,:,j); - end + % Calculate posteriors based on old parameters + [prior,likelihood,marginal,posterior] = probability3(StateMu,StateSqrtP,StateWeights,X,X_weights); + e = sum(log(marginal)); + if (n > 1 && abs((e - eold)/eold) < crit) + return; + else + eold = e; end - end -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',0); + StateSqrtP(:,:,j) = tcov'; + if check + if min(abs(diag(StateSqrtP(:,:,j)))) < MIN_COVAR_SQRT + StateSqrtP(:,:,j) = init_covars(:,:,j); + end + end + end +end diff --git a/nonlinear-filters/src/gaussian_densities.m b/nonlinear-filters/src/gaussian_densities.m index ae3afb130..2246e45e1 100644 --- a/nonlinear-filters/src/gaussian_densities.m +++ b/nonlinear-filters/src/gaussian_densities.m @@ -1,6 +1,6 @@ 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,ThreadsOptions) % -% Elements to calculate the importance sampling ratio +% Elements to calculate the importance sampling ratio % % INPUTS % reduced_form_model [structure] Matlab's structure describing the reduced form model. @@ -36,11 +36,11 @@ function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sq % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -% 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 +% 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,ThreadsOptions) ; eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ; yt_t_1 = sum(yt_t_1_i*weigths1,2) ; @@ -48,5 +48,5 @@ 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 ; +likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ; IncrementalWeights = likelihood.*prior./proposal ; diff --git a/nonlinear-filters/src/gaussian_filter.m b/nonlinear-filters/src/gaussian_filter.m index 2cee18ee7..d7e82df1b 100644 --- a/nonlinear-filters/src/gaussian_filter.m +++ b/nonlinear-filters/src/gaussian_filter.m @@ -112,18 +112,18 @@ for t=1:sample_size if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented StateParticles = bsxfun(@plus,StateVectorMean,StateVectorVarianceSquareRoot*nodes2') ; IncrementalWeights = ... - gaussian_densities(Y(:,t),StateVectorMean,... - StateVectorVarianceSquareRoot,PredictedStateMean,... - PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,... - weights2,weights_c2,ReducedForm,ThreadsOptions) ; + gaussian_densities(Y(:,t),StateVectorMean,... + StateVectorVarianceSquareRoot,PredictedStateMean,... + PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,... + weights2,weights_c2,ReducedForm,ThreadsOptions) ; SampleWeights = weights2.*IncrementalWeights ; - else + else 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,ThreadsOptions) ; + gaussian_densities(Y(:,t),StateVectorMean,... + StateVectorVarianceSquareRoot,PredictedStateMean,... + PredictedStateVarianceSquareRoot,StateParticles,H,const_lik,... + 1/number_of_particles,1/number_of_particles,ReducedForm,ThreadsOptions) ; SampleWeights = IncrementalWeights/number_of_particles ; end SampleWeights = SampleWeights + 1e-6*ones(size(SampleWeights,1),1) ; @@ -132,9 +132,9 @@ for t=1:sample_size SampleWeights = SampleWeights./SumSampleWeights ; if not(ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented) if (ParticleOptions.resampling.status.generic && neff(SampleWeights)1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4)) + 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 diff --git a/nonlinear-filters/src/gaussian_mixture_densities.m b/nonlinear-filters/src/gaussian_mixture_densities.m index 06c3eec1e..7efe1ce01 100644 --- a/nonlinear-filters/src/gaussian_mixture_densities.m +++ b/nonlinear-filters/src/gaussian_mixture_densities.m @@ -1,8 +1,8 @@ function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... - StateMuPost,StateSqrtPPost,StateWeightsPost,... - StateParticles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions) + StateMuPost,StateSqrtPPost,StateWeightsPost,... + StateParticles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions) % -% Elements to calculate the importance sampling ratio +% Elements to calculate the importance sampling ratio % % INPUTS % reduced_form_model [structure] Matlab's structure describing the reduced form model. @@ -38,11 +38,11 @@ function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,State % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -% Compute the density of particles under the prior distribution -[ras,ras,prior] = probability(StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateParticles) ; +% Compute the density of particles under the prior distribution +[ras,ras,prior] = probability(StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateParticles) ; prior = prior' ; -% Compute the density of particles under the proposal distribution -[ras,ras,proposal] = probability(StateMuPost,StateSqrtPPost,StateWeightsPost,StateParticles) ; +% Compute the density of particles under the proposal distribution +[ras,ras,proposal] = probability(StateMuPost,StateSqrtPPost,StateWeightsPost,StateParticles) ; proposal = proposal' ; % Compute the density of the current observation conditionally to each particle yt_t_1_i = measurement_equations(StateParticles,ReducedForm,ThreadsOptions) ; @@ -52,6 +52,5 @@ 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 ; +likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ; IncrementalWeights = likelihood.*prior./proposal ; - diff --git a/nonlinear-filters/src/gaussian_mixture_filter.m b/nonlinear-filters/src/gaussian_mixture_filter.m index 0dc257c8c..eaab0070f 100644 --- a/nonlinear-filters/src/gaussian_mixture_filter.m +++ b/nonlinear-filters/src/gaussian_mixture_filter.m @@ -63,15 +63,15 @@ 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_structural_innovations = length(ReducedForm.Q); - G = ParticleOptions.mixture_state_variables; % number of GM components in state - number_of_particles = ParticleOptions.number_of_particles; - init_flag = 1; + 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_structural_innovations = length(ReducedForm.Q); + G = ParticleOptions.mixture_state_variables; % number of GM components in state + number_of_particles = ParticleOptions.number_of_particles; + init_flag = 1; end % compute gaussian quadrature nodes and weights on states and shocks @@ -104,14 +104,14 @@ else end Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)'; -% Initialize mixtures +% Initialize mixtures StateWeights = ones(1,G)/G ; StateMu = ReducedForm.StateVectorMean ; StateSqrtP = zeros(number_of_state_variables,number_of_state_variables,G) ; temp = reduced_rank_cholesky(ReducedForm.StateVectorVariance)' ; StateMu = bsxfun(@plus,StateMu,bsxfun(@times,diag(temp),(-(G-1)/2:1:(G-1)/2))/10) ; for g=1:G - StateSqrtP(:,:,g) = temp/sqrt(G) ; + StateSqrtP(:,:,g) = temp/sqrt(G) ; end % if ParticleOptions.mixture_structural_shocks==1 @@ -135,11 +135,11 @@ end % for i=1:I % StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky/sqrt(StructuralShocksWeights(i)) ; % end -% -% if ParticleOptions.mixture_measurement_shocks==1 +% +% if ParticleOptions.mixture_measurement_shocks==1 % ObservationShocksMu = zeros(1,number_of_observed_variables) ; % ObservationShocksWeights = 1 ; -% else +% else % if ParticleOptions.proposal_approximation.cubature % [ObservationShocksMu,ObservationShocksWeights] = spherical_radial_sigma_points(number_of_observed_variables); % ObservationShocksWeights = ones(size(ObservationShocksMu,1),1)*ObservationShocksWeights; @@ -150,7 +150,7 @@ end % error('Estimation: This approximation for the proposal is not implemented or unknown!') % end % end -% end +% end % J = size(ObservationShocksWeights,1) ; % ObservationShocksMu = H_lower_triangular_cholesky*(ObservationShocksMu') ; % ObservationShocksSqrtP = zeros(number_of_observed_variables,number_of_observed_variables,J) ; @@ -180,7 +180,7 @@ elseif ParticleOptions.mixture_structural_shocks==1 StructuralShocksMu = Q_lower_triangular_cholesky*(StructuralShocksMu') ; StructuralShocksSqrtP = zeros(number_of_structural_innovations,number_of_structural_innovations,I) ; for i=1:I - StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky ; + StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky ; end else if ParticleOptions.proposal_approximation.cubature @@ -197,7 +197,7 @@ else StructuralShocksMu = Q_lower_triangular_cholesky*(StructuralShocksMu') ; StructuralShocksSqrtP = zeros(number_of_structural_innovations,number_of_structural_innovations,I) ; for i=1:I - StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky/sqrt(StructuralShocksWeights(i)) ; + StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky/sqrt(StructuralShocksWeights(i)) ; end end @@ -208,14 +208,14 @@ ObservationShocksMu = H_lower_triangular_cholesky*(ObservationShocksMu') ; ObservationShocksSqrtP = zeros(number_of_observed_variables,number_of_observed_variables,J) ; ObservationShocksSqrtP(:,:,1) = H_lower_triangular_cholesky ; -% if ParticleOptions.mixture_measurement_shocks==0 +% if ParticleOptions.mixture_measurement_shocks==0 % ObservationShocksMu = zeros(1,number_of_observed_variables) ; % ObservationShocksWeights = 1 ; % J = 1 ; % ObservationShocksMu = H_lower_triangular_cholesky*(ObservationShocksMu') ; % ObservationShocksSqrtP = zeros(number_of_observed_variables,number_of_observed_variables,J) ; % ObservationShocksSqrtP(:,:,1) = H_lower_triangular_cholesky ; -% elseif ParticleOptions.mixture_measurement_shocks==1 +% elseif ParticleOptions.mixture_measurement_shocks==1 % if ParticleOptions.proposal_approximation.cubature % [ObservationShocksMu,ObservationShocksWeights] = spherical_radial_sigma_points(number_of_observed_variables); % ObservationShocksWeights = ones(size(ObservationShocksMu,1),1)*ObservationShocksWeights; @@ -232,7 +232,7 @@ ObservationShocksSqrtP(:,:,1) = H_lower_triangular_cholesky ; % for j=1:J % ObservationShocksSqrtP(:,:,j) = H_lower_triangular_cholesky ; % end -% else +% else % if ParticleOptions.proposal_approximation.cubature % [ObservationShocksMu,ObservationShocksWeights] = spherical_radial_sigma_points(number_of_observed_variables); % ObservationShocksWeights = ones(size(ObservationShocksMu,1),1)*ObservationShocksWeights; @@ -277,10 +277,10 @@ for t=1:sample_size gsecond = gprime + (j-1)*Gprime ; [StateMuPrior(:,gprime),StateSqrtPPrior(:,:,gprime),StateWeightsPrior(1,gprime),... StateMuPost(:,gsecond),StateSqrtPPost(:,:,gsecond),StateWeightsPost(1,gsecond)] =... - gaussian_mixture_filter_bank(ReducedForm,Y(:,t),StateMu(:,g),StateSqrtP(:,:,g),StateWeights(g),... - StructuralShocksMu(:,i),StructuralShocksSqrtP(:,:,i),StructuralShocksWeights(i),... - ObservationShocksMu(:,j),ObservationShocksSqrtP(:,:,j),ObservationShocksWeights(j),... - H,H_lower_triangular_cholesky,const_lik,ParticleOptions,ThreadsOptions) ; + gaussian_mixture_filter_bank(ReducedForm,Y(:,t),StateMu(:,g),StateSqrtP(:,:,g),StateWeights(g),... + StructuralShocksMu(:,i),StructuralShocksSqrtP(:,:,i),StructuralShocksWeights(i),... + ObservationShocksMu(:,j),ObservationShocksSqrtP(:,:,j),ObservationShocksWeights(j),... + H,H_lower_triangular_cholesky,const_lik,ParticleOptions,ThreadsOptions) ; end end end @@ -293,8 +293,8 @@ for t=1:sample_size for i=1:Gsecond StateParticles = bsxfun(@plus,StateMuPost(:,i),StateSqrtPPost(:,:,i)*nodes') ; IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... - StateMuPost,StateSqrtPPost,StateWeightsPost,... - StateParticles,H,const_lik,weights,weights_c,ReducedForm,ThreadsOptions) ; + StateMuPost,StateSqrtPPost,StateWeightsPost,... + StateParticles,H,const_lik,weights,weights_c,ReducedForm,ThreadsOptions) ; SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights) ; end SumSampleWeights = sum(SampleWeights) ; @@ -311,9 +311,9 @@ for t=1:sample_size % Sample particle in the proposal distribution, ie the posterior state GM StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles) ; IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... - StateMuPost,StateSqrtPPost,StateWeightsPost,... - StateParticles,H,const_lik,1/number_of_particles,... - 1/number_of_particles,ReducedForm,ThreadsOptions) ; + StateMuPost,StateSqrtPPost,StateWeightsPost,... + StateParticles,H,const_lik,1/number_of_particles,... + 1/number_of_particles,ReducedForm,ThreadsOptions) ; SampleWeights = IncrementalWeights/number_of_particles ; SumSampleWeights = sum(SampleWeights,1) ; SampleWeights = SampleWeights./SumSampleWeights ; diff --git a/nonlinear-filters/src/gaussian_mixture_filter_bank.m b/nonlinear-filters/src/gaussian_mixture_filter_bank.m index 5bcf702cb..017192712 100644 --- a/nonlinear-filters/src/gaussian_mixture_filter_bank.m +++ b/nonlinear-filters/src/gaussian_mixture_filter_bank.m @@ -1,12 +1,12 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPPost,StateWeightsPost] =... - gaussian_mixture_filter_bank(ReducedForm,obs,StateMu,StateSqrtP,StateWeights,... - StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,... - ObservationShocksMu,ObservationShocksSqrtP,ObservationShocksWeights,... - H,H_lower_triangular_cholesky,normfactO,ParticleOptions,ThreadsOptions) + gaussian_mixture_filter_bank(ReducedForm,obs,StateMu,StateSqrtP,StateWeights,... + StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,... + ObservationShocksMu,ObservationShocksSqrtP,ObservationShocksWeights,... + H,H_lower_triangular_cholesky,normfactO,ParticleOptions,ThreadsOptions) % % Computes the proposal with a gaussian approximation for importance -% sampling -% This proposal is a gaussian distribution calculated à la Kalman +% sampling +% This proposal is a gaussian distribution calculated à la Kalman % % INPUTS % reduced_form_model [structure] Matlab's structure describing the reduced form model. @@ -43,8 +43,8 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPP persistent init_flag2 mf0 mf1 %nodes3 weights3 weights_c3 -persistent number_of_state_variables number_of_observed_variables -persistent number_of_structural_innovations +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; @@ -55,8 +55,8 @@ 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)) + 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 diff --git a/nonlinear-filters/src/importance_sampling.m b/nonlinear-filters/src/importance_sampling.m index 403b6beb3..b4a6daf4a 100644 --- a/nonlinear-filters/src/importance_sampling.m +++ b/nonlinear-filters/src/importance_sampling.m @@ -17,13 +17,12 @@ function State_Particles = importance_sampling(StateMuPost,StateSqrtPPost,StateW % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -[Xdim,Gsecond] = size(StateMuPost) ; +[Xdim,Gsecond] = size(StateMuPost) ; u = rand(numP,1); -[Nc,comp] = histc(u, cumsum([0; StateWeightsPost])); +[Nc,comp] = histc(u, cumsum([0; StateWeightsPost])); State_Particles = zeros(Xdim,numP); for k=1:Gsecond - idx = bsxfun(@eq,comp,k*ones(size(comp))) ; - State_Particles(:,idx) = StateSqrtPPost(:,:,k)*randn(Xdim,Nc(k)); + idx = bsxfun(@eq,comp,k*ones(size(comp))) ; + State_Particles(:,idx) = StateSqrtPPost(:,:,k)*randn(Xdim,Nc(k)); end -State_Particles= State_Particles + StateMuPost(:,comp); - +State_Particles= State_Particles + StateMuPost(:,comp); diff --git a/nonlinear-filters/src/measurement_equations.m b/nonlinear-filters/src/measurement_equations.m index d38c59aea..7cf4d9320 100644 --- a/nonlinear-filters/src/measurement_equations.m +++ b/nonlinear-filters/src/measurement_equations.m @@ -1,4 +1,4 @@ -function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions) +function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions) % Copyright (C) 2013 Dynare Team % @@ -28,6 +28,3 @@ state_variables_steady_state = ReducedForm.state_variables_steady_state; number_of_structural_innovations = length(ReducedForm.Q); yhat = bsxfun(@minus,StateVectors,state_variables_steady_state) ; measure = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,size(yhat,2)),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); - - - diff --git a/nonlinear-filters/src/multivariate_smooth_resampling.m b/nonlinear-filters/src/multivariate_smooth_resampling.m index 09ca287b9..e62111cc1 100644 --- a/nonlinear-filters/src/multivariate_smooth_resampling.m +++ b/nonlinear-filters/src/multivariate_smooth_resampling.m @@ -63,10 +63,10 @@ number_of_states = size(particles,2); [P,D] = eig(particles'*(bsxfun(@times,1/number_of_particles,particles))) ; D = diag(D) ; vectors = bsxfun(@times,P,sqrt(D)') ; -orthogonalized_particles = bsxfun(@rdivide,particles*vectors,D') ; +orthogonalized_particles = bsxfun(@rdivide,particles*vectors,D') ; new_particles = zeros(number_of_particles,number_of_states) ; for j=1:number_of_states - tout = sortrows( [orthogonalized_particles(:,j) weights],1) ; - new_particles(:,j) = univariate_smooth_resampling(tout(:,2),tout(:,1),number_of_particles) ; + tout = sortrows( [orthogonalized_particles(:,j) weights],1) ; + new_particles(:,j) = univariate_smooth_resampling(tout(:,2),tout(:,1),number_of_particles) ; end new_particles = new_particles*(vectors') ; diff --git a/nonlinear-filters/src/mykmeans.m b/nonlinear-filters/src/mykmeans.m index 56af5390b..4506aca25 100644 --- a/nonlinear-filters/src/mykmeans.m +++ b/nonlinear-filters/src/mykmeans.m @@ -1,4 +1,4 @@ -function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod) +function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod) % Copyright (C) 2013 Dynare Team % @@ -20,34 +20,34 @@ function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod) [n,m] = size(x) ; indold = zeros(1,m) ; if cod==0 - d = transpose(sum(bsxfun(@power,bsxfun(@minus,x,mean(x)),2))); - d = sortrows( [transpose(1:m) d],2) ; - d = d((1+(0:1:g-1))*m/g,1) ; - c = x(:,d); + d = transpose(sum(bsxfun(@power,bsxfun(@minus,x,mean(x)),2))); + d = sortrows( [transpose(1:m) d],2) ; + d = d((1+(0:1:g-1))*m/g,1) ; + c = x(:,d); else - c = init ; -end -for iter=1:300 - dist = zeros(g,m) ; - for i=1:g - dist(i,:) = sum(bsxfun(@power,bsxfun(@minus,x,c(:,i)),2)); - end - [rien,ind] = min(dist) ; - if isequal(ind,indold) - break ; - end - indold = ind ; - for i=1:g - lin = bsxfun(@eq,ind,i.*ones(1,m)) ; - h = x(:,lin) ; - c(:,i) = mean(h,2) ; - end + c = init ; end -SqrtVariance = zeros(n,n,g) ; -Weights = zeros(1,g) ; +for iter=1:300 + dist = zeros(g,m) ; + for i=1:g + dist(i,:) = sum(bsxfun(@power,bsxfun(@minus,x,c(:,i)),2)); + end + [rien,ind] = min(dist) ; + if isequal(ind,indold) + break ; + end + indold = ind ; + for i=1:g + lin = bsxfun(@eq,ind,i.*ones(1,m)) ; + h = x(:,lin) ; + c(:,i) = mean(h,2) ; + end +end +SqrtVariance = zeros(n,n,g) ; +Weights = zeros(1,g) ; for i=1:g - temp = x(:,bsxfun(@eq,ind,i*ones(1,m))) ; - u = bsxfun(@minus,temp,mean(temp,2)); %temp-mean(temp,1)' ; - SqrtVariance(:,:,i) = chol( (u*u')/size(temp,2) )' ; - Weights(i) = size(temp,2)/m ; + temp = x(:,bsxfun(@eq,ind,i*ones(1,m))) ; + u = bsxfun(@minus,temp,mean(temp,2)); %temp-mean(temp,1)' ; + SqrtVariance(:,:,i) = chol( (u*u')/size(temp,2) )' ; + Weights(i) = size(temp,2)/m ; end diff --git a/nonlinear-filters/src/nonlinear_kalman_filter.m b/nonlinear-filters/src/nonlinear_kalman_filter.m index dd61a5e50..608fbacec 100644 --- a/nonlinear-filters/src/nonlinear_kalman_filter.m +++ b/nonlinear-filters/src/nonlinear_kalman_filter.m @@ -10,9 +10,9 @@ function [LIK,lik] = nonlinear_kalman_filter(ReducedForm, Y, start, ParticleOpti % from the resulting nodes/particles and allows to generate new distributions at the % following observation. % Pros: The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles -% filters since it treats a lesser number of particles. -% Cons: 1. Application a linear projection formulae in a nonlinear context. -% 2. Parameter estimations may be biaised if the model is truly non-gaussian since predictive and +% filters since it treats a lesser number of particles. +% Cons: 1. Application a linear projection formulae in a nonlinear context. +% 2. Parameter estimations may be biaised if the model is truly non-gaussian since predictive and % filtered densities are unimodal. % % INPUTS @@ -47,7 +47,7 @@ function [LIK,lik] = nonlinear_kalman_filter(ReducedForm, Y, start, ParticleOpti % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -persistent init_flag mf0 mf1 nodes weights weights_c +persistent init_flag mf0 mf1 nodes weights weights_c persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations % Set default @@ -64,8 +64,8 @@ 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)) + 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 @@ -100,12 +100,12 @@ end if ParticleOptions.distribution_approximation.montecarlo set_dynare_seed('default'); -end +end % Get covariance matrices H = ReducedForm.H; -H_lower_triangular_cholesky = chol(H)' ; -Q_lower_triangular_cholesky = chol(ReducedForm.Q)'; +H_lower_triangular_cholesky = chol(H)' ; +Q_lower_triangular_cholesky = chol(ReducedForm.Q)'; % Get initial condition for the state vector. StateVectorMean = ReducedForm.StateVectorMean; diff --git a/nonlinear-filters/src/online_auxiliary_filter.m b/nonlinear-filters/src/online_auxiliary_filter.m index 8018ca69c..a14548559 100644 --- a/nonlinear-filters/src/online_auxiliary_filter.m +++ b/nonlinear-filters/src/online_auxiliary_filter.m @@ -1,5 +1,5 @@ function [xparam,std_param,lb_95,ub_95,median_param] = online_auxiliary_filter(xparam1,DynareDataset,dataset_info,DynareOptions,Model,EstimatedParameters,BayesInfo,bounds,DynareResults) - + % Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters. % % INPUTS @@ -47,8 +47,8 @@ second_resample = 0 ; variance_update = 1 ; % initialization of state particles -[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... - solve_model_for_online_filter(1,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ; +[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... + solve_model_for_online_filter(1,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ; % Set persistent variables. if isempty(init_flag) @@ -61,11 +61,11 @@ 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 ; - start_param = xparam1 ; + start_param = xparam1 ; init_flag = 1; end -% Get initial conditions for the state particles +% Get initial conditions for the state particles StateVectorMean = ReducedForm.StateVectorMean; StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)'; state_variance_rank = size(StateVectorVarianceSquareRoot,2); @@ -73,13 +73,13 @@ StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_r if pruning StateVectors_ = StateVectors; end - -% parameters for the Liu & West filter + +% parameters for the Liu & West filter h_square = (3*liu_west_delta-1)/(2*liu_west_delta) ; h_square = 1-h_square*h_square ; small_a = sqrt(1-h_square) ; -% Initialization of parameter particles +% Initialization of parameter particles xparam = zeros(number_of_parameters,number_of_particles) ; stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/100 ; stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/50 ; @@ -107,14 +107,14 @@ std_xparam = zeros(number_of_parameters,sample_size) ; lb95_xparam = zeros(number_of_parameters,sample_size) ; ub95_xparam = zeros(number_of_parameters,sample_size) ; -%% The Online filter +%% The Online filter for t=1:sample_size if t>1 fprintf('\nSubsample with %s first observations.\n\n', int2str(t)) else fprintf('\nSubsample with only the first observation.\n\n', int2str(t)) end - % Moments of parameters particles distribution + % Moments of parameters particles distribution m_bar = xparam*(weights') ; temp = bsxfun(@minus,xparam,m_bar) ; sigma_bar = (bsxfun(@times,weights,temp))*(temp') ; @@ -124,8 +124,8 @@ for t=1:sample_size % Prediction (without shocks) tau_tilde = zeros(1,number_of_particles) ; for i=1:number_of_particles - % model resolution - [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... + % model resolution + [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... solve_model_for_online_filter(t,xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ; steadystate = ReducedForm.steadystate; state_variables_steady_state = ReducedForm.state_variables_steady_state; @@ -136,7 +136,7 @@ for t=1:sample_size ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; - % particle likelihood contribution + % particle likelihood contribution yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state); if pruning yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state); @@ -151,7 +151,7 @@ for t=1:sample_size tau_tilde(i) = weights(i).*(tpdf(z,3*ones(size(z)))+1e-99) ; %tau_tilde(i) = weights(i).*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError),1))) ; end - % particles selection + % particles selection tau_tilde = tau_tilde/sum(tau_tilde) ; indx = resample(0,tau_tilde',DynareOptions.particle); StateVectors = StateVectors(:,indx) ; @@ -160,7 +160,7 @@ for t=1:sample_size end xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam(:,indx)) ; w_stage1 = weights(indx)./tau_tilde(indx) ; - % draw in the new distributions + % draw in the new distributions wtilde = zeros(1,number_of_particles) ; i = 1 ; while i<=number_of_particles @@ -179,9 +179,9 @@ for t=1:sample_size ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; - % Get covariance matrices and structural shocks + % Get covariance matrices and structural shocks epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations,1) ; - % compute particles likelihood contribution + % compute particles likelihood contribution yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state); if pruning yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state); @@ -196,13 +196,13 @@ for t=1:sample_size i = i+1 ; end end - % normalization + % normalization weights = wtilde/sum(wtilde); if (variance_update==1) && (neff(weights)=0.025 && pass1==1 - lb95_xparam(i,t) = temp(j,1) ; - pass1 = 2 ; - end - if cumulated_weights(j)>=0.5 && pass2==1 - median_xparam(i,t) = temp(j,1) ; - pass2 = 2 ; - end - if cumulated_weights(j)>=0.975 && pass3==1 - ub95_xparam(i,t) = temp(j,1) ; - pass3 = 2 ; - end - end + temp = sortrows([xparam(i,:)' weights'],1) ; + cumulated_weights = cumsum(temp(:,2)) ; + pass1=1 ; + pass2=1 ; + pass3=1 ; + for j=1:number_of_particles + if cumulated_weights(j)>=0.025 && pass1==1 + lb95_xparam(i,t) = temp(j,1) ; + pass1 = 2 ; + end + if cumulated_weights(j)>=0.5 && pass2==1 + median_xparam(i,t) = temp(j,1) ; + pass2 = 2 ; + end + if cumulated_weights(j)>=0.975 && pass3==1 + ub95_xparam(i,t) = temp(j,1) ; + pass3 = 2 ; + end + end end end str = sprintf(' Lower Bound (95%%) \t Mean \t\t\t Upper Bound (95%%)'); @@ -262,7 +262,7 @@ lb_95 = lb95_xparam(:,sample_size) ; ub_95 = ub95_xparam(:,sample_size) ; median_param = median_xparam(:,sample_size) ; -%% Plot parameters trajectory +%% Plot parameters trajectory TeX = DynareOptions.TeX; [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters); @@ -322,7 +322,7 @@ end %% Plot Parameter Densities number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two. bandwidth = 0; % Rule of thumb optimal bandwidth parameter. -kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation. +kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation. for plt = 1:nbplt, if TeX NAMES = []; @@ -366,5 +366,4 @@ for plt = 1:nbplt, fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,' \n'); end -end - +end diff --git a/nonlinear-filters/src/probability.m b/nonlinear-filters/src/probability.m index 5d76e49ee..94f09b9cb 100644 --- a/nonlinear-filters/src/probability.m +++ b/nonlinear-filters/src/probability.m @@ -17,23 +17,23 @@ function [prior,likelihood,C,posterior] = probability(mu,sqrtP,prior,X) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -[dim,nov] = size(X); +[dim,nov] = size(X); M = size(mu,2) ; if nargout>1 - likelihood = zeros(M,nov); - normfact = (2*pi)^(dim/2); - for k=1:M - XX = bsxfun(@minus,X,mu(:,k)); - S = sqrtP(:,:,k); - foo = S \ XX; - likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S)))); - end + likelihood = zeros(M,nov); + normfact = (2*pi)^(dim/2); + for k=1:M + XX = bsxfun(@minus,X,mu(:,k)); + S = sqrtP(:,:,k); + foo = S \ XX; + likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S)))); + end end likelihood = likelihood + 1e-99; if nargout>2 - C = prior*likelihood + 1e-99; + C = prior*likelihood + 1e-99; end if nargout>3 - posterior = bsxfun(@rdivide,bsxfun(@times,prior',likelihood),C) + 1e-99 ; - posterior = bsxfun(@rdivide,posterior,sum(posterior,1)); + posterior = bsxfun(@rdivide,bsxfun(@times,prior',likelihood),C) + 1e-99 ; + posterior = bsxfun(@rdivide,posterior,sum(posterior,1)); end diff --git a/nonlinear-filters/src/probability2.m b/nonlinear-filters/src/probability2.m index a7486d1b8..dfc71ff24 100644 --- a/nonlinear-filters/src/probability2.m +++ b/nonlinear-filters/src/probability2.m @@ -1,7 +1,7 @@ -function [density] = probability2(mu,S,X) +function [density] = probability2(mu,S,X) % % Multivariate gaussian density -% +% % INPUTS % n [integer] scalar, number of variables. % @@ -10,10 +10,10 @@ function [density] = probability2(mu,S,X) % weigths [double] associated weigths % % REFERENCES -% +% % % NOTES -% +% % Copyright (C) 2009-2012 Dynare Team % % This file is part of Dynare. @@ -31,7 +31,7 @@ function [density] = probability2(mu,S,X) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -dim = size(X,1) ; -normfact = bsxfun(@power,(2*pi),(dim/2)) ; +dim = size(X,1) ; +normfact = bsxfun(@power,(2*pi),(dim/2)) ; foo = S\(bsxfun(@minus,X,mu)) ; density = exp(-0.5*sum(foo.*foo)')./abs((normfact*prod(diag(S)))) + 1e-99 ; \ No newline at end of file diff --git a/nonlinear-filters/src/probability3.m b/nonlinear-filters/src/probability3.m index 49a14b429..5e2a9b4a8 100644 --- a/nonlinear-filters/src/probability3.m +++ b/nonlinear-filters/src/probability3.m @@ -17,23 +17,23 @@ function [prior,likelihood,C,posterior] = probability3(mu,sqrtP,prior,X,X_weight % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -[dim,nov] = size(X); +[dim,nov] = size(X); M = size(mu,2) ; if nargout>1 - likelihood = zeros(M,nov); - normfact = (2*pi)^(dim/2); - for k=1:M - XX = bsxfun(@minus,X,mu(:,k)); - S = sqrtP(:,:,k); - foo = S \ XX; - likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S)))); - end + likelihood = zeros(M,nov); + normfact = (2*pi)^(dim/2); + for k=1:M + XX = bsxfun(@minus,X,mu(:,k)); + S = sqrtP(:,:,k); + foo = S \ XX; + likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S)))); + end end wlikelihood = bsxfun(@times,X_weights,likelihood) + 1e-99; if nargout>2 - C = prior*wlikelihood + 1e-99; + C = prior*wlikelihood + 1e-99; end if nargout>3 - posterior = bsxfun(@rdivide,bsxfun(@times,prior',wlikelihood),C) + 1e-99 ; - posterior = bsxfun(@rdivide,posterior,sum(posterior,1)); + posterior = bsxfun(@rdivide,bsxfun(@times,prior',wlikelihood),C) + 1e-99 ; + posterior = bsxfun(@rdivide,posterior,sum(posterior,1)); end diff --git a/nonlinear-filters/src/residual_resampling.m b/nonlinear-filters/src/residual_resampling.m index c2f09e3b8..9a7d8cce8 100644 --- a/nonlinear-filters/src/residual_resampling.m +++ b/nonlinear-filters/src/residual_resampling.m @@ -76,33 +76,33 @@ iWEIGHTS = fix(WEIGHTS); number_of_trials = number_of_particles-sum(iWEIGHTS); if number_of_trials - WEIGHTS = (WEIGHTS-iWEIGHTS)/number_of_trials; - EmpiricalCDF = cumsum(WEIGHTS); - if kitagawa_resampling - u = (transpose(1:number_of_trials)-1+noise(:))/number_of_trials; - else - u = fliplr(cumprod(noise(1:number_of_trials).^(1./(number_of_trials:-1:1)))); - end - j=1; - for i=1:number_of_trials - while (u(i)>EmpiricalCDF(j)) - j=j+1; + WEIGHTS = (WEIGHTS-iWEIGHTS)/number_of_trials; + EmpiricalCDF = cumsum(WEIGHTS); + if kitagawa_resampling + u = (transpose(1:number_of_trials)-1+noise(:))/number_of_trials; + else + u = fliplr(cumprod(noise(1:number_of_trials).^(1./(number_of_trials:-1:1)))); end - iWEIGHTS(j)=iWEIGHTS(j)+1; - if kitagawa_resampling==0 - j=1; + j=1; + for i=1:number_of_trials + while (u(i)>EmpiricalCDF(j)) + j=j+1; + end + iWEIGHTS(j)=iWEIGHTS(j)+1; + if kitagawa_resampling==0 + j=1; + end end - end end k=1; for i=1:number_of_particles - if (iWEIGHTS(i)>0) - for j=k:k+iWEIGHTS(i)-1 - indx(j) = jndx(i); + if (iWEIGHTS(i)>0) + for j=k:k+iWEIGHTS(i)-1 + indx(j) = jndx(i); + end end - end - k = k + iWEIGHTS(i); + k = k + iWEIGHTS(i); end if particles==0 diff --git a/nonlinear-filters/src/solve_model_for_online_filter.m b/nonlinear-filters/src/solve_model_for_online_filter.m index 62ed99cb5..1ea2e67a4 100644 --- a/nonlinear-filters/src/solve_model_for_online_filter.m +++ b/nonlinear-filters/src/solve_model_for_online_filter.m @@ -185,18 +185,18 @@ if EstimatedParameters.ncx Q(k2,k1) = Q(k1,k2); end % Try to compute the cholesky decomposition of Q (possible iff Q is positive definite) -% [CholQ,testQ] = chol(Q); -% if testQ - % The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty. -% a = diag(eig(Q)); -% k = find(a < 0); -% if k > 0 -% fval = objective_function_penalty_base+sum(-a(k)); -% exit_flag = 0; -% info = 43; -% return -% end -% end + % [CholQ,testQ] = chol(Q); + % if testQ + % The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty. + % a = diag(eig(Q)); + % k = find(a < 0); + % if k > 0 + % fval = objective_function_penalty_base+sum(-a(k)); + % exit_flag = 0; + % info = 43; + % return + % end + % end offset = offset+EstimatedParameters.ncx; end @@ -210,18 +210,18 @@ if EstimatedParameters.ncn H(k2,k1) = H(k1,k2); end % Try to compute the cholesky decomposition of H (possible iff H is positive definite) -% [CholH,testH] = chol(H); -% if testH - % The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty. -% a = diag(eig(H)); -% k = find(a < 0); -% if k > 0 -% fval = objective_function_penalty_base+sum(-a(k)); -% exit_flag = 0; -% info = 44; -% return -% end -% end + % [CholH,testH] = chol(H); + % if testH + % The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty. + % a = diag(eig(H)); + % k = find(a < 0); + % if k > 0 + % fval = objective_function_penalty_base+sum(-a(k)); + % exit_flag = 0; + % info = 44; + % return + % end + % end offset = offset+EstimatedParameters.ncn; end @@ -244,7 +244,7 @@ Model.H = H; %disp(info) if info(1) ~= 0 - ReducedForm = 0 ; + ReducedForm = 0 ; exit_flag = 55; return end @@ -312,7 +312,7 @@ if DynareOptions.order>1 ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:); ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:); ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx); -else +else ReducedForm.ghxx = zeros(size(restrict_variables_idx,1),size(dr.kstate,2)); ReducedForm.ghuu = zeros(size(restrict_variables_idx,1),size(dr.ghu,2)); ReducedForm.ghxu = zeros(size(restrict_variables_idx,1),size(dr.ghx,2)); @@ -325,7 +325,7 @@ ReducedForm.mf0 = mf0; ReducedForm.mf1 = mf1; % Set initial condition for t=1 -if observation_number==1 +if observation_number==1 switch DynareOptions.particle.initialization case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model. StateVectorMean = ReducedForm.constant(mf0); diff --git a/nonlinear-filters/src/spherical_radial_sigma_points.m b/nonlinear-filters/src/spherical_radial_sigma_points.m index 754caacfa..db71f6119 100644 --- a/nonlinear-filters/src/spherical_radial_sigma_points.m +++ b/nonlinear-filters/src/spherical_radial_sigma_points.m @@ -1,7 +1,7 @@ -function [nodes,weights] = spherical_radial_sigma_points(n) +function [nodes,weights] = spherical_radial_sigma_points(n) % -% Computes nodes and weigths from a third-degree spherical-radial cubature -% rule. +% Computes nodes and weigths from a third-degree spherical-radial cubature +% rule. % INPUTS % n [integer] scalar, number of variables. % @@ -11,10 +11,10 @@ function [nodes,weights] = spherical_radial_sigma_points(n) % % REFERENCES % -% Arasaratnam & Haykin 2008,2009. +% Arasaratnam & Haykin 2008,2009. % % NOTES -% +% % Copyright (C) 2009-2012 Dynare Team % % This file is part of Dynare. diff --git a/nonlinear-filters/src/traditional_resampling.m b/nonlinear-filters/src/traditional_resampling.m index a9d86eb36..87f550368 100644 --- a/nonlinear-filters/src/traditional_resampling.m +++ b/nonlinear-filters/src/traditional_resampling.m @@ -75,7 +75,7 @@ c = cumsum(weights); % Draw a starting point. if kitagawa_resampling randvec = (transpose(1:number_of_particles)-1+noise(:))/number_of_particles ; -else +else randvec = fliplr(cumprod(noise.^(1./(number_of_particles:-1:1)))); end diff --git a/nonlinear-filters/src/univariate_smooth_resampling.m b/nonlinear-filters/src/univariate_smooth_resampling.m index 8238902a1..e89b2cc6d 100644 --- a/nonlinear-filters/src/univariate_smooth_resampling.m +++ b/nonlinear-filters/src/univariate_smooth_resampling.m @@ -61,16 +61,16 @@ lambda_tilde = [ (.5*weights(1)) ; lambda_bar = cumsum(lambda_tilde) ; u = rand(1,1) ; new_particles = zeros(number_of_new_particles,1) ; -rj = 0 ; +rj = 0 ; i = 1 ; j = 1 ; while i<=number_of_new_particles u_j = ( i-1 + u)/number_of_new_particles ; while u_j>lambda_bar(j) - rj = j ; + rj = j ; j = j+1 ; end - if rj==0 + if rj==0 new_particles(i) = particles(1) ; elseif rj==M new_particles(i) = particles(M) ; diff --git a/nonlinear-filters/src/unscented_sigma_points.m b/nonlinear-filters/src/unscented_sigma_points.m index 38771758b..d1cba6eb5 100644 --- a/nonlinear-filters/src/unscented_sigma_points.m +++ b/nonlinear-filters/src/unscented_sigma_points.m @@ -1,6 +1,6 @@ -function [nodes,W_m,W_c] = unscented_sigma_points(n,ParticleOptions) +function [nodes,W_m,W_c] = unscented_sigma_points(n,ParticleOptions) % -% Computes nodes and weigths for a scaled unscented transform cubature +% Computes nodes and weigths for a scaled unscented transform cubature % INPUTS % n [integer] scalar, number of variables. % @@ -10,10 +10,10 @@ function [nodes,W_m,W_c] = unscented_sigma_points(n,ParticleOptions) % % REFERENCES % -% +% % % NOTES -% +% % Copyright (C) 2009-2012 Dynare Team % % This file is part of Dynare. @@ -31,12 +31,10 @@ function [nodes,W_m,W_c] = unscented_sigma_points(n,ParticleOptions) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -lambda = (ParticleOptions.unscented.alpha^2)*(n+ParticleOptions.unscented.kappa) - n ; +lambda = (ParticleOptions.unscented.alpha^2)*(n+ParticleOptions.unscented.kappa) - n ; nodes = [ zeros(n,1) ( sqrt(n+lambda).*([ eye(n) -eye(n)]) ) ]' ; W_m = lambda/(n+lambda) ; W_c = W_m + (1-ParticleOptions.unscented.alpha^2+ParticleOptions.unscented.beta) ; temp = ones(2*n,1)/(2*(n+lambda)) ; -W_m = [W_m ; temp] ; -W_c = [W_c ; temp] ; - - +W_m = [W_m ; temp] ; +W_c = [W_c ; temp] ;