%compute the steady state of NK Model 

function [ys,params,check]=BGG_financial_3_p45_phi2_steadystate(~,~,M_,~)
global M_
%% DO NOT CHANGE THIS PART.
%%
%% Here we load the values of the deep parameters in a loop.
%%
% NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
% for i = 1:NumberOfParameters                                  % Loop...
%   paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
%   disp([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
% end                                                           % End of the loop.  
check = 0;
%%

gy = M_.params(1);
alp = M_.params(2);
delta = M_.params(3);
zzetap = M_.params(4);
zzetaw = M_.params(5);
sigma_c = M_.params(6);
chi = M_.params(7);
phiX = M_.params(8);
xip = M_.params(9);
xiw = M_.params(10);
gammap = M_.params(11);
gammaw = M_.params(12);
PIE_bar = M_.params(13);
g_bar = M_.params(14);
Rn_bar = M_.params(15);
rho_r = M_.params(16);
theta_pie = M_.params(17);
theta_y = M_.params(18);
theta_dy = M_.params(19);
rhoA = M_.params(20);
rhoG = M_.params(21);
rhoMCS = M_.params(22);
rhoMRSS = M_.params(23);
rhoMPS = M_.params(24);
rhoRPS = M_.params(25);
rhoIS = M_.params(26);
ptarget = M_.params(27);
rhotarget = M_.params(28);
phitarget = M_.params(29);
CEtarget = M_.params(30);
sigma_s = M_.params(31);
wr = M_.params(32);
apsi = M_.params(33);

%% THIS BLOCK IS MODEL SPECIFIC.
%%
%% Here the user has to define the steady state.
%%

H0        = 0.8994;
KYW0      = 7.7627;
sigmaE0   = 0.99;
sigmapsi0 = 1.2729;
psi0      = 0.0406;
muFF0     = 0.1276;
xiE0      = 0.0019;
% x0        = [log(H0); log(KYW0); log(sigmaE0/(1-sigmaE0)) ;  log(sigmapsi0); log(psi0);log(muFF0/(1-muFF0)); log(xiE0/(1-xiE0))];
x0        = [log(H0); log(KYW0); log(sigmapsi0); log(psi0);log(muFF0/(1-muFF0)); log(xiE0/(1-xiE0))];
options = optimset('TolFun',1e-9,'TolX',1e-9,'MaxIter', 5000, 'MaxFunEvals', 5000, 'Display','off');
x       = fsolve(@ss_fun_BGG3,x0,options);

H        = exp(x(1));
KYW      = exp(x(2));
% sigmaE   = exp(x(3))/(1+exp(x(3)));
% sigmapsi = exp(x(4));
% psiFF    = exp(x(5));
% muFF     = exp(x(6))/(1+exp(x(6))); 
% xiE      = exp(x(7))/(1+exp(x(7)));
sigmapsi = exp(x(3));
psiFF    = exp(x(4));
muFF     = exp(x(5))/(1+exp(x(5))); 
xiE      = exp(x(6))/(1+exp(x(6)));
sigmaE   = 0.99;

g       = g_bar/100;
betta   = (PIE_bar/100+1)/(Rn_bar/100+1)/(1+g)^(-sigma_c); % (273)

A       = 1;
MCS     = 1;
MRSS    = 1;
IS      = 1;
MPS     = 1;
RPS     = 1;

PIE     = PIE_bar/100+1;

S       = 0; % (462)
Sdash   = 0; % (463)
PIEw    = (1+g)*PIE; % (464)
Q       = 1; % (465)
LAMBDA  = betta*(1+g)^(-sigma_c); % (466)
PO      = ((1-xip*PIE^((gammap-1)*(1-zzetap)))/(1-xip))^(1/(1-zzetap)); % (467)
Deltap  = (1-xip)/(1-xip*PIE^(zzetap*(1-gammap)))*PO^(-zzetap); % (468)
MC      = (1-1/zzetap)*1/MCS*(1-xip*(1+g)*LAMBDA*PIE^(zzetap*(1-gammap)))/(1-xip*(1+g)*LAMBDA*PIE^((zzetap-1)*(1-gammap)))*PO; % (469)
PWP     = MC; % (470)
WO      = ((1-xiw*PIE^(gammaw*(1-zzetaw))*PIEw^(zzetaw-1))/(1-xiw))^(1/(1-zzetaw)); % (471)
Deltaw  = (1-xiw)/(1-xiw*PIEw^(zzetaw)/PIE^(gammaw*zzetaw))*WO^(-zzetaw); % (472)
Hd      = H/Deltaw; % (473)
rK      = (1-alp)*(1+g)/KYW*PWP; % (474)
YW      = A*Hd*(KYW/(1+g))^((1-alp)/alp); % (475)
K       = YW*KYW; % (476)
Y       = YW/Deltap; % (477)
I       = K/IS*(g+delta)/(1+g); % (478)
G       = gy*Y; % (479)
RK      = rK+(1-delta); % (480)
JJw     = Hd/(1-xiw*LAMBDA*PIEw^zzetaw*PIE^(gammaw*(1-zzetaw)-1)); % (481)
W       = alp*PWP*YW/Hd; % (482)
Jw      = WO*W*JJw; % (483)
Wh      = (1-xiw*(1+g)*LAMBDA*PIEw^zzetaw/PIE^(gammaw*zzetaw))*(1-1/zzetaw)*Jw/(Hd*MRSS); % (484)
mupsi          = -0.5*(sigmapsi)^2;
p              = logncdf(psiFF,mupsi,sigmapsi);
Gofpsi         = normcdf(log(psiFF)/sigmapsi-sigmapsi/2);
Gammaofpsi     = Gofpsi+psiFF*(1-p);
Gammadashofpsi = 1-p;
Gdashofpsi     = 1/(psiFF*sigmapsi*sqrt(2*(4*atan(1))))*exp(-(log(psiFF)+mupsi)^2/(2*(sigmapsi)^2));
CE        = (1-sigmaE)*(1-xiE)*(1-Gammaofpsi)*RK*Q*K/(1+g); % (485)
rhoofpsi  = Gammadashofpsi/((Gammaofpsi-muFF*Gofpsi)*Gammadashofpsi+(1-Gammaofpsi)*(Gammadashofpsi-muFF*Gdashofpsi)); % (486)
C         = Y-G-I-CE-muFF*Gofpsi*RK*Q*K/(1+g); % (487)
NE        = (sigmaE+xiE)*(1-Gammaofpsi)*RK*Q*K/(1+g); % (488)
L         = Q*K-NE; % (489)
Rl        = psiFF*RK*Q*K*PIE/L; % (490)
R         = 1/(LAMBDA); % (498)

U       = (C*(1-chi/(1+g)))^(1-sigma_c)/(1-sigma_c)*exp((sigma_c-1)*H^(1+apsi)/(1+apsi)); % (494)
V       = U/(1-betta*(1+g)^(1-sigma_c)); % (495)
Vmodified = V; % (495a)
X       = 1+g; % (496)
lambda  = (1-sigma_c)*U/(C-chi*C/(1+g))-betta*chi*(1-sigma_c)*U/(1+g)^(1-sigma_c)/((1+g)*C-chi*C); % (497)
Rn      = R*PIE; % (499)
spread  = RK-R; % (500)
Jp      = Y*MC*MCS/((1-1/zzetap)*(1-xip*(1+g)*LAMBDA*PIE^(zzetap*(1-gammap)))); % (501)
JJp     = Jp/PO; % (502)
CEquiv  = (1.01^(1-sigma_c)-1)*U/(1-betta*(1+g)^(1-sigma_c)); % (503)
UH      = -(1-sigma_c)*U*H^apsi;

dyobs   = log((1+g));
dcobs   = log((1+g));
dinvobs = log((1+g));
dwobs   = log((1+g));
labobs  = 0;
robs    = (Rn-1);
pinfobs = (log(PIE));
spreadobs = spread;


%% END OF THE MODEL SPECIFIC BLOCK.

%% DO NOT CHANGE THIS PART.
%%
%% Here we define the steady state vZNues of the endogenous variables of
%% the model.
%%
NumberOfEndogenousVariables = M_.endo_nbr;                    % Number of endogenous variables.
ys = zeros(NumberOfEndogenousVariables,1);                    % Initialization of ys (steady state).
% for i = 1:NumberOfEndogenousVariables                         % Loop...
%   varname = deblank(M_.endo_names(i,:));                      %    Get the name of endogenous variable i.                     
%   disp(['ys(' int2str(i) ') = ' varname ';']);                %    Get the steady state vZNue of this variable.
% end                                                           % End of the loop.
ys(1) = V;
ys(2) = Vmodified;
ys(3) = U;
ys(4) = C;
ys(5) = H;
ys(6) = K;
ys(7) = S;
ys(8) = X;
ys(9) = I;
ys(10) = IS;
ys(11) = Sdash;
ys(12) = lambda;
ys(13) = Wh;
ys(14) = rK;
ys(15) = RPS;
ys(16) = R;
ys(17) = Q;
ys(18) = Rn;
ys(19) = PIE;
ys(20) = PIEw;
ys(21) = W;
ys(22) = Jw;
ys(23) = Hd;
ys(24) = MRSS;
ys(25) = JJw;
ys(26) = WO;
ys(27) = Deltaw;
ys(28) = YW;
ys(29) = A;
ys(30) = PWP;
ys(31) = MC;
ys(32) = Jp;
ys(33) = Y;
ys(34) = MCS;
ys(35) = JJp;
ys(36) = PO;
ys(37) = Deltap;
ys(38) = MPS;
ys(39) = G;
ys(40) = RK;
ys(41) = spread;
ys(42) = CEquiv;
ys(43) = rhoofpsi;
ys(44) = psiFF;
ys(45) = NE;
ys(46) = CE;
ys(47) = L;
ys(48) = Rl;
ys(49) = p;
ys(50) = Gofpsi;
ys(51) = Gammaofpsi;
ys(52) = Gammadashofpsi;
ys(53) = Gdashofpsi;
ys(54) = sigmaE;
ys(55) = xiE;
ys(56) = muFF;
ys(57) = sigmapsi;
ys(58) = dyobs;
ys(59) = dcobs;
ys(60) = pinfobs;
ys(61) = robs;
ys(62) = labobs;
ys(63) = spreadobs;
ys(64) = dinvobs;
ys(65) = dwobs;
ys(66) = UH;


% Set params
params = M_.params;
