%compute the steady state of NK Model 

function [ys,params,check]=GK_equity_MPR_steadystate(ys,exo,M_,options_);

%%
%% DO NOT CHANGE THIS PART.
%%
%% Here we load the values of the deep parameters in a loop.
%%
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;
%%


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

rhoH_0=0.2;
rhoHF_0=0.2;
theta_0=2;
x0=[log(rhoH_0/(1-rhoH_0)); log(theta_0); log(5.0); log(0.47); log(0.0025); log(rhoHF_0/(1-rhoHF_0));  ];
options = optimset('TolFun',1e-9,'TolX',1e-9,'MaxIter', 5000, 'MaxFunEvals', 5000, 'Display','on');
[x,~] = fsolve(@ss_fun_GK_equity_MPR,x0,options);
rhoH=exp(x(1)/(1+exp(x(1))));
theta=1+exp(x(2));
KYW=exp(x(3));
thetaFF=exp(x(4));
xiB=exp(x(5));
rhoHF=exp(x(6)/(1+exp(x(6))));


g=g_bar;
PIE=PIE_bar;

PIEw=PIE_bar*(1+g);

if index==1
%
%with complete indexing
%
PIEtilde=1;
PIEwtilde=1;
PIEwptilde=1;
else
PIEtilde=PIE^(1-gammap);
PIEwtilde=PIEw/PIE^gammaw;
PIEwptilde=PIE^(1-gammaw);
end
A=A_bar;
H=H_bar;
kappa=rhoH/H^theta;
bettag=betta*(1+g)^(-sigma_c);
R=1.0/bettag;


PO=((1-xip*PIEtilde^(zzeta-1))/(1-xip))^(1/(1-zzeta));
MC=(1-1/zzeta)*PO*(1-bettag*(1+g)*xip*PIEtilde^zzeta)/(1-bettag*(1+g)*xip*PIEtilde^(zzeta-1));

if wage>0
WO=((1-xiw*PIEwtilde^(muw-1))/(1-xiw))^(1/(1-muw));
Deltaw=(1-xiw)/(1-xiw*PIEwtilde^muw)*WO^(-muw);
Invmarkup=(1-1/muw)*WO*(1-bettag*xiw*(1+g)*PIEwtilde^muw)/(1-bettag*xiw*PIEwtilde^muw*PIEwptilde^(-1));

else
%flexi wages
%
WO=1;
Deltaw=1;
Invmarkup=1;   
end  




PWP=MC;
Deltap=((1-xip)*(((1-xip*PIEtilde^(zzeta-1))/(1-xip))^(1/(1-zzeta)))^(-zzeta))/(1-xip*PIEtilde^zzeta);

IY=(delta+g)*KYW/(1+g)*Deltap;

CY=1-IY-gy;
Hd=H/Deltaw;
YW=A*Hd*(KYW/(1+g))^((1-alp)/alp);
Y=YW/Deltap;
K=KYW*YW;
I=IY*Y;
G=gy*Y;
W=alp*PWP*YW/Hd;
Wh=W*Invmarkup;
C=CY*Y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JR Preferences Block with Corrected Frisch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Frisch (corrected NEW)
%
F=(1-rhoH*gam);
F1=(1-rhoH);
GG=rhoH;
AA=-gam+sigma_c*F/F1;
B=(sigma_c*theta*GG+(theta-1)*F1)/(sigma_c*F-gam*F1);
CCC=theta*GG*(sigma_c*F-gam*F1)/(sigma_c*F^2+gam*(gam-1)*GG*F1*(1+1/(1-bettag*(1-gam))));
Frisch=AA*(B-CCC);

%
%Check Leisure is a normal good if checknormalgood>0 where
%
%checknormalgood=kappa*H_bar^theta+theta-1; removed
%
Xi=C*(1-chi/(1+g))/(1+g)^((1-gam)/gam);
U=((C*(1-chi/(1+g))-kappa*(H^theta)*Xi)^(1-sigma_c)-1)/(1-sigma_c);
CEquiv=(((1.01*(C*(1-chi/(1+g)))-kappa*(H^theta)*Xi)^(1-sigma_c)-1)/(1-sigma_c)-U)/(1-bettag*(1+g));
Omega=U/(1-bettag*(1+g));
UC=(C*(1-chi/(1+g))-kappa*H^theta*Xi)^(-sigma_c);
UH=-theta*kappa*H^(theta-1)*Xi*UC;
UXi=-kappa*H^theta*UC;
%
%mu=-UXi/(1-betta*(1-gam));
mu=-UXi/(1-bettag*(1+g)*(1-gam));
%
if  habit==1
%
%external habit
%
lam=UC-gam*mu*Xi/(C*(1-chi/(1+g)));
%
else
%
%internal habit
%
%lam=UC*(1-betta*chi)-gam*mu*Xi*(1-betta*chi)/(C*(1-chi/(1+g)));
lam=UC*(1-bettag*chi)-gam*mu*Xi*(1-bettag*chi)/(C*(1-chi/(1+g)));
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


