From da159038e6130f10dd9f7a0f4b548fe56af80e46 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?=
 <stephane.adjemian@univ-lemans.fr>
Date: Fri, 21 Sep 2012 22:59:03 +0200
Subject: [PATCH] Added check for BK and steady state in the optimization of
 the prior density.

---
 matlab/get_prior_info.m             |  9 ++-
 matlab/maximize_prior_density.m     |  6 +-
 matlab/minus_logged_prior_density.m | 96 ++++++++++++++++++++++++++++-
 3 files changed, 103 insertions(+), 8 deletions(-)

diff --git a/matlab/get_prior_info.m b/matlab/get_prior_info.m
index 4498caeed..b355d4fc0 100644
--- a/matlab/get_prior_info.m
+++ b/matlab/get_prior_info.m
@@ -130,6 +130,7 @@ end
 
 if info==2% Prior optimization.
           % Initialize to the prior mode if possible
+    oo_.dr=set_state_space(oo_.dr,M_,options_);
     k = find(~isnan(bayestopt_.p5));
     xparam1(k) = bayestopt_.p5(k);
     % Pertubation of the initial condition.
@@ -139,7 +140,11 @@ if info==2% Prior optimization.
     while look_for_admissible_initial_condition
         xinit = xparam1+scale*randn(size(xparam1));
         if all(xinit(:)>bayestopt_.p3) && all(xinit(:)<bayestopt_.p4)
-            look_for_admissible_initial_condition = 0;
+            M_ = set_all_parameters(xinit,estim_params_,M_);
+            [dr,INFO,M_,options_,oo_] = resol(0,M_,options_,oo_);
+            if ~INFO(1)
+                look_for_admissible_initial_condition = 0;
+            end
         else
             if iter == 2000;
                 scale = scale/1.1;
@@ -155,7 +160,7 @@ if info==2% Prior optimization.
                                bayestopt_.p6, ...
                                bayestopt_.p7, ...
                                bayestopt_.p3, ...
-                               bayestopt_.p4);
+                               bayestopt_.p4,options_,M_,estim_params_,oo_);
     % Display the results.
     disp(' ')
     disp(' ')
diff --git a/matlab/maximize_prior_density.m b/matlab/maximize_prior_density.m
index e9a683201..cee0cc394 100644
--- a/matlab/maximize_prior_density.m
+++ b/matlab/maximize_prior_density.m
@@ -1,5 +1,5 @@
 function [xparams,lpd,hessian] = ...
-    maximize_prior_density(iparams, prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound)
+    maximize_prior_density(iparams, prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound,DynareOptions,DynareModel,EstimatedParams,DynareResults)
 % Maximizes the logged prior density using Chris Sims' optimization routine.
 % 
 % INPUTS 
@@ -41,7 +41,7 @@ gradient_method = 2;
 gradient_epsilon = 1e-6;
 
 [lpd,xparams,grad,hessian,itct,fcount,retcodehat] = ...
-    csminwel1('minus_logged_prior_density',iparams,H0,[],crit,nit,gradient_method, gradient_epsilon, ... 
-              prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound);
+    csminwel1('minus_logged_prior_density',iparams,H0,[],crit,nit,gradient_method, gradient_epsilon, ...
+              prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound,DynareOptions,DynareModel,EstimatedParams,DynareResults);
 
 lpd = -lpd;
diff --git a/matlab/minus_logged_prior_density.m b/matlab/minus_logged_prior_density.m
index 9fd6bdd67..ef883e005 100644
--- a/matlab/minus_logged_prior_density.m
+++ b/matlab/minus_logged_prior_density.m
@@ -1,4 +1,4 @@
-function [f,fake_1, fake_2, fake_3] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4)
+function [fval,fake_1, fake_2, exit_flag ] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,DynareOptions,DynareModel,EstimatedParams,DynareResults)
 % Evaluates minus the logged prior density.
 % 
 % INPUTS 
@@ -28,8 +28,98 @@ function [f,fake_1, fake_2, fake_3] = minus_logged_prior_density(xparams,pshape,
 %
 % 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
+
 fake_1 = 1;
 fake_2 = 1;
-fake_3 = 1;
 
-f = - priordens(xparams,pshape,p6,p7,p3,p4);
\ No newline at end of file
+exit_flag = 1;
+info = 0;
+
+%------------------------------------------------------------------------------
+% 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 ~isequal(DynareOptions.mode_compute,1) && any(xparams<p3)
+    k = find(xparams<p3);
+    fval = objective_function_penalty_base+sum((p3(k)-xparams(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 ~isequal(DynareOptions.mode_compute,1) && any(xparams>p4)
+    k = find(xparams>p4);
+    fval = objective_function_penalty_base+sum((xparams(k)-p4(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).
+DynareModel = set_all_parameters(xparams,EstimatedParams,DynareModel);
+
+Q = DynareModel.Sigma_e;
+H = DynareModel.H;
+
+% Test if Q is positive definite.
+if EstimatedParams.ncx
+    % 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
+end
+
+% Test if H is positive definite.
+if EstimatedParams.ncn
+    % 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
+end
+
+
+%-----------------------------
+% 2. Check BK and steady state
+%-----------------------------
+
+M_ = set_all_parameters(xparams,EstimatedParams,DynareModel);
+[dr,info,DynareModel,DynareOptions,DynareResults] = resol(0,DynareModel,DynareOptions,DynareResults);
+
+% 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 || info(1) == 19
+    fval = objective_function_penalty_base+1;
+    info = info(1);
+    exit_flag = 0;
+    return
+elseif info(1) == 3 || info(1) == 4 || info(1)==6 || info(1) == 20 || info(1) == 21  || info(1) == 23
+    fval = objective_function_penalty_base+info(2);
+    info = info(1);
+    exit_flag = 0;
+    return
+end
+
+
+fval = - priordens(xparams,pshape,p6,p7,p3,p4);
\ No newline at end of file