function [ys,check] = test_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_ options_

% 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
betta = 1/(beta_new/100+1);
gammma = gamma_new/100;
PI_ss = PI_new/100+1;
tau_ss = tau_new/100*4;
gov_ss = gov_new/100;
s_b_ss = s_b_new/100*4;

r_ss = exp(gammma)/betta-(1-delta);
gamma1 = r_ss;
gamma2 = psi_uticost/(1-psi_uticost)*gamma1;
R_ss = PI_ss*exp(gammma)/betta;
Psi_ss = (epsilon_ss-1)/epsilon_ss;
wtilde_ss = (Psi_ss*(1-alppha)^(1-alppha)*alppha^alppha*r_ss^(-alppha))^(1/(1-alppha));
ktilde_L_ss = alppha/(1-alppha)*wtilde_ss/r_ss*exp(gammma);
ydtilde_L_ss = exp(-alppha*gammma)*ktilde_L_ss^alppha;
xtilde_L_ss = (1-(1-delta)*exp(-gammma))*ktilde_L_ss;
ctilde_L_ss = (1-gov_ss)*ydtilde_L_ss-xtilde_L_ss;
lambda_L_ss = 1/ctilde_L_ss*(exp(gammma)-h*betta)/(exp(gammma)-h);
s_ss = s_b_ss*(1-1/betta)+tau_ss-gov_ss;

c = 0;           %//1. consumption
x = 0;           %//2. investment
k = 0;           %//3. capital
yd = 0;          %//4. aggregate output
Psi = 0;         %//5. marginal costs
l = 0;           %//6. aggregate labor bundle
ld = 0;          %//7. aggregate labor demand
lambda = 0;      %//9. Lagrange multiplier
q = 0;           %//10. Tobin's marginal q
v = 0;           %//11. capacity utilization
r = 0;           %//12. rental rate of capital
wstar = 0;       %//13. optimal real wage
w = 0;           %//14. real wage
DeltaW = 0;      %//15. wage dispersion term
PIstarw = 0;     %//16. optimal wage inflation
PIstar = 0;      %//17. optimal price inflation
PI = 0;          %//18. inflation
R = 0;           %//19. nominal Interest rate
z_shock = 0;     %//21. shock to trend growth rate of the economy
mu_shock = 0;    %//22. shock to growth rate of investment-specific technology growth
Ud = 0;          %//23. preference shock
eta_shock = 0;
epsilon_shock = 0;
Utau = 0;
Ugov = 0;         %//25. transformed government spending shock
tau = 0;
gov = 0;
s_b = 0;
s = 0;
dc = 100*gammma;             %consumption growth rate (observable 1)
dinve = 100*gammma;          %investment growth rate (observable 2)
labobs = 0;                  %labor hours (observable 4)
pinfobs = 100*log(PI_ss);    %inflation rate (observable 5)
dw = 100*gammma;             %real wage growth rate (observable 6)
robs = 100*log(R_ss);        %interest rate (observable 7)
gshare = log(gov_ss);
bshare = log(s_b_ss);
tshare = log(tau_ss);
%% 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