function [ys,check] = paper1_steadystate(ys,exo)
% function [ys,check] = NK_baseline_steadystate(ys,exo)
% computes the steady state for the NK_baseline.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
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;
options=optimset('display','iter');

eps_gs=1;eps_f=0;
eps_c=1;eps_s=1;
eps_H=1;eps_l=1; 
eps_Fs=1;eps_ls=1;
PI=PIbar;PI_s=PIsbar;
PI_H=PI;PI_Fs=PI_s;
R=PI/bet;de=PI/PI_s;
R_s=R/(de*(1-phif));
Rbar=R;Rsbar=R_s;
r=(1/bet)-(1-delta);
r_s=(R_s/PI_s)-(1-deltas);

%options=optimset('display','iter');

x=[5.9 2.6 3.1 1.5 48 7 29];

x=fsolve(@(x) [gammmah*((phi_H/(phi_H-1))*((1-bet*thetah*PI^(phi_H-1))/ ...
(1-bet*thetah*PI^phi_H))*((1-thetah)/(1-thetah*PI^ ...
((delta_H-1)*(1-phi_H))))^(1/(1-phi_H))*alp^(-alp)* ...
(1-alp)^(-(1-alp))*r^alp*x(1)^(1-alp))^(1-theta)+ ...
(1-gammmah)*((x(4)^phi_s*x(1)/(x(3)^phi*x(2)))* ...
((phi_Fs/(phi_Fs-1))*((1-bet_s*thetaf*PI_s^(phi_Fs-1))/ ...
(1-bet_s*thetaf*PI_s^phi_Fs))*((1-thetaf)/(1-thetaf*PI_s^ ...
((delta_F-1)*(1-phi_Fs))))^(1/(1-phi_Fs))*alps^(-alps)* ...
(1-alps)^(-(1-alps))*r_s^alps*x(2)^(1-alps)))^(1-theta)-1;

gammmahs*((phi_Fs/(phi_Fs-1))*((1-bet_s*thetaf*PI_s^(phi_Fs-1))/ ...
(1-bet_s*thetaf*PI_s^phi_Fs))*((1-thetaf)/(1-thetaf*PI_s^ ...
((delta_F-1)*(1-phi_Fs))))^(1/(1-phi_Fs))*alps^(-alps)* ...
(1-alps)^(-(1-alps))*r_s^alps*x(2)^(1-alps))^(1-thetas)+ ...
(1-gammmahs)*(((phi_H/(phi_H-1))*((1-bet*thetah*PI^(phi_H-1))/ ...
(1-bet*thetah*PI^phi_H))*((1-thetah)/(1-thetah*PI^ ...
((delta_H-1)*(1-phi_H))))^(1/(1-phi_H))*alp^(-alp)* ...
(1-alp)^(-(1-alp))*r^alp*x(1)^(1-alp))/(x(4)^phi_s*x(1)/ ...
(x(3)^phi*x(2))))^(1-thetas)-1;

gammmah*((phi_H/(phi_H-1))*((1-bet*thetah*PI^(phi_H-1))/ ...
(1-bet*thetah*PI^phi_H))*((1-thetah)/(1-thetah*PI^ ...
((delta_H-1)*(1-phi_H))))^(1/(1-phi_H))*alp^(-alp)* ...
(1-alp)^(-(1-alp))*r^alp*x(1)^(1-alp))^(-theta)*x(5)+ ...
(1-gammmahs)*(((phi_H/(phi_H-1))*((1-bet*thetah*PI^(phi_H-1))/ ...
(1-bet*thetah*PI^phi_H))*((1-thetah)/(1-thetah*PI^ ...
((delta_H-1)*(1-phi_H))))^(1/(1-phi_H))*alp^(-alp)* ...
(1-alp)^(-(1-alp))*r^alp*x(1)^(1-alp))/(x(4)^phi_s* ...
x(1)/(x(3)^phi*x(2))))^(-thetas)*x(6)-((alp*x(1)/((1-alp)*r)))^alp*x(3);

(1-gammmah)*((x(4)^phi_s*x(1)/(x(3)^phi*x(2)))*((phi_Fs/(phi_Fs-1))* ...
((1-bet_s*thetaf*PI_s^(phi_Fs-1))/ ...
(1-bet_s*thetaf*PI_s^phi_Fs))*((1-thetaf)/(1-thetaf*PI_s^ ...
((delta_F-1)*(1-phi_Fs))))^(1/(1-phi_Fs))*alps^(-alps)* ...
(1-alps)^(-(1-alps))*r_s^alps*x(2)^(1-alps)))^(-theta)*x(5)+gammmahs* ...
((phi_Fs/(phi_Fs-1))*((1-bet_s*thetaf*PI_s^(phi_Fs-1))/ ...
(1-bet_s*thetaf*PI_s^phi_Fs))*((1-thetaf)/(1-thetaf*PI_s^ ...
((delta_F-1)*(1-phi_Fs))))^(1/(1-phi_Fs))*alps^(-alps)* ...
(1-alps)^(-(1-alps))*r_s^alps*x(2)^(1-alps))^(-thetas)*x(6)- ...
((alps*x(2)/((1-alps)*r_s)))^alps*x(4);

(1-gammmahs)*((phi_H/(phi_H-1))*((1-bet*thetah*PI^(phi_H-1))/ ...
(1-bet*thetah*PI^phi_H))*((1-thetah)/(1-thetah*PI^ ...
((delta_H-1)*(1-phi_H))))^(1/(1-phi_H))*alp^(-alp)* ...
(1-alp)^(-(1-alp))*r^alp*x(1)^(1-alp))*(((phi_H/(phi_H-1))* ...
((1-bet*thetah*PI^(phi_H-1))/(1-bet*thetah*PI^phi_H))* ...
((1-thetah)/(1-thetah*PI^((delta_H-1)*(1-phi_H))))^ ...
(1/(1-phi_H))*alp^(-alp)*(1-alp)^(-(1-alp))*r^alp*x(1)^ ...
(1-alp))/(x(4)^phi_s*x(1)/(x(3)^phi*x(2))))^(-thetas)*x(6)- ...
((x(4)^phi_s*x(1)/(x(3)^phi*x(2)))*((phi_Fs/(phi_Fs-1))* ...
((1-bet_s*thetaf*PI_s^(phi_Fs-1))/ ...
(1-bet_s*thetaf*PI_s^phi_Fs))*((1-thetaf)/(1-thetaf*PI_s^ ...
((delta_F-1)*(1-phi_Fs))))^(1/(1-phi_Fs))*alps^(-alps)* ...
(1-alps)^(-(1-alps))*r_s^alps*x(2)^(1-alps)))^(1-theta)*(1-gammmah)*x(5);

((x(3)^(-phi)*(1-h)^(-sigmau)*x(1))^(1/sigmau))+ ...
(delta*(alp*x(1)/((1-alp)*r))*x(3))+x(7)-x(5);

((x(4)^(-phi_s)*(1-h_s)^(-sigma_s)*x(2))^(1/sigma_s))+ ...
(deltas*(alps*x(2)/((1-alps)*r_s))*x(4))+eps_gs-x(6)],x);

