function [ys,check,param_update] = TNT_Big_Jan20_steadystate(ys,exe,params1)
% function [ys,check] = TNT_v10_options_steadystate(ys,exe)
%% Do note change this part, from here...
global M_ lgy_ 

if nargin>2
    M_.params=params1;
end

if isfield(M_,'param_nbr') == 1
NumberOfParameters = M_.param_nbr;
for iter = 1:NumberOfParameters
    paramname = deblank(M_.param_names(iter,:));
    eval([ paramname ' = M_.params(' int2str(iter) ');']);
end
check = 0;
end

% ..to here

%% Steady state equations
% Assigning calibrated staedy state values
Rstr_ss     = Rw_ss*xiR_ss;
aX_ss       = a_ss^((2*KGAM_X-1)/KGAM_X);
aCo_ss      = a_ss^((2*KGAM_CO-1)/KGAM_CO);
ppiNFE_ss   = BAR_PPI;
ppi_ss      = BAR_PPI;
BET         = a_ss^(SIG)*ppi_ss/R_ss;  
ppiS_ss     = a_ss^(SIG)*ppi_ss/(Rstr_ss*BET);
ppistr_ss   = ppi_ss/ppiS_ss;
ppiX_ss     = ppi_ss;
ppiM_ss     = ppi_ss;
ppiN_ss     = ppi_ss;
pXst_ss     = 1;
pMst_ss     = 1;
pNst_ss     = 1;
wXst_ss     = 1;
wNst_ss     = 1;
deltX_ss    = 1;
deltM_ss    = 1;
deltN_ss    = 1;
deltwX_ss   = 1;
deltwN_ss   = 1;
mcX_ss      = (EPS_X-1)/EPS_X;
mcM_ss      = (EPS_M-1)/EPS_M;
mcN_ss      = (EPS_N-1)/EPS_N;
hXd_ss      = hX_ss;
hNd_ss      = hN_ss;
muwX_ss     = (EPS_W-1)/EPS_W;
muwN_ss     = (EPS_W-1)/EPS_W;

%% Relative prices
pTI_pI_ss   = (pX_pI_ss)^(GAM_TI)*(pM_pI_ss)^(1-GAM_TI);
pN_pI_ss    = ((1-(1-GAM_I)*pTI_pI_ss^(1-VRHO_I))/GAM_I)^(1/(1-VRHO_I));
pT_pI_ss    = (pX_pI_ss)^GAM_T*(pM_pI_ss)^(1-GAM_T);

muX_pI_ss   = 1/u_ss;
muN_pI_ss   = 1/u_ss; 

RX_ss       = muX_pI_ss*(1-BET*a_ss^(-SIG)*(1-DELT))/(BET*a_ss^(-SIG)*pX_pI_ss);
RN_ss       = muN_pI_ss*(1-BET*a_ss^(-SIG)*(1-DELT))/(BET*a_ss^(-SIG)*pN_pI_ss); 

pm_pI_ss    = mcM_ss*pM_pI_ss;
rer_pI_ss   = pm_pI_ss/pmstr_ss;

pME_pI_ss   = pM_pI_ss^(1-GAM_EM)*pT_pI_ss^(GAM_EM);

mcvX_ss     = ((mcX_ss*pX_pI_ss^(1-GAM_X)*(1-GAM_X)^(1-GAM_X)*GAM_X^GAM_X)/pME_pI_ss^(1-GAM_X))^(1/GAM_X);
mcvN_ss     = ((mcN_ss*pN_pI_ss^(1-GAM_N)*(1-GAM_N)^(1-GAM_N)*GAM_N^GAM_N)/pME_pI_ss^(1-GAM_N))^(1/GAM_N); 

wN_pI_ss     = ((mcvN_ss*zN_ss*a_ss^(1-ALP_N)*pN_pI_ss*(1-ALP_N)^(1-ALP_N)*ALP_N^ALP_N)/(pN_pI_ss*RN_ss)^ALP_N)^(1/(1-ALP_N)); 
xih_N_pI_ss  = muwN_ss*wN_pI_ss/(hN_ss^(VARPHI));
xih_X_pI_ss  = xih_N_pI_ss;
wX_pI_ss     = (xih_X_pI_ss*hX_ss^(VARPHI))/muwX_ss;

zX_ss       = ((pX_pI_ss*RX_ss)^ALP_X*wX_pI_ss^(1-ALP_X))/(mcvX_ss*aX_ss^(1-ALP_X)*pX_pI_ss*(1-ALP_X)^(1-ALP_X)*ALP_X^ALP_X);

vX_ss       = hX_ss*zX_ss*aX_ss^(1-ALP_X)*(ALP_X*wX_pI_ss/((1-ALP_X)*pX_pI_ss*RX_ss))^ALP_X;
vN_ss       = hN_ss*zN_ss*a_ss^(1-ALP_N)*(ALP_N*wN_pI_ss/((1-ALP_N)*pN_pI_ss*RN_ss))^ALP_N;

