Skip to content

GitLab

Commit 3c69ee04 authored by Frédéric Karamé's avatar Frédéric Karamé
Browse files

new file that returns only the index for resampling particles

parent 22a29ab0
Hide whitespace changes
Inline Side-by-side
matlab/particle/index_resample.m 0 → 100644
View file @ 3c69ee04
function resampling_index = index_resample(particles,weights,DynareOptions)
% Resamples particles.
%@info:
%! @deftypefn {Function File} {@var{indx} =} resample (@var{weights}, @var{method})
%! @anchor{particle/resample}
%! @sp 1
%! Resamples particles.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item weights
%! n*1 vector of doubles, particles' weights.
%! @item method
%! string equal to 'residual' or 'traditional'.
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item indx
%! n*1 vector of intergers, indices.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{particle/sequantial_importance_particle_filter}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{residual_resampling}, @ref{traditional_resampling}
%! @sp 2
%! @end deftypefn
%@eod:
% Copyright (C) 2011, 2012 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% 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')
resampling_index = residual_resampling(particles,weights,rand);
elseif strcmpi(DynareOptions.particle.resampling.method2,'stratified')
resampling_index = residual_resampling(particles,weights,rand(size(weights)));
else
error('particle::resample: Unknown method!')
end
case 'traditional'
if strcmpi(DynareOptions.particle.resampling.method2,'kitagawa')
resampling_index = traditional_resampling(particles,weights,rand);
elseif strcmpi(DynareOptions.particle.resampling.method2,'stratified')
resampling_index = traditional_resampling(particles,weights,rand(size(weights)));
else
error('particle::resample: Unknown method!')
end
otherwise
error('particle::resample: Unknown method!')
end
\ No newline at end of file
matlab/particle/online_auxiliary_filter.m 0 → 100644
View file @ 3c69ee04
function [xparam,std_param,lb_95,ub_95,median_param] = online_auxiliary_filter(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% Carvalho & Lopes particle filter = auxiliary particle filter including Liu & West filter on parameters.
%
% INPUTS
% ReducedForm [structure] Matlab's structure describing the reduced form model.
% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors.
% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors.
% ReducedForm.state.dr [structure] output of resol.m.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% mf [integer] pp*1 vector of indices.
% number_of_particles [integer] scalar.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2013 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent Y init_flag mf0 mf1 bounds number_of_particles number_of_parameters liu_west_delta liu_west_chol_sigma_bar
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;
% 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) ;
% Set persistent variables.
if isempty(init_flag)
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_particles = DynareOptions.particle.number_of_particles;
number_of_parameters = size(xparam1,1) ;
Y = DynareDataset.data ;
sample_size = size(Y,2);
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 ;
bounds = [BayesInfo.lb BayesInfo.ub] ;
init_flag = 1;
end
% Get initial conditions for the state particles
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
StateVectors_ = StateVectors;
end
% Initialization of parameter particles
h_square = (3*liu_west_delta-1)/(2*liu_west_delta) ;
h_square = 1-h_square*h_square ;
small_a = sqrt(1-h_square) ;
xparam = zeros(number_of_parameters,number_of_particles) ;
%candidate = prior_draw(1) ;
i = 1 ;
while i<=number_of_particles
candidate = start_param + liu_west_chol_sigma_bar*rand(number_of_parameters,1) ;
%candidate = prior_draw(0) ;
if all(candidate(:) >= bounds(:,1)) && all(candidate(:) <= bounds(:,2))
xparam(:,i) = candidate(:) ;
i = i+1 ;
end
end
% Initialization of the weights of particles.
weights = ones(1,number_of_particles)/number_of_particles ;
% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables;
mean_xparam = zeros(number_of_parameters,sample_size) ;
median_xparam = zeros(number_of_parameters,sample_size) ;
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
for t=1:sample_size
disp(t)
% Moments of parameters particles distribution
m_bar = xparam*(weights') ;
temp = bsxfun(@minus,xparam,m_bar) ;
sigma_bar = (bsxfun(@times,weights,temp))*(temp') ;
chol_sigma_bar = chol(h_square*sigma_bar)' ;
% Prediction (without shocks)
ObservedVariables = zeros(number_of_observed_variables,number_of_particles) ;
for i=1:number_of_particles
% 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;
% Set local state space model (second-order approximation).
constant = ReducedForm.constant;
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
% particle likelihood contribution
yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
if pruning
yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
[tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,1),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
else
tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,1),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
end
ObservedVariables(:,i) = tmp(mf1,:) ;
end
PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ;
PredictionError = bsxfun(@minus,Y(:,t),ObservedVariables);
dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean);
PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ;
wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
% unormalized weights and observation likelihood contribution
tau_tilde = weights.*wtilde ;
sum_tau_tilde = sum(tau_tilde) ;
% particles selection
tau_tilde = tau_tilde/sum_tau_tilde ;
indx = index_resample(0,tau_tilde',DynareOptions);
StateVectors = StateVectors(:,indx) ;
if pruning
StateVectors_ = StateVectors_(:,indx) ;
end
xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ;
xparam = xparam(:,indx) ;
wtilde = wtilde(indx) ;
% draw in the new distributions
i = 1 ;
while i<=number_of_particles
candidate = xparam(:,i) + chol_sigma_bar*rand(number_of_parameters,1) ;
if all(candidate >= bounds(:,1)) && all(candidate <= bounds(:,2))
xparam(:,i) = candidate ;
% model resolution for new parameters particles
[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;
% Set local state space model (second order approximation).
constant = ReducedForm.constant;
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
ghxx = ReducedForm.ghxx;
ghuu = ReducedForm.ghuu;
ghxu = ReducedForm.ghxu;
% Get covariance matrices and structural shocks
epsilon = chol(ReducedForm.Q)'*rand(number_of_structural_innovations,1) ;
% compute particles likelihood contribution
yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
if pruning
yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
[tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
StateVectors_(:,i) = tmp_(mf0,:) ;
else
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
end
StateVectors(:,i) = tmp(mf0,:) ;
ObservedVariables(:,i) = tmp(mf1,:) ;
i = i+1 ;
end
end
PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ;
PredictionError = bsxfun(@minus,Y(:,t),ObservedVariables);
dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean);
PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ;
lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1)));
% importance ratio
wtilde = lnw./wtilde ;
% normalization
weights = wtilde/sum(wtilde);
% final resampling (advised)
indx = index_resample(0,weights,DynareOptions);
StateVectors = StateVectors(:,indx) ;
if pruning
StateVectors_ = StateVectors_(:,indx) ;
end
xparam = xparam(:,indx) ;
weights = ones(1,number_of_particles)/number_of_particles ;
%
mean_xparam(:,t) = mean(xparam,2);
mat_var_cov = bsxfun(@minus,xparam,mean_xparam(:,t)) ;
mat_var_cov = (mat_var_cov*mat_var_cov')/(number_of_particles-1) ;
std_xparam(:,t) = sqrt(diag(mat_var_cov)) ;
for i=1:number_of_parameters
temp = sortrows(xparam(i,:)') ;
lb95_xparam(i,t) = temp(0.025*number_of_particles) ;
median_xparam(i,t) = temp(0.5*number_of_particles) ;
ub95_xparam(i,t) = temp(0.975*number_of_particles) ;
end
end
xparam = mean_xparam(:,sample_size) ;
std_param = std_xparam(:,sample_size) ;
lb_95 = lb95_xparam(:,sample_size) ;
ub_95 = ub95_xparam(:,sample_size) ;
median_param = median_xparam(:,sample_size) ;
%% Plot parameters trajectory
TeX = DynareOptions.TeX;
[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters);
if TeX
fidTeX = fopen([Model.fname '_param_traj.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
z = 1:1:sample_size ;
for plt = 1:nbplt,
if TeX
NAMES = [];
TeXNAMES = [];
end
hh = dyn_figure(DynareOptions,'Name','Parameters Trajectories');
for k=1:min(nstar,length(xparam)-(plt-1)*nstar)
subplot(nr,nc,k)
kk = (plt-1)*nstar+k;
[name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions);
if TeX
if isempty(NAMES)
NAMES = name;
TeXNAMES = texname;
else
NAMES = char(NAMES,name);
TeXNAMES = char(TeXNAMES,texname);
end
end
y = [mean_xparam(kk,:)' median_xparam(kk,:)' lb95_xparam(kk,:)' ub95_xparam(kk,:)' xparam(kk)*ones(sample_size,1)] ;
plot(z,y);
hold on
title(name,'interpreter','none')
hold off
axis tight
drawnow
end
dyn_saveas(hh,[ Model.fname '_param_traj' int2str(plt) ],DynareOptions);
if TeX
% TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:min(nstar,length(x)-(plt-1)*nstar)
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParamTraj%s}\n',Model.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Parameters trajectories.}');
fprintf(fidTeX,'\\label{Fig:ParametersPlots:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
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.
for plt = 1:nbplt,
if TeX
NAMES = [];
TeXNAMES = [];
end
hh = dyn_figure(DynareOptions,'Name','Parameters Densities');
for k=1:min(nstar,length(xparam)-(plt-1)*nstar)
subplot(nr,nc,k)
kk = (plt-1)*nstar+k;
[name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions);
if TeX
if isempty(NAMES)
NAMES = name;
TeXNAMES = texname;
else
NAMES = char(NAMES,name);
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,...
number_of_particles,optimal_bandwidth,kernel_function);
plot(density(:,1),density(:,2));
hold on
title(name,'interpreter','none')
hold off
axis tight
drawnow
end
dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions);
if TeX
% TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:min(nstar,length(x)-(plt-1)*nstar)
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParametersDensities%s}\n',Model.fname,int2str(plt));
fprintf(fidTeX,'\\caption{ParametersDensities.}');
fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
end
matlab/particle/solve_model_for_online_filter.m 0 → 100644
View file @ 3c69ee04
function [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = solve_model_for_online_filter(observation_number,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% solve the dsge model for an particular parameters set.
%@info:
%! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults}] =} non_linear_dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults})
%! @anchor{dsge_likelihood}
%! @sp 1
%! Evaluates the posterior kernel of a dsge model using a non linear filter.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item xparam1
%! Vector of doubles, current values for the estimated parameters.
%! @item DynareDataset
%! Matlab's structure describing the dataset (initialized by dynare, see @ref{dataset_}).
%! @item DynareOptions
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%! @item Model
%! Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
%! @item EstimatedParamemeters
%! Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
%! @item BayesInfo
%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
%! @item DynareResults
%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item fval
%! Double scalar, value of (minus) the likelihood.
%! @item exit_flag
%! Integer scalar, equal to zero if the routine return with a penalty (one otherwise).
%! @item ys
%! Vector of doubles, steady state level for the endogenous variables.
%! @item trend_coeffs
%! Matrix of doubles, coefficients of the deterministic trend in the measurement equation.
%! @item info
%! Integer scalar, error code.
%! @table @ @code
%! @item info==0
%! No error.
%! @item info==1
%! The model doesn't determine the current variables uniquely.
%! @item info==2
%! MJDGGES returned an error code.
%! @item info==3
%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
%! @item info==4
%! Blanchard & Kahn conditions are not satisfied: indeterminacy.
%! @item info==5
%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
%! @item info==6
%! The jacobian evaluated at the deterministic steady state is complex.
%! @item info==19
%! The steadystate routine thrown an exception (inconsistent deep parameters).
%! @item info==20
%! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations).
%! @item info==21
%! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state.
%! @item info==22
%! The steady has NaNs.
%! @item info==23
%! M_.params has been updated in the steadystate routine and has complex valued scalars.
%! @item info==24
%! M_.params has been updated in the steadystate routine and has some NaNs.
%! @item info==30
%! Ergodic variance can't be computed.
%! @item info==41
%! At least one parameter is violating a lower bound condition.
%! @item info==42
%! At least one parameter is violating an upper bound condition.
%! @item info==43
%! The covariance matrix of the structural innovations is not positive definite.
%! @item info==44
%! The covariance matrix of the measurement errors is not positive definite.
%! @item info==45
%! Likelihood is not a number (NaN).
%! @item info==45
%! Likelihood is a complex valued number.
%! @end table
%! @item Model
%! Matlab's structure describing the model (initialized by dynare, see @ref{M_}).
%! @item DynareOptions
%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%! @item BayesInfo
%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
%! @item DynareResults
%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{dynare_estimation_1}, @ref{mode_check}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{priordens}
%! @end deftypefn
%@eod:
% Copyright (C) 2010-2012 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
% frederic DOT karame AT univ DASH lemans DOT fr
%global objective_function_penalty_base
% Declaration of the penalty as a persistent variable.
persistent init_flag
persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1
persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
% Initialization of the returned arguments.
fval = [];
ys = [];
trend_coeff = [];
exit_flag = 1;
% Set the number of observed variables
nvobs = DynareDataset.info.nvobs;
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
%if (DynareOptions.mode_compute~=1) && any(xparam1<BayesInfo.lb)
% k = find(xparam1(:) < BayesInfo.lb);
% fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
% exit_flag = 0;
% info = 41;
% return
%end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
%if (DynareOptions.mode_compute~=1) && any(xparam1>BayesInfo.ub)
% k = find(xparam1(:)>BayesInfo.ub);
% fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
% exit_flag = 0;
% info = 42;
% return
%end
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
Q = Model.Sigma_e;
H = Model.H;
for i=1:EstimatedParameters.nvx<