fixed problem with penalty in estimation. Created a new global scalar:

objective_function_penalty_base. It is the only simple way that I
found to keep csminwel1.m to be able to handle general functions.

matlab/DsgeVarLikelihood.m

@@ -34,11 +34,10 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihoo
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+global objective_function_penalty_base
 % Declaration of the persistent variables.
 persistent dsge_prior_weight_idx
 
-penalty = BayesInfo.penalty;
-
 grad=[];
 hess=[];
 exit_flag = [];
@@ -85,7 +84,7 @@ exit_flag = 1;
 % Return, with endogenous penalty, if some dsge-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 = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
+    fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
     exit_flag = 0;
     info = 41;
     return;
@@ -94,7 +93,7 @@ end
 % Return, with endogenous penalty, if some dsge-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 = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
+    fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
     exit_flag = 0;
     info = 42;
     return;
@@ -131,7 +130,7 @@ dsge_prior_weight = Model.params(dsge_prior_weight_idx);
 
 % Is the dsge prior proper?
 if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.info.ntobs;
-    fval = penalty+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
+    fval = objective_function_penalty_base+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
     exit_flag = 0;
     info = 51;
     return
@@ -148,12 +147,12 @@ end
 % Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
 if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
             == 8 || info(1) == 22 || info(1) == 24
-    fval = penalty+1;
+    fval = objective_function_penalty_base+1;
     info = info(1);
     exit_flag = 0;
     return
 elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21
-    fval = penalty+info(2);
+    fval = objective_function_penalty_base+info(2);
     info = info(1);
     exit_flag = 0;
     return
@@ -228,7 +227,7 @@ if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model
     if ~ispd(SIGMAu)
         v = diag(SIGMAu);
         k = find(v<0);
-        fval = penalty + sum(v(k).^2);
+        fval = objective_function_penalty_base + sum(v(k).^2);
         info = 52;
         exit_flag = 0;
         return;

matlab/csminwel1.m

@@ -40,7 +40,7 @@ function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,m
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-% global bayestopt_
+global objective_function_penalty_base
 
 fh = [];
 xh = [];
@@ -91,7 +91,9 @@ f=f0;
 H=H0;
 cliff=0;
 while ~done
-    varargin{5}.penalty = f;
+    % penalty for dsge_likelihood and DsgeVarLikelihood
+    objective_function_penalty_base = f;
+    
     g1=[]; g2=[]; g3=[];
     %addition fj. 7/6/94 for control
     disp('-----------------')

matlab/dsge_likelihood.m

@@ -130,7 +130,7 @@ function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,Bayes
 
 % AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR
 
-penalty = BayesInfo.penalty;
+global objective_function_penalty_base
 
 % Initialization of the returned variables and others...
 fval        = [];
@@ -177,7 +177,7 @@ end
 % Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
 if ~isequal(DynareOptions.mode_compute,1) && any(xparam1<BayesInfo.lb)
     k = find(xparam1<BayesInfo.lb);
-    fval = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
+    fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
     exit_flag = 0;
     info = 41;
     if analytic_derivation,
@@ -189,7 +189,7 @@ end
 % Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
 if ~isequal(DynareOptions.mode_compute,1) && any(xparam1>BayesInfo.ub)
     k = find(xparam1>BayesInfo.ub);
-    fval = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
+    fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
     exit_flag = 0;
     info = 42;
     if analytic_derivation,
@@ -213,7 +213,7 @@ if EstimatedParameters.ncx
         a = diag(eig(Q));
         k = find(a < 0);
         if k > 0
-            fval = penalty+sum(-a(k));
+            fval = objective_function_penalty_base+sum(-a(k));
             exit_flag = 0;
             info = 43;
             return
@@ -230,7 +230,7 @@ if EstimatedParameters.ncn
         a = diag(eig(H));
         k = find(a < 0);
         if k > 0
-            fval = penalty+sum(-a(k));
+            fval = objective_function_penalty_base+sum(-a(k));
             exit_flag = 0;
             info = 44;
             return
@@ -249,7 +249,7 @@ end
 % Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
 if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
             == 8 || info(1) == 22 || info(1) == 24
-    fval = penalty+1;
+    fval = objective_function_penalty_base+1;
     info = info(1);
     exit_flag = 0;
     if analytic_derivation,
@@ -257,7 +257,7 @@ if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
     end
     return
 elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21  || info(1) == 23
-    fval = penalty+info(2);
+    fval = objective_function_penalty_base+info(2);
     info = info(1);
     exit_flag = 0;
     if analytic_derivation,

matlab/dynare_estimation_init.m

