function [ys,check] = myexample_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
  
mutilde_ss = exp(mu_new/100); %steady state gross growth rate of investmentment-specific technology
Atilde_ss = exp(A_new/100); %steady state gross growth rate of neutral technology
PI_ss = PI_new/100+1; %steady state gross inflation rate

ztilde_ss = exp((log(Atilde_ss)+alppha*log(mutilde_ss))/(1-alppha)); %s.s. growth rate of the economy
rtilde_ss = ztilde_ss*mutilde_ss/betta-(1-delta); %s.s. real rental rate of capital (gamma1)
gammma1 = rtilde_ss;
R_ss = PI_ss*ztilde_ss/betta; %s.s. nominal interest rate
PIstar_ss = ((1-theta_p*PI_ss^((1-chi)*(epsilon-1)))/(1-theta_p))^(1/(1-epsilon)); %s.s. ratio of optimal price to aggregate price
Psi_ss = ((epsilon-1)/epsilon)*PIstar_ss*(1-betta*theta_p*PI_ss^((1-chi)*epsilon))/(1-betta*theta_p*PI_ss^((1-chi)*(epsilon-1))); %s.s. real marginal cost
PIstarw_ss = ((1-theta_w*PI_ss^((1-chi_w)*(eta-1))*ztilde_ss^(eta-1))/(1-theta_w))^(1/(1-eta)); %s.s. ratio of optimal wage to aggregate wage
wtilde_ss = (1-alppha)*(Psi_ss*(alppha/rtilde_ss)^alppha)^(1/(1-alppha)); %s.s. real aggregate wage
wstartilde_ss = wtilde_ss*PIstarw_ss; %s.s. real optimal wage
DeltaP_ss = (1-theta_p)*PIstar_ss^(-epsilon)/(1-theta_p*PI_ss^((1-chi)*epsilon)); %s.s. relative price dispersion
DeltaW_ss = (1-theta_w)*PIstarw_ss^(-eta)/(1-theta_w*PI_ss^((1-chi_w)*eta)*ztilde_ss^eta); %s.s. relative wage dispersion
OMEGA_ss = (alppha/(1-alppha))*wtilde_ss*ztilde_ss*mutilde_ss/rtilde_ss; %s.s. capital-labor ratio

% Solve for s.s. L_d
LHS = (1-betta*theta_w*ztilde_ss^(eta*(1+gamma_l))*PI_ss^(eta*(1-chi_w)*(1+gamma_l)))/...
      (1-betta*theta_w*ztilde_ss^(eta-1)*PI_ss^((1-chi_w)*(eta-1)));
coeff1 = Atilde_ss*OMEGA_ss^alppha/(ztilde_ss*DeltaP_ss)-(ztilde_ss*mutilde_ss-(1-delta))*OMEGA_ss/(ztilde_ss*mutilde_ss);
x0 = 0.5; %initial guess of L_d
%options = optimset('Display','off');   % Option to display output
[ld_ss,Fval,exitflag] = fsolve(@(x)LHS-varpsi*PIstarw_ss^(-eta*gamma_l)*x^gamma_l/...
    (((eta-1)/eta)*wstartilde_ss*(1-h/ztilde_ss)^(-1)*(coeff1*x-Phi/DeltaP_ss)^(-1)),x0,options); %solve for s.s. labor demand
%To check if fsolve converges to a solution
if exitflag ~= 1 
   check=1; %set failure indicator
   return; %return without updating steady states
end  
    
ktilde_ss = OMEGA_ss*ld_ss; %s.s. capital stock
ydtilde_ss = (Atilde_ss*ktilde_ss^alppha*ld_ss^(1-alppha)/ztilde_ss-Phi)/DeltaP_ss; %s.s. output
xtilde_ss = (ztilde_ss*mutilde_ss-(1-delta))*ktilde_ss/(ztilde_ss*mutilde_ss); %s.s. investment
ctilde_ss = coeff1*ld_ss-Phi/DeltaP_ss; %s.s. consumption
lambdatilde_ss = (1-h/ztilde_ss)^(-1)*ctilde_ss^(-1); %s.s. shadow price
l_ss = DeltaW_ss*ld_ss; %s.s. labor supply
numtilde_ss = lambdatilde_ss*(1-betta/(ztilde_ss*PI_ss)); %s.s. ratio related to money demand

omega_ss = ztilde_ss*PI_ss; %steady state money growth rate

ss = [ztilde_ss,...         %1
      rtilde_ss,...         %2
      R_ss,...              %3
      PIstar_ss,...         %4
      Psi_ss,...            %5
      PIstarw_ss,...        %6
      wtilde_ss,...         %7
      wstartilde_ss,...     %8
      DeltaP_ss,...         %9
      DeltaW_ss,...         %10
      OMEGA_ss,...          %11
      ld_ss,...             %12
      ktilde_ss,...         %13
      ydtilde_ss,...        %14
      xtilde_ss,...         %15
      ctilde_ss,...         %16
      lambdatilde_ss,...    %17
      l_ss,...              %18
      numtilde_ss,...       %19
      gammma1,...           %20
      mutilde_ss,...        %21
      Atilde_ss,...         %22
      PI_ss,...             %23
      omega_ss];            %24

%To check if s.s. values are complex 
N_ss = size(ss,2);
flag_clxss = zeros(1,N_ss); 
for i = 1:N_ss
   if imag(ss(i)) ~= 0
      flag_clxss(i) = 1;
   end
end
invalid_ss = find(flag_clxss,1);
if ~isempty(invalid_ss)
    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
  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
