diff --git a/matlab/estimation/online/online_auxiliary_filter.m b/matlab/estimation/online/online_auxiliary_filter.m
index ebe1394ff1d2e18c98bfd4e4e69aa3f5f848e7f9..f4945f491f57cf98ff706546f32519b8e976d73f 100644
--- a/matlab/estimation/online/online_auxiliary_filter.m
+++ b/matlab/estimation/online/online_auxiliary_filter.m
@@ -41,7 +41,7 @@ SimulationFolder = CheckPath('online', M_.dname);
bounds = prior_bounds(bayestopt_, options_.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
% initialization of state particles
-[~, M_, options_, oo_, ReducedForm] = solve_model_for_online_filter(true, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
+[~, ~, ReducedForm] = solve_model_for_online_filter(true, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
order = options_.order;
mf0 = ReducedForm.mf0;
@@ -80,7 +80,7 @@ for i=1:number_of_particles
info = 12042009;
while info
candidate = Prior.draw();
- [info, M_, options_, oo_] = solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
+ [info] = solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if ~info
xparam(:,i) = candidate(:);
end
@@ -117,8 +117,8 @@ for t=1:sample_size
tau_tilde = zeros(1,number_of_particles);
for i=1:number_of_particles
% model resolution
- [info, M_, options_, oo_, ReducedForm] = ...
- solve_model_for_online_filter(false, fore_xparam(:,i), dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
+ [info, M_, ReducedForm] = ...
+ solve_model_for_online_filter(false, fore_xparam(:,i), dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if ~info(1)
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
@@ -211,8 +211,8 @@ for t=1:sample_size
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
% model resolution for new parameters particles
- [info, M_, options_, oo_, ReducedForm] = ...
- solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_) ;
+ [info, M_, ReducedForm] = ...
+ solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if ~info(1)
xparam(:,i) = candidate ;
steadystate = ReducedForm.steadystate;
@@ -359,12 +359,7 @@ for t=1:sample_size
end
end
save(sprintf('%s%sparameters_particles-%u.mat', SimulationFolder, filesep(), t), 'xparam');
-
- str = sprintf(' Lower Bound (95%%) \t Mean \t\t\t Upper Bound (95%%)');
- for l=1:size(xparam,1)
- str = sprintf('%s\n %5.4f \t\t %7.5f \t\t %5.4f', str, lb95_xparam(l,t), mean_xparam(l,t), ub95_xparam(l,t));
- end
- disp(str)
+ dyntable(options_,'', {'Parameter'; 'Lower Bound (95%)';'Mean';'Upper Bound (95%)'}, bayestopt_.name, [lb95_xparam(:,t), mean_xparam(:,t), ub95_xparam(:,t)], cellofchararraymaxlength(bayestopt_.name)+2, 10, 6)
disp('')
end
diff --git a/matlab/nonlinear-filters/auxiliary_initialization.m b/matlab/nonlinear-filters/auxiliary_initialization.m
deleted file mode 100644
index 3f4f03147685089ed3d8bd84f2d65c68261784e0..0000000000000000000000000000000000000000
--- a/matlab/nonlinear-filters/auxiliary_initialization.m
+++ /dev/null
@@ -1,116 +0,0 @@
-function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,ParticleOptions,ThreadsOptions)
-
-% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling.
-
-% Copyright © 2011-2022 Dynare Team
-%
-% This file is part of Dynare (particles module).
-%
-% 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 particles module 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 <https://www.gnu.org/licenses/>.
-
-persistent init_flag mf1 number_of_particles
-persistent number_of_observed_variables number_of_structural_innovations
-
-% Set default
-if isempty(start)
- start = 1;
-end
-
-% Set flag for prunning
-%pruning = ParticleOptions.pruning;
-
-% Get steady state and mean.
-%steadystate = ReducedForm.steadystate;
-constant = ReducedForm.constant;
-ss = ReducedForm.ys;
-state_variables_steady_state = ReducedForm.state_variables_steady_state;
-
-% Set persistent variables.
-if isempty(init_flag)
- mf1 = ReducedForm.mf1;
- number_of_observed_variables = length(mf1);
- number_of_structural_innovations = length(ReducedForm.Q);
- number_of_particles = ParticleOptions.number_of_particles;
- init_flag = 1;
-end
-
-order = options_.order;
-
-% Set local state space model (first order approximation).
-ghx = ReducedForm.ghx;
-ghu = ReducedForm.ghu;
-% Set local state space model (second order approximation).
-ghxx = ReducedForm.ghxx;
-ghuu = ReducedForm.ghuu;
-ghxu = ReducedForm.ghxu;
-ghs2 = ReducedForm.ghs2;
-if (order == 3)
- ghxxx = ReducedForm.ghxxx;
- ghuuu = ReducedForm.ghuuu;
- ghxxu = ReducedForm.ghxxu;
- ghxuu = ReducedForm.ghxuu;
- ghxss = ReducedForm.ghxss;
- ghuss = ReducedForm.ghuss;
-end
-
-% Get covariance matrices
-Q = ReducedForm.Q;
-H = ReducedForm.H;
-if isempty(H)
- H = 0;
-end
-
-% Get initial condition for the state vector.
-StateVectorMean = ReducedForm.StateVectorMean;
-StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
-state_variance_rank = size(StateVectorVarianceSquareRoot,2);
-%Q_lower_triangular_cholesky = chol(Q)';
-%if pruning
-% StateVectorMean_ = StateVectorMean;
-% StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
-%end
-
-% Set seed for randn().
-options_=set_dynare_seed_local_options(options_,'default');
-
-% Initialization of the likelihood.
-const_lik = log(2*pi)*number_of_observed_variables;
-
-% Initialization of the weights across particles.
-weights = ones(1,number_of_particles)/number_of_particles ;
-StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
-%if pruning
-% StateVectors_ = StateVectors;
-%end
-yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
-%if pruning
-% 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
-if (order == 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);
-elseif (order == 3)
- tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ss, options_.threads.local_state_space_iteration_3, false);
-else
- error('Orders > 3 not allowed');
-end
-%end
-PredictedObservedMean = weights*(tmp(mf1,:)');
-PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
-dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean');
-PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + H;
-wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
-tau_tilde = weights.*wtilde ;
-tau_tilde = tau_tilde/sum(tau_tilde);
-initial_distribution = resample(StateVectors',tau_tilde',ParticleOptions)' ;
diff --git a/matlab/nonlinear-filters/sequential_importance_particle_filter.m b/matlab/nonlinear-filters/sequential_importance_particle_filter.m
index e84351686867612490f41c2dd95039cb45b8053b..0ac02e8f7d6cfe60977cff726b18dd95b9963779 100644
--- a/matlab/nonlinear-filters/sequential_importance_particle_filter.m
+++ b/matlab/nonlinear-filters/sequential_importance_particle_filter.m
@@ -19,11 +19,6 @@ function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,P
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-persistent init_flag
-persistent mf0 mf1
-persistent number_of_particles number_of_state_variables
-persistent sample_size number_of_observed_variables number_of_structural_innovations
-
% Set default value for start
if isempty(start)
start = 1;
@@ -39,17 +34,11 @@ state_variables_steady_state = ReducedForm.state_variables_steady_state;
order = options_.order;
-% Set persistent variables (if needed).
-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);
- number_of_particles = ParticleOptions.number_of_particles;
- init_flag = 1;
-end
+sample_size = size(Y,2);
+number_of_state_variables = length(ReducedForm.mf0);
+number_of_observed_variables = length(ReducedForm.mf1);
+number_of_structural_innovations = length(ReducedForm.Q);
+number_of_particles = ParticleOptions.number_of_particles;
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
@@ -89,10 +78,6 @@ lik = NaN(sample_size,1);
% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';%reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
-if pruning
- StateVectorMean_ = StateVectorMean;
- StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
-end
% Get the rank of StateVectorVarianceSquareRoot
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
@@ -110,11 +95,11 @@ if pruning
if order == 2
StateVectors_ = StateVectors;
state_variables_steady_state_ = state_variables_steady_state;
- mf0_ = mf0;
+ mf0_ = ReducedForm.mf0;
elseif order == 3
StateVectors_ = repmat(StateVectors,3,1);
state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
- mf0_ = repmat(mf0,1,3);
+ mf0_ = repmat(ReducedForm.mf0,1,3);
mask2 = number_of_state_variables+1:2*number_of_state_variables;
mask3 = 2*number_of_state_variables+1:3*number_of_state_variables;
mf0_(mask2) = mf0_(mask2)+size(ghx,1);
@@ -151,7 +136,7 @@ for t=1:sample_size
end
end
%PredictedObservedMean = tmp(mf1,:)*transpose(weights);
- PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
+ PredictionError = bsxfun(@minus,Y(:,t),tmp(ReducedForm.mf1,:));
%dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
%PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H;
%PredictedObservedVariance = H;
@@ -167,15 +152,15 @@ for t=1:sample_size
weights = wtilde/sum(wtilde);
if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
if pruning
- temp = resample([tmp(mf0,:)' tmp_(mf0_,:)'],weights',ParticleOptions);
+ temp = resample([tmp(ReducedForm.mf0,:)' tmp_(mf0_,:)'],weights',ParticleOptions);
StateVectors = temp(:,1:number_of_state_variables)';
StateVectors_ = temp(:,number_of_state_variables+1:end)';
else
- StateVectors = resample(tmp(mf0,:)',weights',ParticleOptions)';
+ StateVectors = resample(tmp(ReducedForm.mf0,:)',weights',ParticleOptions)';
end
weights = ones(1,number_of_particles)/number_of_particles;
elseif ParticleOptions.resampling.status.none
- StateVectors = tmp(mf0,:);
+ StateVectors = tmp(ReducedForm.mf0,:);
if pruning
StateVectors_ = tmp_(mf0_,:);
end
diff --git a/matlab/nonlinear-filters/solve_model_for_online_filter.m b/matlab/nonlinear-filters/solve_model_for_online_filter.m
index 6305f214e613ac8c5ee63a055e9bc9297ed12dfe..f93eb112508fa101c10894a911559ec4a9d97039 100644
--- a/matlab/nonlinear-filters/solve_model_for_online_filter.m
+++ b/matlab/nonlinear-filters/solve_model_for_online_filter.m
@@ -1,5 +1,7 @@
-function [info, M_, options_, oo_, ReducedForm] = ...
- solve_model_for_online_filter(setinitialcondition, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_)
+function [info, M_, ReducedForm] = ...
+ solve_model_for_online_filter(setinitialcondition, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
+% [info, M_, ReducedForm] = ...
+% solve_model_for_online_filter(setinitialcondition, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, dr , endo_steady_state, exo_steady_state, exo_det_steady_state)
% Solves the dsge model for an particular parameters set.
%
@@ -11,16 +13,18 @@ function [info, M_, options_, oo_, ReducedForm] = ...
% - M_ [struct] Model description.
% - estim_params_ [struct] Estimated parameters.
% - bayestopt_ [struct] Prior definition.
-% - oo_ [struct] Dynare results.
+% - bounds [struct] Prior bounds.
+% - dr [structure] Reduced form model.
+% - endo_steady_state [vector] steady state value for endogenous variables
+% - exo_steady_state [vector] steady state value for exogenous variables
+% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUTS
% - info [integer] scalar, nonzero if any problem occur when computing the reduced form.
% - M_ [struct] M_ description.
-% - options_ [struct] Dynare options.
-% - oo_ [struct] Dynare results.
% - ReducedForm [struct] Reduced form model.
-% Copyright © 2013-2023 Dynare Team
+% Copyright © 2013-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -37,83 +41,25 @@ function [info, M_, options_, oo_, ReducedForm] = ...
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-persistent init_flag restrict_variables_idx state_variables_idx mf0 mf1 number_of_state_variables
-
-info = 0;
+info = zeros(4,1);
%----------------------------------------------------
% 1. Get the structural parameters & define penalties
%----------------------------------------------------
-% Test if some parameters are smaller than the lower bound of the prior domain.
-if any(xparam1<bounds.lb)
- info = 41;
- return
-end
-
-% Test if some parameters are greater than the upper bound of the prior domain.
-if any(xparam1>bounds.ub)
- info = 42;
- return
-end
-
-% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
-Q = M_.Sigma_e;
-H = M_.H;
-for i=1:estim_params_.nvx
- k =estim_params_.var_exo(i,1);
- Q(k,k) = xparam1(i)*xparam1(i);
-end
-offset = estim_params_.nvx;
-if estim_params_.nvn
- for i=1:estim_params_.nvn
- H(i,i) = xparam1(i+offset)*xparam1(i+offset);
- end
- offset = offset+estim_params_.nvn;
-else
- H = zeros(size(dataset_.data, 2));
-end
-
-% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
-if estim_params_.ncx
- for i=1:estim_params_.ncx
- k1 =estim_params_.corrx(i,1);
- k2 =estim_params_.corrx(i,2);
- Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
- Q(k2,k1) = Q(k1,k2);
- end
- % Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
- [~, testQ] = chol(Q);
- if testQ
- % The variance-covariance matrix of the structural innovations is not definite positive.
- info = 43;
- return
- end
- offset = offset+estim_params_.ncx;
+% Ensure that xparam1 is a column vector.
+% (Don't do the transformation if xparam1 is empty, otherwise it would become a
+% 0×1 matrix, which create issues with older MATLABs when comparing with [] in
+% check_bounds_and_definiteness_estimation)
+if ~isempty(xparam1)
+ xparam1 = xparam1(:);
end
-% Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite.
-if estim_params_.ncn
- corrn_observable_correspondence = estim_params_.corrn_observable_correspondence;
- for i=1:estim_params_.ncn
- k1 = corrn_observable_correspondence(i,1);
- k2 = corrn_observable_correspondence(i,2);
- H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
- H(k2,k1) = H(k1,k2);
- end
- % Try to compute the cholesky decomposition of H (possible iff H is positive definite)
- [~, testH] = chol(H);
- if testH
- % The variance-covariance matrix of the measurement errors is not definite positive.
- info = 44;
- return
- end
- offset = offset+estim_params_.ncn;
-end
+M_ = set_all_parameters(xparam1,estim_params_,M_);
-% Update estimated structural parameters in Mode.params.
-if estim_params_.np > 0
- M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
+[~,info,~,Q,H]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, bounds);
+if info(1)
+ return
end
% Update M_.Sigma_e and M_.H.
@@ -125,33 +71,27 @@ M_.H = H;
%------------------------------------------------------------------------------
warning('off', 'MATLAB:nearlySingularMatrix')
-[oo_.dr, info, M_.params] = ...
- compute_decision_rules(M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+[dr, info, M_.params] = ...
+ compute_decision_rules(M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
warning('on', 'MATLAB:nearlySingularMatrix')
if info(1)~=0
- if nargout==5
+ if nargout==3
ReducedForm = 0;
end
return
end
-% Get decision rules and transition equations.
-dr = oo_.dr;
-
% Set persistent variables (first call).
-if isempty(init_flag)
- mf0 = bayestopt_.mf0;
- mf1 = bayestopt_.mf1;
- restrict_variables_idx = dr.restrict_var_list;
- state_variables_idx = restrict_variables_idx(mf0);
- number_of_state_variables = length(mf0);
- init_flag = true;
-end
+mf0 = bayestopt_.mf0;
+mf1 = bayestopt_.mf1;
+restrict_variables_idx = dr.restrict_var_list;
+state_variables_idx = restrict_variables_idx(mf0);
+number_of_state_variables = length(mf0);
% Return reduced form model.
-if nargout>4
+if nargout>2
ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
@@ -174,7 +114,6 @@ if nargout>4
ReducedForm.ghuu = zeros(size(restrict_variables_idx,1),n_shocks^2);
ReducedForm.ghxu = zeros(size(restrict_variables_idx,1),n_states*n_shocks);
ReducedForm.constant = ReducedForm.steadystate;
-% ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:);
end
ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
ReducedForm.Q = Q;
@@ -194,16 +133,11 @@ if setinitialcondition
StateVectorVariance = StateVectorVariance(mf0,mf0);
case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
- old_DynareOptionsperiods = options_.periods;
options_.periods = 5000;
- old_DynareOptionspruning = options_.pruning;
options_.pruning = options_.particle.pruning;
y_ = simult(dr.ys, dr, M_, options_);
y_ = y_(dr.order_var(state_variables_idx),2001:options_.periods);
StateVectorVariance = cov(y_');
- options_.periods = old_DynareOptionsperiods;
- options_.pruning = old_DynareOptionspruning;
- clear('old_DynareOptionsperiods','y_');
case 3% Initial state vector covariance is a diagonal matrix.
StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
StateVectorVariance = options_.particle.initial_state_prior_std*eye(number_of_state_variables);
diff --git a/meson.build b/meson.build
index 13e0c01d97d6c770d1b0d8657cfa3570a796b2c3..be322ea4804b91aa15c4fd8309f6a61a52493f0a 100644
--- a/meson.build
+++ b/meson.build
@@ -1772,6 +1772,7 @@ mod_and_m_tests = [
'particle/benchmark.m',
'particle/mysample.m'] },
{ 'test' : [ 'particle/first_spec.mod',
+ 'particle/first_spec_hssmc.mod',
'particle/local_state_space_iteration_k_test.mod',
'particle/local_state_space_iteration_3_test.mod' ],
'extra' : [ 'particle/first_spec_common.inc' ] },
diff --git a/tests/particle/first_spec_hssmc.mod b/tests/particle/first_spec_hssmc.mod
new file mode 100644
index 0000000000000000000000000000000000000000..0343ec263ae533175b3ec968a28e96bdfb38c2e4
--- /dev/null
+++ b/tests/particle/first_spec_hssmc.mod
@@ -0,0 +1,26 @@
+//File testing error message if initial state vector is not positive definite
+
+@#include "first_spec_common.inc"
+
+varobs q ca;
+
+shocks;
+var eeps = 0.04^2;
+var nnu = 0.03^2;
+var q = 0.01^2;
+var ca = 0.01^2;
+end;
+
+stoch_simul(order=3,periods=200, irf=0);
+send_endogenous_variables_to_workspace;
+save('my_data.mat','q','ca');
+
+estimation(datafile='my_data.mat',order=2,mode_compute=0,mh_replic=0,filter_algorithm=sis,nonlinear_filter_initialization=2
+ ,cova_compute=0, %tell program that no covariance matrix was computed
+ posterior_sampling_method='hssmc',
+ posterior_sampler_options=('steps',5,
+ 'lambda',2,
+ 'particles', 200,
+ 'scale',.5,
+ 'target', .25)
+);
\ No newline at end of file