From 268095276c6fe1266577e6e2e75ba028c12203b4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?=
 <stephane.adjemian@univ-lemans.fr>
Date: Thu, 20 Jun 2013 12:59:01 +0200
Subject: [PATCH] Factorized code using ispd routine.

---
 matlab/DsgeVarLikelihood.m |  7 +++----
 matlab/dsge_likelihood.m   | 40 +++++++++++++-------------------------
 2 files changed, 17 insertions(+), 30 deletions(-)

diff --git a/matlab/DsgeVarLikelihood.m b/matlab/DsgeVarLikelihood.m
index 87d54ed52..9b57abcd2 100644
--- a/matlab/DsgeVarLikelihood.m
+++ b/matlab/DsgeVarLikelihood.m
@@ -224,10 +224,9 @@ if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model
     tmp1 = dsge_prior_weight*DynareDataset.info.ntobs*GYX + mYX;
     tmp2 = inv(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX);
     SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0');
-    if ~ispd(SIGMAu)
-        v = diag(SIGMAu);
-        k = find(v<0);
-        fval = objective_function_penalty_base + sum(v(k).^2);
+    [SIGMAu_is_positive_definite, penalty] = ispd(SIGMAu)
+    if ~SIGMAu_is_positive_definite
+        fval = objective_function_penalty_base + penalty;
         info = 52;
         exit_flag = 0;
         return;
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 225993af0..5911ced69 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -206,36 +206,24 @@ Q = Model.Sigma_e;
 H = Model.H;
 
 % Test if Q is positive definite.
-if EstimatedParameters.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
+if ~issquare(Q) && EstimatedParameters.ncx
+    [Q_is_positive_definite, penalty] = ispd(Q);
+    if ~Q_is_positive_definite
+        fval = objective_function_penalty_base+penalty;
+        exit_flag = 0;
+        info = 43;
+        return
     end
 end
 
 % Test if H is positive definite.
-if EstimatedParameters.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
+if ~issquare(H) && EstimatedParameters.ncn
+    [H_is_positive_definite, penalty] = ispd(H);
+    if ~H_is_positive_definite
+        fval = objective_function_penalty_base+penalty;
+        exit_flag = 0;
+        info = 44;
+        return
     end
 end