function [ys,params, check] = main_steadystate(ys,exo,M_, options_)

%initialize constraint parameter for passing through
occbin_borrcon_bind=M_.params(strcmp('occbin_borrcon_bind',M_.param_names));

% Pick up the parameters defined in param_core.m
param_core

% all shocks in the steady state are equal to 0
a = 1;
z_cb = 1;
u_g = 1;
u_tr = 1;

pi = 0.0196;
pi_im = pi;

R = 1 / bbeta;
R_foreign = 1 / bbeta;

exch = 1/1.08;
mc = p_SS * ((eeta - 1)/eeta);
p_im = (eeta_im / (eeta_im - 1)) * p_foreign_SS * exch;

%p_c = (1 / ((aalpha_c)^(aalpha_c) * (1 - aalpha_c)^(1 - aalpha_c))) * p_SS^(aalpha_c) * p_im^(1 - aalpha_c);
p_I = (1 / ((aalpha_I)^(aalpha_I) * (1 - aalpha_I)^(1 - aalpha_I))) * p_SS^(aalpha_I) * p_im^(1 - aalpha_I);
p_Ex = (1 / ((aalpha_Ex)^(aalpha_Ex) * (1 - aalpha_Ex)^(1 - aalpha_Ex))) * p_SS^(aalpha_Ex) * p_im^(1 - aalpha_Ex);


mup = (ddelta / u_constr) - 1;
llambda_K = mup + 1;
K_L = ((aalpha * p_SS * a) / ((llambda_K / bbeta) - llambda_K * ((1 - ddelta) + mup * u_constr)))^(1 / (1 - aalpha));

w = (1 - aalpha) * a * mc * (K_L)^aalpha;

x_0 = [0.1, 1];
options=optimset('Display', 'off', 'Algorithm', 'levenberg-marquardt'); % set options for numerical solver
fun = @(x) find_l(x, K_L, aalpha, aalpha_I, aalpha_Ex, ddelta, p_SS, p_im, p_I, p_Ex, thetta, w, pphi, cchi, a, nu_Im, nu_Ex);
xx = fsolve(fun,x_0,options);
l_0 = real(xx(1));
aalpha_c = real(xx(2));

k = K_L * l_0;
y = a * k^aalpha * (l_0)^(1 - aalpha);
I = ddelta * k;

p_c = (1 / ((aalpha_c)^(aalpha_c) * (1 - aalpha_c)^(1 - aalpha_c))) * p_SS^(aalpha_c) * p_im^(1 - aalpha_c);

Im    = nu_Im * p_SS * y / p_im;
y_im = Im;
Ex  = nu_Ex * p_SS * y / p_Ex;

c = (y * p_SS + Im * p_im - (I * p_I + Ex * p_Ex)) / p_c;


llambda = ((c * (1 - thetta))^(-cchi)) / p_c;
ll = (llambda * w)^(1 / pphi);

c_im = ((1 - aalpha_c) * p_c * c) / p_im;
c_d = (aalpha_c * p_c * c) / p_SS;
I_im = ((1 - aalpha_I) * p_I * I) / p_im;
I_d = (aalpha_I * p_I * I) / p_SS;
Ex_im = ((1 - aalpha_Ex) * p_Ex * Ex) / p_im;
Ex_d = (aalpha_Ex * p_Ex * Ex) / p_SS;

B_foreign = nu_B_foreign * (p_SS * y) / exch;
%Ex = (B_foreign * (R_foreign - 1) + Im * p_foreign_SS) * (exch /p_Ex);
B_foreign = (Im * p_foreign_SS - (Ex * p_Ex / exch)) / (1 - R_foreign);


Ex_SS = Ex * (1 / exch)^(eeta_Ex);
loan = llambda_K * k;
r_real = R - pi;


p_foreign_SS = log(p_foreign_SS);
y_im = log(y_im);
B_foreign = log(B_foreign);
Ex_SS = log(Ex_SS);
Ex = log(Ex);
Im = log(Im);
c_im = log(c_im);
c_d = log(c_d);
I_im = log(I_im);
I_d = log(I_d);
Ex_im = log(Ex_im);
Ex_d = log(Ex_d);
loan = log(loan);
p_im = log(p_im);
p_c = log(p_c);
p_I = log(p_I);
p_Ex = log(p_Ex);
a = log(a);
u_tr = log(u_tr);
%z_cb = log(z_cb);
c = log(c);
y = log(y);
I = log(I);
k = log(k);
mc = log(mc);
w = log(w);
ll = log(ll);
p_SS = log(p_SS);


%% Leave it unchanged

check = 0;
options=optimset();

% export parameters
nparams = size(M_.param_names,1);

params=NaN(nparams,1);
for icount = 1:nparams
    eval(['params(' num2str(icount) ') = ',M_.param_names{icount},';'])
end

nvars = M_.endo_nbr;

ys = zeros(nvars,1);

for i_indx = 1:nvars
   eval(['ys(i_indx)=',M_.endo_names{i_indx},';']) 
end
