function [ys,check] =BayesianEstimation0311_steadystate(ys,exe)
% 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;



interest=(1/betta)-1;
X=exp(x_ss);
x=x_ss;
T=exp(tau_ss);
tau=tau_ss;
A1=1;
A2=1;
a1=0;
a2=0;
vp=1;
Pi=0; 
PiStar=0; 
%C2=Rbar; 
%Y2=Rbar;
%D2 =C2/delta_g;

N=Nbar;
w=(nu-1)/nu;
mc=(nu-1)/nu;
% N1=N*0.98;
% To calibrate the the N1/N=0.98, I did the fsolve and ge the nu.
options = optimset('TolFun',10e-99,'TolX',10e-99,'Display','iter');
C1=fsolve(@(C1) C1-(gamma_L*(1-N)^(-eta)/w)^(-1/sigmma),0.9,options);
%gamma_L=fsolve(@(gamma_L) 0.91391-(gamma_L*(1-N)^(-eta)/w)^(-1/sigmma),0.9,options)
%gamma_L is calibrated so that C1=N1=Y1=mean(y_obs)=0.027957
%nu=fsolve(@(nu) w-(nu-1)/nu, 1.5,options);
%nu=1.5884
N1=C1;
N2=N-N1;
m=((1-betta)*C1^(-sigmma)/gamma_M)^(-1/delta); 
%gamma_M=fsolve(@(gamma_M) 0.12342-((1-betta)*C1^(-sigmma)/gamma_M)^(-1/delta),0.5);
% to match the data moment mean(M2_obs)=2.3
Y1=N1;
% pg=((kappa*(Y2/delta_g)^(-sigmma))/(C1^(-sigmma)));
% R=(Y2^(1/(1-alpha2)))/(N2^(alpha2/(1-alpha2)));
% i*pg-((i*w*N2)/(alpha2*Y2))=((1-alpha2)*w*N2)/(alpha2*R);
% Y2=fsolve(@Y2 (kappa*(Y2/delta_g)^(-sigmma))/(C1^(-sigmma))+((w*N2)/(alpha2*Y2))*((Y2*(1-alpha2)/((Y^(1/(1-alpha2)))/(N2^(alpha2/(1-alpha2)))))-1)-(1+i)*(((kappa*(Y2/delta_g)^(-sigmma))/(C1^(-sigmma)))-((w*N2)/(alpha2*Y2)))

%Y2=fsolve(@(Y2) (kappa*(Y2/delta_g)^(-sigmma))/(C1^(-sigmma))+((w*N2)/(alpha2*Y2))*((Y2*(1-alpha2)/((Y2^(1/(1-alpha2)))/(N2^(alpha2/(1-alpha2)))))-1)-(1+interest)*(((kappa*(Y2/delta_g)^(-sigmma))/(C1^(-sigmma)))-((w*N2)/(alpha2*Y2))),0.0000045,options)
%interest*((kappa*(Y2/delta_g)^(-sigmma))/(C1^(-sigmma)))-interest*w*N2/(alpha2*Y2)-((1-alpha2)*w*N2)*(N2^(alpha2/(1-alpha2)))/(alpha2*(Y2^(1/(1-alpha2))))
%Y2=fsolve(@(Y2) interest*((kappa*((Y2+IM)/delta_g)^(-sigmma))/(C1^(-sigmma)))-interest*w*N2/(alpha2*Y2)-(1-alpha2)*w*N2^(1/(1-alpha2))/(alpha2*(Y2^(1/(1-alpha2)))),0.0027417,options)
%Y2=fsolve(@(Y2) interest*((kappa*(Y2/delta_g)^(-sigmma))/(C1^(-sigmma)))-interest*w*N2/(alpha2*Y2)-(1-alpha2)*w*N2^(1/(1-alpha2))/(alpha2*(Y2^(1/(1-alpha2)))),0.0027417,options);
%delta_g=fsolve(@(delta_g) interest*((kappa*(0.003068067*C1/delta_g)^(-sigmma))/(C1^(-sigmma)))-interest*w*N2/(alpha2*Y2)-(1-alpha2)*w*N2^(1/(1-alpha2))/(alpha2*(Y2^(1/(1-alpha2)))),0.0027417,options);
Y2=fsolve(@(Y2) interest*((kappa*((Y2)/delta_g)^(-sigmma))/(C1^(-sigmma)))-interest*w*N2/(alpha2*Y2)-(1-alpha2)*w*N2^(1/(1-alpha2))/(alpha2*(Y2^(1/(1-alpha2)))),0.0027417,options);

%EX=fsolve(@(EX) interest*((kappa*((0.003*C1-EX)/delta_g)^(-sigmma))/(C1^(-sigmma)))-interest*w*N2/(alpha2*0.003*C1)-(1-alpha2)*w*N2^(1/(1-alpha2))/(alpha2*((0.003*C1)^(1/(1-alpha2)))),0.002,options)
% so EX_ss=0.0025577 net export
%x_ss=Y2
%{
 % Octave
 % x0 = [0.000005;40;0.000001];
 % options = optimset('TolFun',10e-99,'TolX',10e-99,'Display','iter');
 % [x,fval,info] = fsolve(@goldnonlinearproblemOctave,x0,options)
%}
%it doesnt really matter how big is the pg, cuz u wont know what is the unit of good C2
Y=Y1+Y2;
pg=((kappa*(Y2/delta_g)^(-sigmma))/(C1^(-sigmma)));
R=((Y2^(1/(1-alpha2)))/(N2^(alpha2/(1-alpha2))));
C2=Y2;
D2 =C2/delta_g;
G=m*T/pg;
x1=(C1^(-sigmma)*mc*Y1)/(1-phi*betta);
x2=(C1^(-sigmma)*Y1)/(1-phi*betta); 
Y_obs=Y;
PIY_obs=Pi;
rr_obs=interest;
M2_obs=m;
%% 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

