function [ys,params,check] = aob_policy_regime_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = aob_policy_regime_steadystate(ys,exo,M_,options_)
% computes the steady state for the aob_policy_regime_steadystate.mod and uses a numerical
% solver to do so
% 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:
%   - 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
%                    1 of not (allows to impose restrictions on parameters)

% Copyright (C) 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

% 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) ');']);
end
% initialize indicator
check = 0;

%% Enter model equations here

% Inflation (1)
PI=PIbar;
% capital utilization (2)
iota=1;
% Tobin's Marginal q (3)
tq=1;
% preference shifter (4)
d=0;
% UB shifter (5)
eta=0;
% SS vacancy filling rate Q (6) % fixed following Christiano et.al.(2016)
Q=0.7;
% SS employment rate l (7) % fixed following Christiano et.al.(2016)
l=0.945;

% LambdaYd (p4)
LambdaYd=alppha*LambdaPsi/(1-alppha)+LambdaZ;
% SS growth til_Phi (8)
til_Phi=exp(LambdaYd);
% SS growth til_Psi (9)
til_Psi=exp(LambdaPsi);
% SS growth til_Z (10)
til_Z=exp(LambdaZ);



%set the parameter LambdaI (p3)
LambdaI=exp(LambdaYd);

%set the parameter gamma1 (p1)
gamma1=til_Phi*til_Psi/betta-(1-delk);

% ren (11)
ren=gamma1+gamma2*(iota-1);
% R (12)
R=PI*exp(LambdaYd)/betta;

%set Rbar (p2)
Rbar=R;

% PIstar (13)
PIstar=((1-thetap*(PI)^(epsilon-1))/(1-thetap))^(1/(1-epsilon));
% Xi (14)
Xi=(1-thetap)*(1-thetap*(PI)^epsilon)^(-1)*PIstar^(-epsilon);
% mc (15)
mc=(epsilon-1)/epsilon*(1-betta*thetap*PI^epsilon)/(1-betta*thetap*PI^(epsilon-1))*PIstar;

% set mc_ss
mc_ss=mc;

% nu (16)
nu=(1-alppha)*(mc*(alppha/ren)^alppha)^(1/(1-alppha));

% set nu_ss
nu_ss=nu;

% nup (17)
nup=(1-(1-rhos)*betta)^(-1)*nu;
% mat (18)
mat=rhos*l;
% ump (19)
ump=1-(1-rhos)*l;
% job finding rate f (20)
f=mat/ump;

% set f_ss
f_ss=f;

% va (21)
va=mat/Q;

% set varatio_ss
varatio_ss=va/l;

% sigmm (p9)
sigmm=mat*ump^(-sigmaa)*va^(sigmaa-1);

% for zeta1-4
delta  = 0.0019;
M  = 60;

zeta1=1-delta+(1-delta)^(M);
zeta2=1-(1-delta)^(M);
zeta3=zeta2*(1-delta)/delta-zeta1;
zeta4=(1-delta)*zeta2/((2-delta)*M)+1-zeta2;

% gammac (p10)
gammac=0.59*nu/M;

%% numerical solution for J, f, wp, V, U, A and Q

if matlab_ver_less_than('9.1')
    options=optimset('TolFun',1e-07,'TolX',1e-07);
else
    options=optimoptions('fsolve','FunctionTolerance',1e-07,'OptimalityTolerance',1e-07);
end
options.MaxFunctionEvaluations = 100e+02;

[outvalue,fval,exitflag]=fsolve(@(par_vec)solve_calibration(par_vec,rhos,betta,nup,nu,zeta1,zeta2,zeta3,zeta4,gammac,f),[270.4;271.2;270.3;0.90],options);
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
    check=1; %set failure indicator
    return; %return without updating steady states
end

U=outvalue(1); % (22)
V=outvalue(2); % (23)
A=outvalue(3); % (24)
wp=outvalue(4); % (25)

% set V_ss
V_ss=V;

% set U_ss
U_ss=U;

% w (26)
w=(1-(1-rhos)*betta)*wp;

% set w_ss
w_ss=w;

% set UBbar (p8)
UBbar=0.4*(1-(1-rhos)*betta)*wp; % 0.4 from replacement ratio UBbar/w

% J (27)
J=nup-wp;

% set J_ss
J_ss=J;

% k (28)
k=alppha/(1-alppha)*nu/ren*til_Phi*til_Psi*l/iota;
% I (29)
I=k*(1-(1-delk)/(til_Phi*til_Psi));

% yd (30)
yd=(iota^alppha*l^(1-alppha)*(k/(til_Phi*til_Psi))^alppha)/(0.42+Xi);

% phif (p11) % phif/yd = 0.42 from estimated value in Christiano et.al.(2016)
phif=0.42*yd;

% set y_ss
y_ss=yd;

% set kappa (p12) (from etah/(etah+etas)=0.94 in Christiano et.al.(2016))
kappa=0.94*J;

% set etah
etah=mat*kappa/yd;

% set s (p13)
s=0.06*J*Q;

% set etas
etas=s*va/yd;


% set Gbar (0.2 from calibrated SS Gbar/yd ratio) (p7)
Gbar = 0.2*yd;

% c (31)
c=yd-I-Gbar-til_Phi^(-1)*til_Psi^(-1)*(gamma1*(iota-1)+gamma2/2*(iota-1)^2)*k-(s/Q+kappa)*mat;
% lambda (32)
lambda=(1-h*betta*til_Phi^(-1))*(1-h*til_Phi^(-1))^(-1)*c^(-1);

% b (0.40 : calibrated SS b/GDP ratio; GDP = c+I+Gbar) (33)
b=0.40*(c+I+Gbar);

%set tilb (p6)
tilb=b;

tau=Gbar+(1-l)*UBbar+b*PI^(-1)*til_Phi^(-1)-b/R; % (34)

%set tiltau (p5)
tiltau=tau;

% g1 (35)
g1=lambda*mc*yd/(1-betta*thetap*PI^epsilon);
% g2 (36)
g2=epsilon/(epsilon-1)*g1;

%% The conditions that need to be satisfied
residuals=NaN(4,1);

residuals(1)=V-wp-A; % (eq 19) 
residuals(2)=(1-(1-rhos)*betta)*A-rhos*betta*(f*V+(1-f)*U); % (eq 20)
residuals(3)=0.4*(1-(1-rhos)*betta)*wp+betta*f*V+(betta*(1-f)-1)*U; % (eq 21) plugged with (p8) at steadystate.m
residuals(4)=wp-(zeta1*nup+zeta2*(U-A)+zeta3*gammac-zeta4*(nu-0.4*(1-(1-rhos)*betta)*wp))/(zeta1+zeta2);% (eq 24) plugged with (p8) at steadystate.m
if max(abs(residuals)>1e-6)
    check=1; %set failure indicator
    return; %return without updating steady states
end


%% end own model equations

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  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 ';']);
end