LAMBDA=bettag;
S=0;
Sdash=0;
Q=1;
X=1+g;

JJp=Y/(1-bettag*(1+g)*xip*PIEtilde^(zzeta-1));
JJtilde=(PIEtilde^(zzeta-1))*JJp;
Jp=(1/(1-1/zzeta))*Y*MC/(1-bettag*(1+g)*xip*PIEtilde^(zzeta));
Jtilde=PIEtilde^zzeta*Jp;

if wage>0
JJw=Hd/(1-bettag*xiw*PIEwtilde^muw/PIEwptilde);
Jw=(1/(1-(1/muw)))*Wh*Hd/(1-bettag*xiw*(1+g)*PIEwtilde^muw);
else
JJw=Hd;
Jw=Hd*Wh;   
end


MCS=1;
MRSS=1;
MPS=1;
RPS=1;

IS=1;
KQ=1;
Rn=1/bettag*PIE;
ER=1/bettag;

%
%Banking Sector
%
Sasset=K;
Q=1;

RK=((1-alp)*PWP*YW/(K/(1+g))+(1-delta)*Q)/Q;
xequity=xequitytarg;
kappaFF=-epsilon/xequitytarg;
qE=xequity*Q*Sasset;
Re=R;
mue=0;
N=((RK*(sigmaB+xiB)-sigmaB*R)*Q*Sasset-sigmaB*(Re-R)*qE)/(1+g-sigmaB*R);
D=Q*Sasset-N-qE;
phiB=Q*Sasset/N;
lever=Q*Sasset/(N+qE);
Thetax=thetaFF*(1+epsilon*xequity+kappaFF*xequity^2/2);
omega=1-sigmaB+Thetax*sigmaB*phiB;
nud=omega;
mus=Thetax-nud/phiB;
q=(1-alp)*PWP*YW/(K/(1+g))/(Re-(1-delta));
lambda=(mus+mue*xequity)/(Thetax-mus-mue*xequity);
DThetax=0;
spread=RK-R;
E=qE/q;


spread=RK-R;
epsAtrendstate=0;

%
%Flexi price and wage bloc from RBC
%
HF=H_bar;
kappaF=rhoHF/HF^theta;
RF=R;
RKF=RF;
KYF=(1+g)*(1-alp)/(R-1+delta);
IYF=(1-alp)*(delta+g)/(R-1+delta);
YF=A*HF*(KYF/(1+g))^((1-alp)/alp);
CYF=1-IYF-gy*Y/YF;


KF=KYF*YF;
IF=IYF*YF;
WF=alp*YF/HF;
CF=CYF*YF;
%
%JR Preferences
%
XiF=CF*(1-chi/(1+g))/(1+g)^((1-gam)/gam);
UCF=(CF*(1-chi/(1+g))-kappaF*HF^theta*XiF)^(-sigma_c);
UHF=-theta*kappaF*HF^(theta-1)*XiF*UCF;
UXiF=-kappaF*HF^theta*UCF;



muF=-UXiF/(1-bettag*(1+g)*(1-gam));
%
%habit
%
if  habit==1
%
%external habit
%
lamF=UCF-gam*muF*XiF/(CF*(1-chi/(1+g)));
%
else
%
%internal habit
%
lamF=UCF*(1-bettag*chi)-gam*muF*XiF*(1-bettag*chi)/(CF*(1-chi/(1+g)));
%
end

OUTGAP=Y/YF;


QF=Q;
LAMBDAF=LAMBDA;
SF=S;
SdashF=Sdash;
XF=X;








%% END OF THE MODEL SPECIFIC BLOCK.

%% DO NOT CHANGE THIS PART.
%%
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    


