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;
% Real profits of intermediate firms (6)
F=0;


% LambdaYd (p4)
LambdaYd=alppha*LambdaPsi/(1-alppha)+LambdaZ;
% SS growth til_Phi (7)
til_Phi=exp(LambdaYd);
% SS growth til_Psi (8)
til_Psi=exp(LambdaPsi);
% SS growth til_Z (9)
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 (10)
ren=gamma1+gamma2*(iota-1);
% R (11)
R=PI*exp(LambdaYd)/betta;

%set Rbar (p2)
Rbar=R;

% PIstar (12)
PIstar=((1-thetap*(PI)^(epsilon-1))/(1-thetap))^(1/(1-epsilon));
% Xi (13)
Xi=(1-thetap)*(1-thetap*(PI)^epsilon)^(-1)*PIstar^(-epsilon);
% mc (14)
mc=(epsilon-1)/epsilon*(1-betta*thetap*PI^epsilon)/(1-betta*thetap*PI^(epsilon-1))*PIstar;
% nu (15)
nu=(1-alppha)*(mc*(alppha/ren)^alppha)^(1/(1-alppha));
% nup (16)
nup=(1-(1-rhos)*betta)^(-1)*nu;

% for zeta1-4, and zeta5-7(beta1-3 in Christiano et.al.(2016))
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;
zeta5  = zeta2/zeta1;
zeta6  = zeta3/zeta1;
zeta7  = zeta4/zeta1;


%% 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,s,kappa,nup,nu,zeta5,zeta6,zeta7,gammac,sigmm,sigmaa),[0.06;0.63;0.90;270.4;271.2;270.3;0.7],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

J=outvalue(1); % (17)
f=outvalue(2); % (18)
wp=outvalue(3); % (19)
U=outvalue(4); % (20)
V=outvalue(5); % (21)
A=outvalue(6); % (22)
Q=outvalue(7); % (23)


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

% Gamma (25) : combined (eq 14), (eq 22) and (eq 24)
Gamma=f/Q;
% l (26) : combined (eq 12), (eq 13), (eq 23), and (eq 24)
l=sigmm*Gamma^(sigmaa)*(rhos+(1-rhos)*sigmm*Gamma^(sigmaa))^(-1);
% ump (27)
ump=1-(1-rhos)*l;
% mat (28) : from (eq 12)
mat=rhos*l;
% va (29)
va=mat/Q;
% k (30)
k=alppha/(1-alppha)*nu/ren*til_Phi*til_Psi*l*iota^(-1);
% I (31)
I=k*(1-(1-delk)*til_Phi^(-1)*til_Psi^(-1));

% yd (From F=0 at SS) (32)
yd=1/(1-alppha)*nu*l;
% set Gbar (0.2 from calibrated SS Gbar/yd ratio) (p7)
Gbar = 0.2*yd;
% phif (phif to ensure F is 0 at SS) (p9)
phif=iota^alppha*l^(1-alppha)*(k*til_Phi^(-1)*til_Psi^(-1))^alppha-Xi*yd;

% c (33)
c=yd-I-Gbar-til_Phi^(-1)*til_Psi^(-1)*(gamma1*(iota-1)+gamma2/2*(iota-1)^2)*k-(s/Q+kappa)*mat;
% lambda (34)
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) (35)
b=0.40*(c+I+Gbar);

%set tilb (p6)
tilb=b;

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

%set tiltau (p5)
tiltau=tau;

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

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

residuals(1)=J-nup+wp; % (eq 15)
residuals(2)=J-(s/Q+kappa); % (eq 18)
residuals(3)=V-wp-A; % (eq 19) 
residuals(4)=(1-(1-rhos)*betta)*A-rhos*betta*(f*V+(1-f)*U); % (eq 20)
residuals(5)=0.4*(1-(1-rhos)*betta)*wp+betta*f*V+(betta*(1-f)-1)*U; % (eq 21) plugged with (p8) at steadystate.m
residuals(6)=J-zeta5*(V-U)-zeta7*nu+(zeta6*gammac+zeta7*0.4*(1-(1-rhos)*betta)*wp); % (eq 25) plugged with (p8) at steadystate.m
residuals(7)=Q-sigmm^(1/(1-sigmaa))*f^(-sigmaa/(1-sigmaa)); % (eq 14), (eq 22), and (eq 23)

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) ') =  0 ;']);
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
