diff --git a/matlab/display_problematic_vars_Jacobian.m b/matlab/display_problematic_vars_Jacobian.m
index b77c60804..1abf0ba71 100644
--- a/matlab/display_problematic_vars_Jacobian.m
+++ b/matlab/display_problematic_vars_Jacobian.m
@@ -36,9 +36,13 @@ function []=display_problematic_vars_Jacobian(problemrow,problemcol,M_,x,type,ca
 if nargin<6
     caller_string='';
 end
-aux_vars_type = [M_.aux_vars.type];
-aux_eq_nbr = max(find(aux_vars_type == 6)); 
-if isempty(aux_eq_nbr)
+if ~isempty(M_.aux_vars)
+    aux_vars_type = [M_.aux_vars.type];
+    aux_eq_nbr = max(find(aux_vars_type == 6)); 
+    if isempty(aux_eq_nbr)
+        aux_eq_nbr=0;
+    end
+else
     aux_eq_nbr=0;
 end
 if strcmp(type,'dynamic')
@@ -84,6 +88,7 @@ if strcmp(type,'dynamic')
     fprintf('\n%s  The problem most often occurs, because a variable with\n',caller_string)
     fprintf('%s  exponent smaller than 1 has been initialized to 0. Taking the derivative\n',caller_string)
     fprintf('%s  and evaluating it at the steady state then results in a division by 0.\n',caller_string)
+    fprintf('%s  If you are using model-local variables (# operator), check their values as well.\n',caller_string)
 elseif strcmp(type,'static')
     for ii=1:length(problemrow)
         if problemcol(ii)<=M_.orig_endo_nbr
@@ -119,6 +124,7 @@ elseif strcmp(type,'static')
     fprintf('\n%s  The problem most often occurs, because a variable with\n',caller_string)
     fprintf('%s exponent smaller than 1 has been initialized to 0. Taking the derivative\n',caller_string)
     fprintf('%s and evaluating it at the steady state then results in a division by 0.\n',caller_string)
+    fprintf('%s  If you are using model-local variables (# operator), check their values as well.\n',caller_string)
 else
     error('Unknown Type')    
 end
\ No newline at end of file
diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m
index 21b7d7d71..83478960b 100644
--- a/matlab/model_diagnostics.m
+++ b/matlab/model_diagnostics.m
@@ -123,6 +123,12 @@ for b=1:nb
         fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the static model contains Inf or NaN. The problem arises from: \n\n')
         display_problematic_vars_Jacobian(infrow,infcol,M,dr.ys,'static','MODEL_DIAGNOSTICS: ')
     end
+    if any(any(~isreal(jacob)))
+        problem_dummy=1;
+        [imagrow,imagcol]=find(abs(imag(jacob))>1e-15);
+        fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the static model contains imaginary parts. The problem arises from: \n\n')
+        display_problematic_vars_Jacobian(imagrow,imagcol,M,dr.ys,'static','MODEL_DIAGNOSTICS: ')
+    end
     try
         rank_jacob = rank(jacob); %can sometimes fail
     catch
@@ -227,6 +233,12 @@ if any(any(isinf(jacobia_) | isnan(jacobia_)))
     fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains Inf or NaN. The problem arises from: \n\n')
     display_problematic_vars_Jacobian(infrow,infcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ')
 end
+if any(any(~isreal(jacobia_)))
+    problem_dummy=1;
+    [imagrow,imagcol]=find(abs(imag(jacobia_))>1e-15);
+    fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n')
+    display_problematic_vars_Jacobian(imagrow,imagcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ')
+end
 if exist('hessian1','var')
     if any(any(isinf(hessian1) | isnan(hessian1)))
         problem_dummy=1;
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index f3a98db70..8e4504905 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -158,6 +158,11 @@ if ~isreal(jacobia_)
     if max(max(abs(imag(jacobia_)))) < 1e-15
         jacobia_ = real(jacobia_);
     else
+        if options_.debug
+            [imagrow,imagcol]=find(abs(imag(jacobia_))>1e-15);
+            fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n')
+            display_problematic_vars_Jacobian(imagrow,imagcol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ')
+        end
         info(1) = 6;
         info(2) = sum(sum(imag(jacobia_).^2));
         return