kX_ss       = (a_ss*vX_ss/(zX_ss*aX_ss^(1-ALP_X)))*(ALP_X*wX_pI_ss/((1-ALP_X)*pX_pI_ss*RX_ss))^(1-ALP_X);
kN_ss       = (a_ss*vN_ss/(zN_ss*a_ss^(1-ALP_N)))*(ALP_N*wN_pI_ss/((1-ALP_N)*pN_pI_ss*RN_ss))^(1-ALP_N);

yX_ss       = vX_ss*(GAM_X*pME_pI_ss/((1-GAM_X)*pX_pI_ss*mcvX_ss))^(-(1-GAM_X));
yN_ss       = vN_ss*(GAM_N*pME_pI_ss/((1-GAM_N)*pN_pI_ss*mcvN_ss))^(-(1-GAM_N));

mX_ss       = yX_ss*(((1-GAM_X)*pX_pI_ss*mcvX_ss)/(GAM_X*pME_pI_ss))^(GAM_X);
mN_ss       = yN_ss*(((1-GAM_N)*pN_pI_ss*mcvN_ss)/(GAM_N*pME_pI_ss))^(GAM_N);

f1X_ss      = xiX_ss*(EPS_X-1)*yX_ss/(EPS_X*(1-BET*a_ss^(1-SIG)*THET_X));
f1N_ss      = xiN_ss*(EPS_N-1)*yN_ss/(EPS_N*(1-BET*a_ss^(1-SIG)*THET_N));

iX_ss       = kX_ss*(1-(1-DELT)/a_ss)/u_ss;
iN_ss       = kN_ss*(1-(1-DELT)/a_ss)/u_ss;

i_ss        = iX_ss + iN_ss;

itilN_ss    = GAM_I*pN_pI_ss^(-VRHO_I)*i_ss;
itilT_ss    = (1-GAM_I)*pTI_pI_ss^(-VRHO_I)*i_ss;
itilX_ss    = GAM_TI*(pX_pI_ss/pTI_pI_ss)^(-1)*itilT_ss;
itilM_ss    = (1-GAM_TI)*(pM_pI_ss/pTI_pI_ss)^(-1)*itilT_ss;

pYgdp_pI_ss = (pX_pI_ss*yX_ss + pN_pI_ss*yN_ss-pM_pI_ss*(mX_ss+mN_ss))/(1-sCo_ss-sM_ss*(pM_pI_ss-pm_pI_ss)/pm_pI_ss);

yCo_ss      = sCo_ss*pYgdp_pI_ss*a_ss/(rer_pI_ss*pCostr_ss*aCo_ss);

yM_ss       = sM_ss*pYgdp_pI_ss/(pm_pI_ss);
g_ss        = sG_ss*pYgdp_pI_ss/(pN_pI_ss);

f1M_ss      = xiM_ss*(EPS_M-1)*yM_ss/(EPS_M*(1-BET*a_ss^(1-SIG)*THET_M));

m_ss        = yM_ss;
cN_ss       = yN_ss - g_ss - itilN_ss;
cM_ss       = yM_ss - itilM_ss - mX_ss - mN_ss;
cT_ss       = (cM_ss/(1-GAM_T))*(pM_pI_ss/pT_pI_ss);
cX_ss       = GAM_T*(pX_pI_ss/pT_pI_ss)^(-1)*cT_ss;

GAM         = pN_pI_ss^(VRHO)*cN_ss/(pT_pI_ss^VRHO*cT_ss+pN_pI_ss^VRHO*cN_ss);

pNFE_pI_ss  = ((1-GAM)*pT_pI_ss^(1-VRHO)+GAM*pN_pI_ss^(1-VRHO))^(1/(1-VRHO));

pI_ss       = (pNFE_pI_ss^(1-GAM_FC-GAM_EC)*pT_pI_ss^(GAM_FC+GAM_EC))^(-1);

pX_ss       = pX_pI_ss*pI_ss;
pM_ss       = pM_pI_ss*pI_ss;
pN_ss       = pN_pI_ss*pI_ss;
pT_ss       = pT_pI_ss*pI_ss;
pTI_ss      = pTI_pI_ss*pI_ss;
pNFE_ss     = pNFE_pI_ss*pI_ss;
pME_ss      = pME_pI_ss*pI_ss;
rer_ss      = rer_pI_ss*pI_ss;
pm_ss       = pm_pI_ss*pI_ss;
muX_ss      = muX_pI_ss*pI_ss;
muN_ss      = muN_pI_ss*pI_ss;

wX_ss       = wX_pI_ss*pI_ss;
wN_ss       = wN_pI_ss*pI_ss;
xih_N_ss    = xih_N_pI_ss*pI_ss;
xih_X_ss    = xih_X_pI_ss*pI_ss;
 
