From ab6b642596683987bd1193696d2417dcf3041554 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Tue, 17 May 2022 12:34:26 +0200
Subject: [PATCH] NK_baseline_steadystate.m: provide proper error handling

---
 examples/NK_baseline_steadystate.m | 28 ++++++++++++++++++----------
 1 file changed, 18 insertions(+), 10 deletions(-)

diff --git a/examples/NK_baseline_steadystate.m b/examples/NK_baseline_steadystate.m
index d0f3629a9..3ea4f0327 100644
--- a/examples/NK_baseline_steadystate.m
+++ b/examples/NK_baseline_steadystate.m
@@ -2,14 +2,14 @@ function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
 % function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
 % computes the steady state for the NK_baseline.mod and uses a numerical
 % solver to do so
-% Inputs: 
+% Inputs:
 %   - ys        [vector] vector of initial values for the steady state of
 %                   the endogenous variables
 %   - exo       [vector] vector of values for the exogenous variables
 %   - M_        [structure] Dynare model structure
 %   - options   [structure] Dynare options structure
 %
-% Output: 
+% Output:
 %   - ys        [vector] vector of steady state values for the the endogenous variables
 %   - params    [vector] vector of parameter values
 %   - check     [scalar] set to 0 if steady state computation worked and to
@@ -35,8 +35,8 @@ function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
 % read out parameters to access them with their name
 NumberOfParameters = M_.param_nbr;
 for ii = 1:NumberOfParameters
-  paramname = M_.param_names{ii};
-  eval([ paramname ' = M_.params(' int2str(ii) ');']);
+    paramname = M_.param_names{ii};
+    eval([ paramname ' = M_.params(' int2str(ii) ');']);
 end
 % initialize indicator
 check = 0;
@@ -65,6 +65,7 @@ Lambdax=mu_z;
 %set the parameter gammma1
 gammma1=mu_z*mu_I/betta-(1-delta);
 if gammma1<0 % parameter violates restriction; Preventing this cannot be implemented via prior restriction as it is a composite of different parameters and the valid prior region has unknown form
+    params=M_.params;
     check=1; %set failure indicator
     return; %return without updating steady states
 end
@@ -86,13 +87,20 @@ vp=(1-thetap)/(1-thetap*PI^((1-chi)*epsilon))*PIstar^(-epsilon);
 vw=(1-thetaw)/(1-thetaw*PI^((1-chiw)*eta)*mu_z^eta)*PIstarw^(-eta);
 tempvaromega=alppha/(1-alppha)*w/r*mu_z*mu_I;
 
-[ld,fval,exitflag]=fzero(@(ld)(1-betta*thetaw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)))/(1-betta*thetaw*mu_z^(eta*(1+gammma))*PI^(eta*(1-chiw)*(1+gammma)))...
--(eta-1)/eta*wstar/(varpsi*PIstarw^(-eta*gammma)*ld^gammma)*((1-h*mu_z^(-1))^(-1)-betta*h*(mu_z-h)^(-1))*...
-((mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alppha-tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi)^(-1),0.25,options);
+try
+    %proper error handling for cases for infeasible initial value, which would result in error instead of valid exitflag
+    [ld,fval,exitflag]=fzero(@(ld)(1-betta*thetaw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)))/(1-betta*thetaw*mu_z^(eta*(1+gammma))*PI^(eta*(1-chiw)*(1+gammma)))...
+        -(eta-1)/eta*wstar/(varpsi*PIstarw^(-eta*gammma)*ld^gammma)*((1-h*mu_z^(-1))^(-1)-betta*h*(mu_z-h)^(-1))*...
+        ((mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alppha-tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi)^(-1),0.25,options);
+catch
+    exitflag = 0;
+end
+
 if exitflag <1
     %indicate the SS computation was not sucessful; this would also be detected by Dynare
     %setting the indicator here shows how to use this functionality to
     %filter out parameter draws
+    params=M_.params;
     check=1; %set failure indicator
     return; %return without updating steady states
 end
@@ -115,11 +123,11 @@ g2=epsilon/(epsilon-1)*g1;
 
 params=NaN(NumberOfParameters,1);
 for iter = 1:length(M_.params) %update parameters set in the file
-  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
+    eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
 end
 
 NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
 for ii = 1:NumberOfEndogenousVariables
-  varname = M_.endo_names{ii};
-  eval(['ys(' int2str(ii) ') = ' varname ';']);
+    varname = M_.endo_names{ii};
+    eval(['ys(' int2str(ii) ') = ' varname ';']);
 end