function [ys,check] = TBDSGE_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;

%% Enter model equations here
options=optimset(); % set options for numerical solver

% Note: steady-state computations quite different from FVRR
PI=log(PIbar);
%d=log(1;
%phi=log(1);
d=1;
phi=1;
%d=0;
%phi=0;
m=0;



%NEW: r
R=log(exp(PI)/bet);
%Rbar=exp(PI)/bet;
%R=log(Rbar);
%r=log(exp(R)/exp(PI)-1+del);
r=log(1/bet-1+del);

Rbar=exp(PI)/bet;
R=log(Rbar);


%set Rbar
%log(Rbar)=log(PIbar/bet);
%log(exp(PI)/bet)=log(Rbar)

PIstar=log(((1-tetp*exp(PI)^(-(1-ela)*(1-chi)))/(1-tetp))^(1/(1-ela)));
PIstarw=log(((1-tetw*exp(PI)^(-(1-chiw)*(1-eta)))/(1-tetw))^(1/(1-eta)));

mc=log((ela-1)/ela*(1-bet*tetp*exp(PI)^((1-chi)*ela))/(1-bet*tetp*exp(PI)^(-(1-ela)*(1-chi)))*exp(PIstar));
w=log((1-alf)*(exp(mc)*(alf/exp(r))^alf)^(1/(1-alf)));
wstar=log(exp(w)*exp(PIstarw));
vp=log((1-tetp)/(1-tetp*exp(PI)^((1-chi)*ela))*exp(PIstar)^(-ela));
vw=log((1-tetw)/(1-tetw*exp(PI)^((1-chiw)*eta))*exp(PIstarw)^(-eta));
tempvaromega=alf/(1-alf)*exp(w)/exp(r);


%no log here because I think it has something to do with what f means
[ld,fval,exitflag]=fsolve(@(ld)(1-bet*tetw*exp(PI)^(-(1-chiw)*(1-eta)))/(1-bet*tetw*exp(PI)^(eta*(1-chiw)*(1+gam)))...
    -(eta-1)/eta*exp(wstar)/(psy*exp(PIstarw)^(-eta*gam)*exp(ld)^(gam))*...
    ((exp(vp)^(-1)*((tempvaromega*exp(ld))^alf)*(exp(ld)^(1-alf)))-del*(tempvaromega*exp(ld)))^(-1),0.25,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
disp(ld)
% ((1-h*mu_z^(-1))^(-1)-bet*h*(mu_z-h)^(-1))
% [ld,fval,exitflag]=fsolve(@(ld)(1-bet*tetw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)))/(1-bet*tetw*mu_z^(eta*(1+gam))*PI^(eta*(1-chiw)*(1+gam)))...
%     -(eta-1)/eta*wstar/(psy*PIstarw^(-eta*gam)*ld^gam)*((1-h*mu_z^(-1))^(-1)-bet*h*(mu_z-h)^(-1))*...
%     ((mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alf-tempvaromega*(1-(1-del)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi)^(-1),0.25,options);

%added back in l=vw*ld
l=log(exp(vw)*exp(ld));
k=log(tempvaromega*exp(ld)); 
x=log(del*exp(k));
yd=log((exp(k)^alf*exp(ld)^(1-alf))/exp(vp));

%mu_yd=exp(Lamyd);
%adding in ydbar for Taylor Rule
%ydbar=1*yd;

c=log((exp(vp)^(-1)*((tempvaromega*exp(ld))^alf)*(exp(ld)^(1-alf)))-del*(tempvaromega*exp(ld)));
lam=log(1/exp(c));

f=log((eta-1)/eta*exp(wstar)*exp(PIstarw)^(-eta)*exp(lam)*exp(ld)/(1-bet*tetw*exp(PI)^(-(1-chiw)*(1-eta))));
%f=log(psy*exp(PIstarw)^(-eta*(1+gam))*exp(ld)^(1+gam)/(1-bet*tetw*(exp(PI)^(eta*(1-chiw)*(1+gam)))));
%the above equation is the original but not sure if it's right?
%f2=log(psy*exp(d)*exp(phi)*exp(PIstarw)^(-eta*(1+gam))*exp(ld)^(1+gam)/(1-bet*tetw*(exp(PI)^chiw/exp(PI))^(-eta*(1+gam))));
%f=log((eta-1)/eta*exp(wstar)*exp(PIstarw)^(-eta)*exp(lam)*exp(ld)/(1-bet*tetw*exp(PI)^(-(1-chiw)*(1-eta))));
f2=log(psy*d*phi*exp(PIstarw)^(-eta*(1+gam))*exp(ld)^(1+gam)/(1-bet*tetw*(exp(PI)^chiw/exp(PI))^(-eta*(1+gam))*(exp(wstar)/exp(wstar))^(eta*(1+gam))));

g1=log(exp(lam)*exp(mc)*exp(yd)/(1-bet*tetp*exp(PI)^((1-chi)*ela)));
g2=log(ela/(ela-1)*exp(g1));

% %Measurement Equations
% %yd_obs=yd;
% %c_obs=c;
yd_obs=0;
c_obs=0;
R_obs=R;
PI_obs=PI;
l_obs=0;
% PI_obs=0;
% %R_obs=R-Rbar;
% %//PI_obs=(1+PI)-(1+PIbar);
% l_obs=l;

%% end own model equations
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