function [ys,params,check] = Anzoategui_2017_Simple_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% 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
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)
 
global M_
    
%alph = [];
%bet  = [];
%gamma = []; % these 3 are needed so they don't clash with functions with the same name - in the future better to use different parameter names: bet, alph, gam

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

% initialize indicator (check=failure indicator -- indicates ss computation
% failed)
check = 0;



%%%% ---------------------- %%%    
%% Solve for Steady State %%%%%
%%% ------------------------%%%

%% Follow own derivations


% Gross interest rate
r = (1 + GAMMA_Y)/BET;

% HHs Stoch. Discount factor
delta = BET;

% Rental rate of capital
d = r - 1;

% Adopted technologies (normalized to one)
ad = a_ss;

% Marginal costs final goods producer
mc = 1/ad^(THETA - 1);

kd_yd_ss = ALPHA*mc/(VARSIGMA*(r - 1 + DELTA));

wdl_yd = (1 - ALPHA)*mc/VARSIGMA;

id_yd = ((GAMMA_Y + DELTA)/(1 + GAMMA_Y))*kd_yd_ss;

cd_yd = 1 - id_yd;

% Exo. process
muw = 1;

% Unskilled Hours Worked
l = l_ss;

% Output
yd = kd_yd_ss^(ALPHA/(1 - ALPHA));

% Capital
kd = kd_yd_ss^(1/(1 - ALPHA));

% Investment
id = (GAMMA_Y + DELTA)/(1 + GAMMA_Y)*kd_yd_ss^(1/(1 - ALPHA));

% Consumption
cd = yd - id;

% Wage rate (unskilled)
wd = (1 - ALPHA)*mc/VARSIGMA*yd;

NU = wd*cd*1/(l^(RHO)*muw);


% (scaled) Value of adopted tech
vad = (1 + GAMMA_A)/(1 + GAMMA_A - PHI*BET)*((VARSIGMA - 1)/VARSIGMA)*yd*mc;

% Adoption probability
lambda = lambda_ss;

% Unadopted technologies
zd = (1 + GAMMA_A - PHI + lambda*PHI)/(lambda*PHI);


% (Scaled) value of unadopted tech.
jzd = ((1 - RHO_LAMBDA)*lambda*PHI*BET)/(1 + GAMMA_A - (1 - lambda + RHO_LAMBDA*lambda)*PHI*BET)*vad*zd;

lsaz_wsd = RHO_LAMBDA*lambda*PHI*BET/(1 + GAMMA_A)*(vad*zd - jzd);

lsr_wsd = BET*jzd/(1 + GAMMA_A)*(1 + GAMMA_A - PHI);

ls_wsd = (1 - 1/zd)*lsaz_wsd + lsr_wsd;

% Skilled Hours worked
ls = ls_ss;

NU_S = ls_wsd*cd*1/muw;


% Wage rate skilled
wsd = (1 - 1/zd)*lsaz_wsd + lsr_wsd;


% Skilled labor in adoption
lsaz = lsaz_wsd/wsd;

% Skilled labor in R&D
lsr = lsr_wsd/wsd;

%% Back out KAPPA_LAMBDA and steady state chi
KAPPA_LAMBDA = lambda/(lsaz)^(RHO_LAMBDA);

chi = (1 + GAMMA_A - PHI)/(lsr^RHO_Z);

chi_ss = log(chi);

% =============================== %
% Map to logs %

r = log(r);
delta = log(delta);
d = log(d);
mc = log(mc);
yd = log(yd);
kd = log(kd);
id = log(id);
cd = log(cd);
wd = log(wd);
vad = log(vad);
zd = log(zd);
jzd = log(jzd);
wsd = log(wsd);
lambda = log(lambda);
l = log(l);
ls = log(ls);
lsaz = log(lsaz);
lsr = log(lsr);
ad = log(ad);

muw = log(muw);
chi = log(chi);



%%%%---------------------------%%%%
%% End Own Steady State Solution %%
%%%%---------------------------%%%%


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

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end

    
  
end


