function [ys,check] = Thesis_nlnbl_steadystate(ys,exo)
% function [ys,check] = Thesis_nlnbl_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

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 

% initialize indicator
check = 0;

%% Enter model equations here

options=optimset(); % set options for numerical solver

% Foreign debt interest rate
r_f    = rstar+pbar;

% CPI inflation equals target
pi_cpi = pi_bar;

% Steady state of r_ib equals target
r_ib   = r_ib_bar;

% Bank interest rates
R_b    = r_ib;
r_b    = R_b+muy_b;

% Set parameter delt_b
delt_b = (muy_b+nuy*r_ib_bar)/nuy;

% exchange rate depreciation
ge = pi_bar;

% exogenous processes
% Technology Shock
a      = 1;

% Foreign output Shock
ystar  = 1;

% Monetary shock
u      = 1;

% Domestic goods & imported goods inflation & relative price
pi_h   = pi_bar;
pi_f   = pi_bar;
pih   = 0;
pif   = 0;

% Domestic & Imported good optimal price
pi_hr   = (((1+pi_h)^(1-epsl)-thetta_h)/(1-thetta_h))^(1/(1-epsl))-1;
pi_fr   = (((1+pi_f)^(1-epsl)-thetta_f)/(1-thetta_f))^(1/(1-epsl))-1;

% Domestic & import price dispersion
%vh = ((1-thetta_h)*((1+pi_h)/(1+pi_hr))^epsl)/(1-thetta_h*(1+pi_h)^(epsl));
vf = ((1-thetta_f)*((1+pi_f)/(1+pi_fr))^epsl)/(1-thetta_f*(1+pi_f)^(epsl));

% Real marginal cost of producers
mc = (1-betta*thetta_h*(1+pi_h)^(epsl+1)/(1+pi_cpi))/(1-betta*thetta_h*(1+pi_h)^(epsl)/(1+pi_cpi))*(1+pi_hr)/(1+pi_h)*(epsl-1)/epsl;

% lopg Real marginal cost of importers
lopg = (1-betta*thetta_f*(1+pi_f)^(epsl+1)/(1+pi_cpi))/(1-betta*thetta_f*(1+pi_f)^(epsl)/(1+pi_cpi))*(1+pi_fr)/(1+pi_f)*(epsl-1)/epsl;

% Foreign demand for domestic goods
c_hstar = gam*lopg^etta*ystar;

% FOC investment Tobin q
q      = mc;

% Real wage
w      = (1-alppha)*(1+pih)*mc*a^(1/(1-alppha))*(alppha/(r_b+delt_k))^(alppha/(1-alppha));

% {solve numerically for n}
[n,fval,exitflag]=fzero(@(n)n^(varrphi)/w-((1-gam)/(a*n*(w*alppha/(r_b+delt_k)/(1-alppha)/(1+pih)/mc)^alppha-c_hstar))^siggma,0.8066,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


% capital of producers
k      = n*w*alppha/(r_b+delt_k)/(1-alppha)/(1+pih)/mc;

% Law of motion of capital of domestic producers
invest  = delt_k*k;

% Production function
y      = a*k^alppha*n^(1-alppha);

% Local demand for domestic goods
c_h     = y-c_hstar;

% Total consumption
c = c_h/(1-gam);

% Local demand for imported goods
c_f     = gam*c;

% Credit market clearing 
b      = k;

% Bank balance sheet
m      = b*(1-nuy);

% Leverage
k_b    = b*nuy;

% Bank profit
j_b = delt_b*k_b;
% j_b can also be determined by j_b = r_b*b-r_ib*m; however, delt_b is not set yet, so these 2 equations together are used to set delt_b as above

% Debt
d = ((r_ib-pi_cpi)/(1+pi_cpi)*m-(1+pih)*mc*(r_b*b+invest)+(1+pih)*c_hstar-lopg*(1+pif)*c_f*vf)/((r_f-pi_cpi)/(1+pi_cpi)*(1+pif)*lopg);

%set the parameter d_bar
d_bar   = d;
%if d_bar<0 % parameter violates restriction; Preventing this cannot be implemented via prior restriction as it is a composite of different parameters and the valid prior region has unknown form
%    check=1; %set failure indicator
%    return; %return without updating steady states
%end


x1=y*mc/(1-betta*thetta_h*(1+pi_h)^(epsl+1)/(1+pi_cpi));

x2=y/(1-betta*thetta_h*(1+pi_h)^(epsl)/(1+pi_cpi));

g1=c_f*lopg/(1-betta*thetta_f*(1+pi_f)^(epsl+1)/(1+pi_cpi));

g2=c_f/(1-betta*thetta_f*(1+pi_f)^(epsl)/(1+pi_cpi));

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



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