@@ -35,6 +35,8 @@ function [dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_, fake] =
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+global objective_function_penalty_base
+
 if isempty(gsa_flag)
     gsa_flag = 0;
 else% Decide if a DSGE or DSGE-VAR has to be estimated.
@@ -191,7 +193,7 @@ else% Yes!
 end
 
 % Set the "size" of penalty.
-bayestopt_.penalty = 1e8;
+objective_function_penalty_base = 1e8;
 
 % Get informations about the variables of the model.
 dr = set_state_space(oo_.dr,M_,options_);

matlab/independent_metropolis_hastings_core.m

@@ -38,7 +38,7 @@ if nargin<4,
 end
 
 
-global bayestopt_ estim_params_ options_  M_ oo_
+global bayestopt_ estim_params_ options_  M_ oo_ objective_function_penalty_base
 
 % Reshape 'myinputs' for local computation.
 % In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
@@ -69,7 +69,7 @@ if whoiam
 end
 
 % (re)Set the penalty.
-bayestopt_.penalty = Inf;
+objective_function_penalty_base = Inf;
 
 MhDirectoryName = CheckPath('metropolis',M_.dname);
 

matlab/newrat.m

@@ -40,6 +40,8 @@ function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ft
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+global objective_function_penalty_base
+
 icount=0;
 nx=length(x);
 xparam1=x;
@@ -106,7 +108,7 @@ while norm(gg)>gtol && check==0 && jit<nit
     jit=jit+1;
     tic
     icount=icount+1;
-    BayesInfo.penalty = fval0(icount);
+    objective_function_penalty_base = fval0(icount);
     disp([' '])
     disp(['Iteration ',num2str(icount)])
     [fval,x0,fc,retcode] = csminit1(func0,xparam1,fval0(icount),gg,0,H,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);

matlab/non_linear_dsge_likelihood.m

@@ -121,22 +121,12 @@ function [fval,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,Dynar
 % 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 penalty
 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 persistent variable.
-if ~nargin || isempty(penalty)
-    penalty = 1e8;
-    if ~nargin, return, end
-end
-if nargin==1
-    penalty = xparam1;
-    return
-end
-
 % Initialization of the returned arguments.
 fval            = [];
 ys              = [];
@@ -153,7 +143,7 @@ nvobs = DynareDataset.info.nvobs;
 % 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 = penalty+sum((BayesInfo.lb(k)-xparam1(k)).^2);
+    fval = objective_function_penalty_base+sum((BayesInfo.lb(k)-xparam1(k)).^2);
     exit_flag = 0;
     info = 41;
     return
@@ -162,7 +152,7 @@ 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 = penalty+sum((xparam1(k)-BayesInfo.ub(k)).^2);
+    fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
     exit_flag = 0;
     info = 42;
     return
@@ -201,7 +191,7 @@ if EstimatedParameters.ncx
         a = diag(eig(Q));
         k = find(a < 0);
         if k > 0
-            fval = penalty+sum(-a(k));
+            fval = objective_function_penalty_base+sum(-a(k));
             exit_flag = 0;
             info = 43;
             return
@@ -225,7 +215,7 @@ if EstimatedParameters.ncn
         a = diag(eig(H));
         k = find(a < 0);
         if k > 0
-            fval = penalty+sum(-a(k));
+            fval = objective_function_penalty_base+sum(-a(k));
             exit_flag = 0;
             info = 44;
             return
@@ -251,11 +241,11 @@ Model.H = H;
 [T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
 
 if info(1) == 1 || info(1) == 2 || info(1) == 5
-    fval = penalty+1;
+    fval = objective_function_penalty_base+1;
     exit_flag = 0;
     return
 elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21
-    fval = penalty+info(2);
+    fval = objective_function_penalty_base+info(2);
     exit_flag = 0;
     return
 end
@@ -357,10 +347,10 @@ ReducedForm.StateVectorVariance = StateVectorVariance;
 DynareOptions.warning_for_steadystate = 0;
 LIK = feval(DynareOptions.particle.algorithm,ReducedForm,Y,[],DynareOptions);
 if imag(LIK)
-    likelihood = penalty;
+    likelihood = objective_function_penalty_base;
     exit_flag  = 0;
 elseif isnan(LIK)
-    likelihood = penalty;
+    likelihood = objective_function_penalty_base;
     exit_flag  = 0;
 else
     likelihood = LIK;

matlab/random_walk_metropolis_hastings.m

@@ -60,7 +60,8 @@ function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv
 % In Metropolis, we set penalty to Inf to as to reject all parameter sets
 % triggering error in target density computation
 
-bayestopt_.penalty = Inf;
+global objective_function_penalty_base
+objective_function_penalty_base = Inf;
 
 %%%%
 %%%% Initialization of the random walk metropolis-hastings chains.