c_ss        = pN_ss^VRHO*cN_ss/GAM;

cXstr_ss    = yX_ss - itilX_ss - cX_ss;
ystr_ss     = cXstr_ss*(pX_ss/rer_ss)^EPS_STR/xiXstr_ss;

gdp_ss      = c_ss + g_ss + i_ss + cXstr_ss + yCo_ss*aCo_ss/a_ss -m_ss;

pY_ss       = pI_ss*pYgdp_pI_ss/gdp_ss;
tb_ss       = rer_ss*pCostr_ss*yCo_ss*aCo_ss/a_ss+pX_ss*cXstr_ss-pm_ss*m_ss;

bstr_ss     = (tb_ss-(1-VARTHETA)*rer_ss*pCostr_ss*yCo_ss*aCo_ss/a_ss)/(rer_ss*(1-Rstr_ss/(ppistr_ss*a_ss)));
BAR_B       = bstr_ss*rer_ss/(pY_ss*gdp_ss);

la_ss       = xibet_ss*c_ss^(-SIG)*(1-PHI_C/a_ss)^(-SIG);

f1wX_ss     = (EPS_W-1)*la_ss*hXd_ss/(EPS_W*(1-BET*a_ss^(1-SIG)*THET_WX));
f1wN_ss     = (EPS_W-1)*la_ss*hNd_ss/(EPS_W*(1-BET*a_ss^(1-SIG)*THET_WN));

% EXTRA VARIABLES
tby_ss      = tb_ss/(pY_ss*gdp_ss);
g_aCo_ss    = a_ss;
expo_ss     = cXstr_ss+yCo_ss*aCo_ss/a_ss;
pNpT_ss      = pN_ss/pT_ss;
TICo_ss     = pCostr_ss/pmstr_ss;
TINoCo_ss   = pX_ss/(rer_ss*pmstr_ss);
TITot_ss    = (cXstr_ss*TINoCo_ss+yCo_ss*aCo_ss/a_ss*TICo_ss)/expo_ss;
h_ss        = (hX_ss^2 + hN_ss^2)^(1/2);
muh_ss      = h_ss^(VARPHI);
w_ss        = (wX_ss*hX_ss+wN_ss*hN_ss)/(hX_ss+hN_ss);
ppiwX_ss    = ppi_ss*a_ss;
ppiwN_ss    = ppi_ss*a_ss;
ppiw_ss     = ppi_ss*a_ss;


% EXON not defined in main file
e_pstr_ss    = 1;
e_pmstr_ss   = 1;
e_pCostr_ss  = 1;
del_F_ss     = ppistr_ss;
f_rel_ss     = 1;
pE_ss        = pT_ss;
pF_ss        = pT_ss;


auxi_ss     = 1;
for j=1:9
    eval(['auxi' num2str(j) '_ss    = 1;']);
end

% Indicate here the name of variables that are not logged
not_to_log={'bstr','tb','tby'};

% Indicate here the name of observables (all are assume dto have zero mean)
observables={'g_pib_obs','g_pib_x2_obs','g_pib_nt2_obs','g_pib_co_obs',...
    'g_cons_obs','g_inv_obs','gov_lt_obs','tby_obs', 'ppiN2_obs','ppiT2_obs',...
    'ppiM_obs','ppi_A_obs','ppi_E_obs', 'ppiWN_cmo_obs','ppiWX_cmo_obs',...
    'ppiS_obs','R_obs','RW_obs','count_prem_obs','ppi_star_obs',...
    'ppi_Co_star_obs','ppi_mstr_obs','pib_star_lt_obs',...
    'ppi_obs'};
    

%% Do note change this part, from here to the end
for iter = 1:length(M_.params)
    eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

if isfield(M_,'param_nbr') == 1
    if isfield(M_,'orig_endo_nbr') == 1
        NumberOfEndogenousVariables = M_.orig_endo_nbr;
    else
        NumberOfEndogenousVariables = M_.endo_nbr;
    end
    ys = zeros(NumberOfEndogenousVariables,1);
    for iter = 1:NumberOfEndogenousVariables
        varname = deblank(M_.endo_names(iter,:));
        aa=0;
        bb=0;
        for j=1: size(observables,2)
            bb=bb+ strcmp(observables{j},varname);
        end
        if bb==0
            for j=1: size(not_to_log,2)
                aa=aa+ strcmp(not_to_log{j},varname);
            end
            if aa==1
                eval(['ys(' int2str(iter) ') = ' varname '_ss;']);
            else
                eval(['ys(' int2str(iter) ') = log(' varname '_ss);']);
            end
        else
            eval(['ys(' int2str(iter) ') = 0;']);
        end
    end
else
    ys=zeros(length(lgy_),1);
    for iter = 1:length(lgy_)
        ys(iter) = eval(lgy_(iter,:));
    end
    check = 0;
end
if nargout>2
    param_update=M_.params;
end
end