From f28a7a62bbc45e278a8d1624687f2c04b63d7ebe Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Thu, 10 Dec 2015 15:13:35 +0100
Subject: [PATCH] Only test non-zero measurement error covariance entries for
 positive definiteness

Otherwise, not having measurement error on one variable is not allowed during estimation
---
 matlab/dsge_likelihood.m            |  4 ++--
 matlab/dynare_estimation_init.m     | 25 +++++++++++++++++++++++++
 matlab/non_linear_dsge_likelihood.m |  4 ++--
 3 files changed, 29 insertions(+), 4 deletions(-)

diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index c07c62881..ade4a0d0d 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -210,7 +210,7 @@ H = Model.H;
 
 % Test if Q is positive definite.
 if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
-    [Q_is_positive_definite, penalty] = ispd(Q);
+    [Q_is_positive_definite, penalty] = ispd(Q(EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness,EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness));
     if ~Q_is_positive_definite
         fval = objective_function_penalty_base+penalty;
         exit_flag = 0;
@@ -231,7 +231,7 @@ end
 
 % Test if H is positive definite.
 if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME')
-    [H_is_positive_definite, penalty] = ispd(H);
+    [H_is_positive_definite, penalty] = ispd(H(EstimatedParameters.H_entries_to_check_for_positive_definiteness,EstimatedParameters.H_entries_to_check_for_positive_definiteness));
     if ~H_is_positive_definite
         fval = objective_function_penalty_base+penalty;
         exit_flag = 0;
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index 81e0b7e04..7e6d90000 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -571,3 +571,28 @@ else
         error('The option "prefilter" is inconsistent with the non-zero mean measurement equations in the model.')
     end
 end
+
+%% get the non-zero rows and columns of Sigma_e and H
+
+H_non_zero_rows=find(~all(M_.H==0,1));
+H_non_zero_columns=find(~all(M_.H==0,2));
+if ~isequal(H_non_zero_rows,H_non_zero_columns')
+    error('Measurement error matrix not symmetric')
+end
+if isfield(estim_params_,'nvn_observable_correspondence')
+    estim_params_.H_entries_to_check_for_positive_definiteness=union(H_non_zero_rows,estim_params_.nvn_observable_correspondence(:,1));
+else
+    estim_params_.H_entries_to_check_for_positive_definiteness=H_non_zero_rows;
+end
+
+Sigma_e_non_zero_rows=find(~all(M_.Sigma_e==0,1));
+Sigma_e_non_zero_columns=find(~all(M_.Sigma_e==0,2));
+if ~isequal(Sigma_e_non_zero_rows,Sigma_e_non_zero_columns')
+    error('Structual error matrix not symmetric')
+end
+if ~isempty(estim_params_.var_exo)
+    estim_params_.Sigma_e_entries_to_check_for_positive_definiteness=union(Sigma_e_non_zero_rows,estim_params_.var_exo(:,1));
+else
+    estim_params_.Sigma_e_entries_to_check_for_positive_definiteness=Sigma_e_non_zero_rows;
+end
+
diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m
index 66cc31570..22e421824 100644
--- a/matlab/non_linear_dsge_likelihood.m
+++ b/matlab/non_linear_dsge_likelihood.m
@@ -166,7 +166,7 @@ Q = Model.Sigma_e;
 H = Model.H;
 
 if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
-    [Q_is_positive_definite, penalty] = ispd(Q);
+    [Q_is_positive_definite, penalty] = ispd(Q(EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness,EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness));
     if ~Q_is_positive_definite
         fval = objective_function_penalty_base+penalty;
         exit_flag = 0;
@@ -187,7 +187,7 @@ if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calib
 end
 
 if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME')
-    [H_is_positive_definite, penalty] = ispd(H);
+    [H_is_positive_definite, penalty] = ispd(H(EstimatedParameters.H_entries_to_check_for_positive_definiteness,EstimatedParameters.H_entries_to_check_for_positive_definiteness));
     if ~H_is_positive_definite
         fval = objective_function_penalty_base+penalty;
         exit_flag = 0;