From d75e7278c22f02030e818906e299c14ffa45233f Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Sat, 27 May 2017 19:26:12 +0200
Subject: [PATCH] Make sure that reverting to order=1 if Hessian is 0 does not
 crash other functions

Closes #1453
---
 matlab/disp_dr.m                     |  5 ++++
 matlab/disp_th_moments.m             |  2 +-
 matlab/irf.m                         | 13 +++++---
 matlab/th_autocovariances.m          |  9 ++++--
 tests/Makefile.am                    |  1 +
 tests/ls2003/ls2003_hessian_zero.mod | 44 ++++++++++++++++++++++++++++
 6 files changed, 67 insertions(+), 7 deletions(-)
 create mode 100644 tests/ls2003/ls2003_hessian_zero.mod

diff --git a/matlab/disp_dr.m b/matlab/disp_dr.m
index f8f4e4b4c..d5fcb84fb 100644
--- a/matlab/disp_dr.m
+++ b/matlab/disp_dr.m
@@ -27,6 +27,11 @@ function disp_dr(dr,order,var_list)
 
 global M_ options_
 
+if M_.hessian_eq_zero && order~=1
+    order = 1;
+    warning('disp_dr: using order = 1 because Hessian is equal to zero');
+end
+    
 nx =size(dr.ghx,2);
 nu =size(dr.ghu,2);
 if options_.block
diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index 014102f5c..cfb8f9eae 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -44,7 +44,7 @@ m(non_stationary_vars) = NaN;
 i1 = find(abs(diag(oo_.gamma_y{1})) > 1e-12);
 s2 = diag(oo_.gamma_y{1});
 sd = sqrt(s2);
-if options_.order == 2
+if options_.order == 2 && ~M_.hessian_eq_zero
     m = m+oo_.gamma_y{options_.ar+3};
 end
 
diff --git a/matlab/irf.m b/matlab/irf.m
index c80b65b1b..2355392d7 100644
--- a/matlab/irf.m
+++ b/matlab/irf.m
@@ -44,11 +44,16 @@ else
 end
 y       = 0;
 
-if iorder == 1
+local_order = iorder;
+if M_.hessian_eq_zero && local_order~=1
+    local_order = 1;
+end
+
+if local_order == 1
     y1 = repmat(dr.ys,1,long);
     ex2 = zeros(long,M_.exo_nbr);
     ex2(1,:) = e1';
-    y2 = simult_(temps,dr,ex2,iorder);
+    y2 = simult_(temps,dr,ex2,local_order);
     y = y2(:,M_.maximum_lag+1:end)-y1;
 else
     % eliminate shocks with 0 variance
@@ -61,8 +66,8 @@ else
         ex1(:,i_exo_var) = randn(long+drop,nxs)*chol_S;
         ex2 = ex1;
         ex2(drop+1,:) = ex2(drop+1,:)+e1';
-        y1 = simult_(temps,dr,ex1,iorder);
-        y2 = simult_(temps,dr,ex2,iorder);
+        y1 = simult_(temps,dr,ex1,local_order);
+        y2 = simult_(temps,dr,ex2,local_order);
         y = y+(y2(:,M_.maximum_lag+drop+1:end)-y1(:,M_.maximum_lag+drop+1:end));
     end
     y=y/replic;
diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m
index 78e6fff4b..3cde22f31 100644
--- a/matlab/th_autocovariances.m
+++ b/matlab/th_autocovariances.m
@@ -67,6 +67,11 @@ if options_.order >= 3
     error('Theoretical moments not implemented above 2nd order')
 end
 
+local_order = options_.order;
+if M_.hessian_eq_zero && local_order~=1
+    local_order = 1;
+end
+
 endo_nbr = M_.endo_nbr;
 exo_names_orig_ord  = M_.exo_names_orig_ord;
 if isoctave
@@ -128,7 +133,7 @@ end
 % Compute stationary variables (before HP filtering),
 % and compute 2nd order mean correction on stationary variables (in case of
 % HP filtering, this mean correction is computed *before* filtering)
-if options_.order == 2 || options_.hp_filter == 0
+if local_order == 2 || options_.hp_filter == 0
     [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
     if options_.block == 0
         iky = inv_order_var(ivar);
@@ -143,7 +148,7 @@ if options_.order == 2 || options_.hp_filter == 0
     end
     aa = ghx(iky,:);
     bb = ghu(iky,:);
-    if options_.order == 2         % mean correction for 2nd order
+    if local_order == 2         % mean correction for 2nd order
         if ~isempty(ikx)
             Ex = (dr.ghs2(ikx)+dr.ghxx(ikx,:)*vx(:)+dr.ghuu(ikx,:)*M_.Sigma_e(:))/2;
             Ex = (eye(n0)-AS(ikx,:))\Ex;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index dbc4f69c7..3775f822d 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -28,6 +28,7 @@ MODFILES = \
 	estimation/fs2000_MCMC_jumping_covariance.mod \
 	ms-sbvar/test_ms_variances_repeated_runs.mod \
 	fs2000/fs2000.mod \
+	ls2003/ls2003_hessian_zero.mod \
 	ep/rbc.mod \
 	estimation/fs2000_with_weibull_prior.mod \
 	estimation/fs2000_initialize_from_calib.mod	\
diff --git a/tests/ls2003/ls2003_hessian_zero.mod b/tests/ls2003/ls2003_hessian_zero.mod
new file mode 100644
index 000000000..7476b5b3d
--- /dev/null
+++ b/tests/ls2003/ls2003_hessian_zero.mod
@@ -0,0 +1,44 @@
+//test whether Dynare correctly reverts to linear approximation if 0 Hessian is detected
+
+var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
+varexo e_R e_q e_ys e_pies e_A;
+
+parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
+
+psi1 = 1.54;
+psi2 = 0.25;
+psi3 = 0.25;
+rho_R = 0.5;
+alpha = 0.3;
+rr = 2.51;
+k = 0.5;
+tau = 0.5;
+rho_q = 0.4;
+rho_A = 0.2;
+rho_ys = 0.9;
+rho_pies = 0.7;
+
+
+model;
+y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
+pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
+pie = de+(1-alpha)*dq+pie_s;
+R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
+dq = rho_q*dq(-1)+e_q;
+y_s = rho_ys*y_s(-1)+e_ys;
+pie_s = rho_pies*pie_s(-1)+e_pies;
+A = rho_A*A(-1)+e_A;
+y_obs = y-y(-1)+A;
+pie_obs = 4*pie;
+R_obs = 4*R;
+end;
+
+shocks;
+var e_R = 1.25^2;
+var e_q = 2.5^2;
+var e_A = 1.89;
+var e_ys = 1.89;
+var e_pies = 1.89;
+end;
+
+stoch_simul(order=2);
\ No newline at end of file