function [ys,params,check] = NK_baseline_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)

% 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 = 0;


%% Enter model equations here

options=optimset(); % set options for numerical solver

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Structural Parameters (quarterly calibration)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Steady State
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Shock processes

% inflation targeting
pi=pibar;   

% Foreign inflation
pistar = pibar;
%pistar = log(1.005);

% tfp
a=log(1);

% Price markup shock
m=log(1);


rstar = 1/betta;

% share of public spending in ss
g=gss;     

ystar=log(1);

% Marginal effiency of investment
z=log(1);

% Shock on discount factor
nu=log(1);

% Shock on disutility of labor
kappa=log(1);

% domestic inflation                     
piH=pi;





% foreign debt
d=dss;                       


% domestic price
pH=log(1);                        

% real FX rate
Q=log(1);                         

% hours of work/labor
n=log(0.33);                       

% real interest rate
rr=log(1/betta);                   

% nominal interest rate
i=iss;                   

% marginal value of investment (in terms of lambda) 
q=log(1);                          

% rental rate of capital
r=log(1/betta-(1-delta)); 

% aggregate prices
piHsharp=log(((exp(piH)^(1-epsilonp)-phip*exp(piH)^(zetap*(1-epsilonp)))/(1-phip))^(1/(1-epsilonp)));

% marginal cost
mc=log(((epsilonp-1)/epsilonp)*(exp(piHsharp)/exp(pi))*((1-phip*betta*exp(pi)^((1-zetap)*epsilonp))/(1-phip*betta*exp(pi)^((zetap-1)*(1-epsilonp))))); 

% Price dispersion
%nup = a*k^allpha*n^(1-allpha);
nup = log(((exp(pi)^(epsilonp)*(1-phip)*exp(piHsharp)^(-epsilonp))/(1-phip*exp(pi)^((1-zetap)*epsilonp))));

% capital
k=log(((allpha*exp(mc)*exp(a)*exp(n)^(1-allpha))/(exp(nup)*exp(r)))^(1/(1-allpha))); 

% domestic output
yH=(log((exp(a)*exp(k)^(allpha)*exp(n)^(1-allpha))/exp(nup)));                     

       

% real wage
w=log((1-allpha)*exp(mc)*exp(yH)/exp(n));         

% investment
inv=log(delta*exp(k));                   

% consumption
c=log(exp(yH)-exp(inv)-exp(g)-d*(1/betta-1));       

% marginal utility of consumption
lambda = log((exp(c)-b*exp(c))^(-sigma) - betta*b*(exp(c)-b*exp(c))^(-sigma));

% aggregate wage
%wsharp = (((1-phiw*pi^((zetaw-1)*(1-epsilonw)))/(1-phiw))^(1/(1-epsilonw)))*w; 
wsharp = log(((1-phiw*exp(pi)^((zetaw-1)*(1-epsilonw)))/(1-phiw))^(1/(1-epsilonw)) *exp(w));


% labor preference parameter
%psi= lambda*((epsilonw-1)/epsilonw)*(wsharp^(2+phi)/(w^(1+phi)*n^phi))*((1-phiw*betta*pi^((1-zetaw)*epsilonw*(1+phi)))/(1-phiw*betta*pi^((zetaw-1)*(1-epsilonw))));
psi= ((1-phiw*exp(pi)^((zetaw-1)*(1-epsilonw)))/(1-phiw))^((1+epsilonw*phi)/(1-epsilonw)) *exp(w) *exp(lambda)*exp(n)^(-phi)*((epsilonw-1)/epsilonw)*((1-phiw*betta*exp(pi)^((1-zetaw)*epsilonw*(1+phi)))/(1-phiw*betta*exp(pi)^((zetaw-1)*(1-epsilonw))));



tot=log(1);

% foreign parameter
gammastar=(exp(yH)-(1-gama)*(exp(c)+exp(inv))-exp(g))/(exp(ystar)*exp(tot)^(eta));   



% Auxiliary variable wage setting
%f1 = ((psi*(w/wsharp)^(epsilonw*(1+phi))*n^(1+phi))/(1-phiw*betta*pi^((1-zetaw)*epsilonw*(1+phi))));
f1 = log(((psi*exp(w)^(epsilonw*(1+phi))*exp(n)^(1+phi))/(1-phiw*betta*exp(pi)^((1-zetaw)*epsilonw*(1+phi)))));

% Auxiliary variable wage setting
%f2 = ((lambda*(w/wsharp)^(epsilonw)*n)/(1-phiw*betta*pi^((zetaw-1)*(1-epsilonw))));
f2 = log(((exp(lambda)*exp(w)^(epsilonw)*exp(n))/(1-phiw*betta*exp(pi)^((zetaw-1)*(1-epsilonw)))));

% Auxiliary variable wage setting
x1 = log(((exp(lambda)*exp(mc)*exp(yH))/(1-phip*betta*exp(piH)^((1-zetap)*epsilonp))));

% Auxiliary variable wage setting
x2 = log(((exp(lambda)*exp(yH))/(1-phip*betta*exp(piH)^(zetap*(1-epsilonp))*exp(pi)^(-1)*exp(piH)^(epsilonp))));

%De=log(exp(pi)/exp(pistar));
De = log(1);
ex=log(exp(pH)*gammastar*exp(ystar)*(exp(tot))^eta);

im=log(gama*exp(Q)^(1-eta)*(exp(c)+exp(inv)));

tb=log(exp(ex)-exp(im)+exp(Q)*kappaD/2*(d-dss)^2);

u=log(1);

gdp=log(exp(pH)*exp(yH));

hatk=k;     

proof=0;

gdp_obs = 0;
c_obs =0;
n_obs=0;
pi_obs=pibar;
i_obs= iss;
inv_obs=0;
ystar_obs=log(1);
rstar_obs = 1/betta;
pistar_obs = pibar;
tot_obs = 0;

%% end own model equations

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
