function [ys,check] = GSAmodel_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,:)); % get parameter name
  eval([ paramname ' = M_.params(' int2str(ii) ');']); % get parameter value
end


% initialize indicator
check = 0;


%% Enter model equations here

options=optimset('MaxFunEvals',400,'MaxIter',10000,'TolFun',1e-5); % set options for numerical solver

% varialbes with closed form solution
del_P=1; % price dispersion home
del_P_star=1; % price dispersion foreign
del_P_tilde=1;
del_P_tilde_star=1;
del_e=1;
P_H_opt_P_H=1; % home opt price / home price
P_F_opt_star_P_F_star=1; % foreign opt price / foreign price
PI=1; % home inflation
PI_star=1; % foreign inflation
PI_H=1; % home producer inflation
PI_F_star=1; % foreign producer inflation
t_r=0; % transfer due to porfolio adjustment costs
t_r_star=0;
markup = (muu/(muu-1));
R=1/bet;
R_star=1/bet;

b_F=o_bar;
b_H_star=o_bar_star;

b_H=(1-n)/n*b_H_star;
b_F_star=n/(1-n)*b_F;

% variables to be solved numerically
% names of variables whose ss have to be solved numerically
global names_x0 
names_x0=strvcat( 'c', 'c_star', 'P_H_P', 'P_F_star_P_star'); %n has to be the last
% read the initial guess for these variables from ys
x0=NaN(size(names_x0,1),1);
for jkl = 1:size(names_x0, 1)
  varname_ = deblank(names_x0(jkl,:));
  posjunk = loc(M_.endo_names, varname_); % find the location of the variable in ys
  x0(jkl)= ys(posjunk); % get initial guess for variables from initval bloc
end

[ss_vec,fval,exitflag ]=fsolve(@find_ss,x0,options); % ss_vec returns the endog. ss vector
[ss_vec,fval,exitflag ]=fsolve(@find_ss,real(ss_vec),options); % ss_vec returns the endog. ss vector

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

% set variables to calculated SS values
for jkl = 1:size(names_x0, 1)
  varname_ = deblank(names_x0(jkl,:)); % get name of variable
  eval([varname_ ' =  ss_vec(' int2str(jkl) ');']); % set variable to calculated steady-state values value
end

lamb = c^(-sig_c);
lamb_star = c_star^(-sig_c);

F = 1/(1-thet*bet)*lamb*P_H_P^(1-eps)*(n*c+(1-n)*c_star);
F_star = 1/(1-thet_star*bet)*lamb_star*P_F_star_P_star^(1-eps)*(n*c+(1-n)*c_star);
K=F;
K_star=F_star;

eps_g =0;
eps_g_star = 0;
eps_t=0;
eps_t_star=0;
eps_R =0;
eps_R_star =0;
A = A_ss;
A_star = A_ss_star;
i=R-1;
i_star=R_star-1;

y=P_H_P^(-eps)*((1-alph)*n*c+alph_star*(1-n)*c_star)/n;
y_star=P_F_star_P_star^(-eps)*(alph*n*c+(1-alph_star)*(1-n)*c_star)/(1-n);



% update parameters (those for which this ss file solves)
R_ss = R;
R_ss_star = R_star;
PI_ss=PI;
PI_ss_star=PI_star;
y_ss = y;
y_ss_star=y_star;
%% 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

% ---------------------------------------------------------------------------------------------------
% ---------------------------------------------------------------------------------------------------



function [fval] = find_ss(x0)

global M_ names_x0

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:)); % get parameter name
  eval([ paramname ' = M_.params(' int2str(ii) ');']); % set param to parameter value
end

% read out initial guesses with their names
for jkl = 1:size(names_x0, 1)
  varname_ = deblank(names_x0(jkl,:)); % get name of variable
  eval([varname_ ' =  x0 (' int2str(jkl) ');']); % set variable to initial guess value
end

markup = (muu/(muu-1));
%markup_w = (etaa/(etaa-1));
R=1/bet;
R_star=1/bet;
b_F=o_bar;
b_H_star=o_bar_star;
b_H=(1-n)/n*b_H_star;
b_F_star=n/(1-n)*b_F;

% -- evaluate ss

jkl=0; clear fval; 
jkl=jkl+1;
fval(jkl,1)= c-(1-bet)*(b_H - b_F)-P_H_P^(1-eps)*(n*c+(1-n)*c_star);     
jkl=jkl+1;
fval(jkl,1)= 1 - (1-alph)*P_H_P^(1-eps) - alph_star*P_F_star_P_star^(1-eps);     
jkl=jkl+1;
fval(jkl,1)= -P_H_P^(-eps*sig_N)*((n*c+(1-n)*c_star)/A_ss)^(1+sig_N) + P_H_P*(n*c+(1-n)*c_star)*c^(-sig_c);
jkl=jkl+1;
fval(jkl,1)= -P_F_star_P_star^(-eps*sig_N)*((n*c+(1-n)*c_star)/A_ss_star)^(1+sig_N) + P_F_star_P_star*(n*c+(1-n)*c_star)*c_star^(-sig_c);