w=x(1);
w_s=x(2);
l=x(3);
l_s=x(4);
y=x(5);
y_s=x(6);
g=x(7);

ph=(phi_H/(phi_H-1))*((1-bet*thetah*PI^(phi_H-1))/ ...
(1-bet*thetah*PI^phi_H))*((1-thetah)/(1-thetah*PI^ ...
((delta_H-1)*(1-phi_H))))^(1/(1-phi_H))*alp^(-alp)* ...
(1-alp)^(-(1-alp))*r^alp*w^(1-alp);

pfs=(phi_Fs/(phi_Fs-1))*((1-bet_s*thetaf*PI_s^(phi_Fs-1))/ ...
(1-bet_s*thetaf*PI_s^phi_Fs))*((1-thetaf)/(1-thetaf*PI_s^ ...
((delta_F-1)*(1-phi_Fs))))^(1/(1-phi_Fs))*alps^(-alps)* ...
(1-alps)^(-(1-alps))*r_s^alps*w_s^(1-alps);

rer=l_s^phi_s*w/(l^phi*w_s);
mc_h=alp^(-alp)*(1-alp)^(-(1-alp))*r^alp*w^(1-alp);
mc_f=alps^(-alps)*(1-alps)^(-(1-alps))*r_s^alps*w_s^(1-alps);
k=(alp*w/((1-alp)*r))*l;
y_H=k^alp*l^(1-alp);
i=delta*k;
c=(l^(-phi)*(1-h)^(-sigmau)*w)^(1/sigmau);
c_s=(l_s^(-phi_s)*(1-h_s)^(-sigma_s)*w_s)^(1/sigma_s);
i_s=deltas*(alps*w_s/((1-alps)*r_s))*l_s;
k_s=i_s/deltas;
mm=((1-bet/PI)*l^phi/x(1))^(-1/sigmam);
m_s=((1-bet_s/PI_s)*l_s^phi_s/w_s)^(-1/sigmams);
ca=0;ca_s=0;      
f=0;f_s=0;
fbar=f;
tb=0;tb_s=0;
gdp=y;gdp_s=y_s;
gdpsbar=gdp_s;
gdpbar=gdp;
y_Fs=k_s^alps*l_s^(1-alps);
q=1;q_s=1;V=PI;
phn=(c-h*c)^(-sigmau)*ph^phi_H*y_H/(1-bet*thetah*PI^(phi_H-1));
phd=(phi_H/(phi_H-1))*((c-h*c)^(-sigmau)*ph^phi_H*y_H*mc_h)/(1-bet*thetah*PI^phi_H);
pfns=(c_s-h_s*c_s)^(-sigma_s)*pfs^phi_Fs*y_Fs/(1-bet_s*thetaf*PI_s^(phi_Fs-1));
pfds=(phi_Fs/(phi_Fs-1))*((c_s-h_s*c_s)^(-sigma_s)* ...
pfs^phi_Fs*y_Fs*mc_f)/(1-bet_s*thetaf*PI_s^phi_Fs);
gbar=g;
d=(R/PI-1)^(-1)*(omega*gdp+(1-1/PI)*mm-g);
dbar=d;Vbar=V;

for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end