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

pi=log(piss); 

a=0;
n=log(0.333);

xi=log(epsy/(epsy-1));

if matlab_ver_less_than('9.1')
    options=optimset('TolFun',1e-10,'TolX',1e-10);
else
    options=optimoptions('fsolve','FunctionTolerance',1e-15,'OptimalityTolerance',1e-15);
end

[z,fval,exitflag]=fsolve(@(par_vec)solve_calibration(par_vec,kai,bbc,bbe,ddk,ddg,aa,ddb,bip,big,bg,bor,teta,car,xi,pi,n,csige,a,csig,etab,etad),[3.1736;-0.3513;0.2302;-0.9567;1.5030;0.7375;3.0901;-0.9434;0.6479;3.8957;-0.0566;1.5041;1.5041;1.8633;2.7081;-0.7985;-0.4620;-1.1590;-3.2702;-1.1552;-0.9163;0;1;1.2;1.45;0.05;0.045;0.045;0.06],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

b=z(1);
psi=z(2);
q=z(3);
ce=z(4);
lae=z(5);
cc=z(6);
d=z(7);
jb=z(8);
kb=z(9);
ke=z(10);
ip=z(11);
y=z(12);
ye=z(13);
w=z(14);
kg=z(15);
ig=z(16);
g=z(17);
lac=z(18);
rk=z(19);
or=z(20);
lv=z(21);
gae=z(22);
jr=z(23);
mkd=z(24);
mkb=z(25);
rd=z(26);
rrs=z(27);
rrb=z(28);
rb=z(29);

yobrrbs=0;
ipobs=0;
igobs=0;

blv=exp(lv);
bpsi=exp(psi);

o1=rd;
o2=rb; 
%% The conditions that need to be satisfied
residuals=NaN(24,1);
residuals(1)=b-log((exp(lv)*exp(pi)*exp(q)*exp(ke)*(1-ddk))/(1+exp(rb)));
residuals(2)=lae - log(exp(ce)^-csige);
residuals(3)=ce-log(exp(ye)/exp(xi)+exp(b)+(1-ddk)*exp(ke)*exp(q)-exp(w)*exp(n)-exp(ke)*exp(q)-(1+exp(rb))*exp(b)/exp(pi));
residuals(4)=cc-log(exp(y)-exp(ce)-exp(g)-exp(ig)-exp(ip));
residuals(5)=q-log(1/(1 - kai /2 * (exp(psi)-1) ^ 2 - kai*(exp(psi)- 1) * exp(psi) + bbc * kai * (exp(psi) -1) * exp(psi))) ;
residuals(6)=d-log(exp(b)-exp(kb));
residuals(7)= jb-log(exp(b)*exp(rb)-exp(rd)*exp(d));
residuals(8)= kb-log(((exp(rb)-exp(rd))*exp(b))/(exp(pi)-1+ddb-exp(rd)));
residuals(9)= w-log((1-aa)*exp(ye)/(exp(n)*exp(xi)));
residuals(10)= ig-log(big*exp(ye));
residuals(11)= kg-log(exp(ig)/ddg);
residuals(12)= ip-log(bip*exp(ye));
residuals(13)= or-log(bor*exp(y));
residuals(14)= ke-log((exp(ye)/(exp(a)*(exp(n)^(1-aa))*(exp(kg)^teta)))^(1/aa));
residuals(15)= psi-log(-(((2/kai)*(1-ddk*exp(ke)/exp(ip)))^(1/2))+1);
residuals(16)=rk-log((aa)*exp(ye)/(exp(ke)*exp(xi)));
residuals(17)= g-log(bg*exp(ye));
residuals(18)= lac - log(exp(cc)^-csig);
residuals(19)= jr - log(exp(y)*(1 - (1/exp(xi))));
residuals(20)=lv-log((exp(q)-bbe*(exp(q)*(1-ddk)+exp(rk)))/(exp(pi)*exp(q)*(1-ddk)*((1-bbe*(1+exp(rb))/exp(pi))/(1+exp(rb)))));
residuals(21)=gae-log((exp(lae)-exp(lae)*bbe*(1+exp(rb))/exp(pi))/(1+exp(rb)));
residuals(22)=car-exp(b)/exp(kb);
residuals(23)=y-ye;
residuals(24)=mkd-log(etad/(etad-1));
residuals(25)=mkb-log(etab/(etab-1));
residuals(26)=rd- log(exp(pi)/bbc - 1);
residuals(27)=rrs - log(exp(rd) /exp(mkd));
residuals(28)=rrb- rrs;
residuals(29)=rb-log(exp(rrb)*exp(mkb));
residuals(30)=(beke-log(bbc*exp(q)*exp(lv)*(1+exp(rd))*(1-ddk)/(1+exp(rb))));


if max(abs(residuals)>1e-6)
    check=1; %set failure indicator
    return; %return without updating steady states
end


%% end own model equations

for iter = 1:length(M_.params) %update parameters set in the file
  paramname = deblank(M_.param_names(iter,:));
  eval(['M_.params(' num2str(iter) ') = ' paramname ';']);
% 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