function [ys,check] = master_steadystate(ys,exo)
% function [ys,check] = master_steadystate(ys,exo)
% computes the steady state for the master.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 for 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;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Enter model equations here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

options=optimset('Display','Final','TolX',1e-10,'TolFun',1e-10, 'MaxIter', 10000); % set options for numerical solver

% Productivity factor for each sector:
A_Tss = 1;
A_Nss = 1;
A_Css = 1;

% Exogenous Prices:
P_Css = Pss_C;

% Interest rate in the steady state:
r =  RPss + r_wss;

% Capital rate in the steady state:
mu_Tss = (1/bbeta) - (1 - ddelta);
mu_Nss = mu_Tss;
mu_Css = mu_Tss;

% Total trade balance in the steady state:
TBss = r*Dss;

% Commdity sector: Computing Capital to Labor ratio: KL_C = K_C/L_C and W_C:
KL_C = ((aalpha_C*P_Css*A_Css)/(mu_Css*(1 + eeta_C*(r/(1+r)))))^(1/(1 - aalpha_C));
W_Css = (((1 - aalpha_C)*P_Css*A_Css)/(1 + eeta_C*(r/(1+r))))*(KL_C^aalpha_C);


% Tradable sector: Computing Capital to Labor ratio: KL_T = Kss_T/Lss_T and Wss_T
theta = ((ggamma_T*mu_Tss)/(P_Css*aalpha_T))^ggamma_T;
KL_T = ((theta*aalpha_T*A_Tss)/(mu_Tss*(1 + eeta_T*(r/(1 + r)))))^(1/(1 - aalpha_T - ggamma_T));
W_Tss = ((1 - aalpha_T - ggamma_T)/(1 + eeta_T*(r/(1+r))))*(theta*A_Tss)*(KL_T^(aalpha_T + ggamma_T));

	
% Nontradable sector: {solve numerically for reer, Kss_N/Lss_N, Lss_N and Pss_N}
[FVAL y exitflag] = fsolve(@(t)(F(t, b, bbeta, vvarphi, ssigma, oomega_T, oomega_N, oomega_C, cchi, ddelta, aalpha_T, aalpha_N, aalpha_C, ggamma_T, pphi_T, pphi_N, pphi_C, v_C, v_D, rho_AC, rho_AT, rho_AN,  Dss, r, mu_Tss, mu_Nss, mu_Css, TBss, P_Css, KL_C, W_Css, KL_T, W_Tss, eeta_N, eeta_T, eeta_C, A_Nss, A_Tss, A_Css, theta)), [.5;2;1;1], options);

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

reer = FVAL(1);
KL_N = FVAL(2);
L_Nss = FVAL(3);
P_Nss = FVAL(4);



% Using labor supply conditions:
L_Css = (((1 - b*bbeta)*W_Css)/reer)^(1/(oomega_C - 1));
L_Tss = (((1 - b*bbeta)*W_Tss)/reer)^(1/(oomega_T - 1));
W_Nss = (reer*(L_Nss^(oomega_N - 1)))/(1 - b*bbeta);

% Capital stock in the steady state is found by
K_Css = KL_C*L_Css;
K_Tss = KL_T*L_Tss;
K_Nss = KL_N*L_Nss;

% Using the definition of optimal commodity demand:
CM_Tss = ((ggamma_T*mu_Tss)/(P_Css*aalpha_T))*K_Tss;

% Production in the steady state:
Y_Tss = A_Tss*(K_Tss^aalpha_T)*(CM_Tss^ggamma_T)*(L_Tss^(1 - ggamma_T - aalpha_T));
Y_Css = A_Css*(K_Css^aalpha_C)*(L_Css^(1 - aalpha_C));
Y_Nss = A_Nss*(K_Nss^aalpha_N)*(L_Nss^(1 - aalpha_N));

% Trade balance for Commodity and Tradable sector:
TB_Css = P_Css*(Y_Css - CM_Tss);
TB_Tss = TBss - TB_Css;

% Investment in the steady state for each sector: 
I_Tss = ddelta*K_Tss;
I_Css = ddelta*K_Css;
I_Nss = ddelta*K_Nss;

% Consumption of Tradable and nontradable goods:
C_Nss = Y_Nss;
C_Tss = C_Nss*(((cchi/(1 - cchi))*P_Nss)^vvarphi);
Css = (cchi*(C_Tss^((vvarphi - 1)/vvarphi)) + (1 - cchi)*(C_Nss^((vvarphi - 1)/vvarphi)))^(vvarphi/(vvarphi-1));


% Debt for NonTradable Sector:
D_Nss = (1/(1+r))*eeta_N*(W_Nss*L_Nss + mu_Nss*K_Nss); 

% Debt for Commodity Sector:
D_Css = (1/(1+r))*eeta_C*(W_Css*L_Css + mu_Css*K_Css);

% Debt for Tradable Sector:
D_Tss = (1/(1+r))*eeta_T*(W_Tss*L_Tss + mu_Tss*K_Tss + Pss_C*CM_Tss);

% Debt for Households: 
D_Hss = Dss - D_Tss - D_Nss - D_Css; 
%D_Hss = (W_Tss*L_Tss + W_Nss*L_Nss + W_Css*L_Css + mu_Tss*K_Tss + mu_Nss*K_Nss + mu_Css*K_Css - C_Tss - P_Nss*C_Nss - I_Tss - I_Nss - I_Css)/r; 


% Total Output in the numeraire tradable units: 
Yss = Y_Tss + P_Nss*Y_Nss + TB_Css;

% Total Investment in tradable unit:
Iss = I_Tss + I_Nss + I_Css;




% 
UtCtss = ((1 - b)*Css - (L_Tss^oomega_T)/oomega_T - (L_Nss^oomega_N)/oomega_N - (L_Css^oomega_C)/oomega_C)^(-ssigma);




% 
PHI_Tss = 0;
PHI_Nss = 0;
PHI_Css = 0;

%
PHI_Ttss = 0;
PHI_Ntss = 0;
PHI_Ctss = 0;



% 
A_CTss = cchi*((Css/C_Tss)^(1/vvarphi));

A_CNss = (1 - cchi)*((Css/C_Nss)^(1/vvarphi));



% lagrange multiplier:
lambdass = (1 - b*bbeta)*A_CTss*UtCtss;







C_T = log(C_Tss);
C_N = log(C_Nss);
I_T = log(I_Tss);  	
I_N = log(I_Nss);
I_C = log(I_Css);
L_T = log(L_Tss);
L_N = log(L_Nss);
L_C = log(L_Css);
D_H = log(D_Hss);
lambda = log(lambdass);

A_N = log(A_Nss);
K_N = log(K_Nss);
Y_N = log(Y_Nss);
%M_N = log(M_Nss);
D_N = log(D_Nss);

A_C = log(A_Css);
K_C = log(K_Css);
Y_C = log(Y_Css);
%M_C = log(M_Css);
D_C = log(D_Css);
TB_C = log(TB_Css);

A_T = log(A_Tss);
K_T = log(K_Tss);
CM_T = log(CM_Tss);
Y_T = log(Y_Tss);
%M_T = log(M_Tss);
D_T = log(D_Tss);
TB_T = log(TB_Tss);

P_N = log(P_Nss);
P_C = log(P_Css);
mu_T = log(mu_Tss);
mu_N = log(mu_Nss);
mu_C = log(mu_Css);
W_T = log(W_Tss);
W_N = log(W_Nss);
W_C = log(W_Css);
r = r;
RP = RPss;
r_w = r_wss;


Y = log(Yss);
C = log(Css);
I = log(Iss);
TB = log(TBss);
D = log(Dss);

UtCt = log(UtCtss);

A_CT = log(A_CTss);
A_CN = log(A_CNss);

PHI_T = log(PHI_Tss);
PHI_N = log(PHI_Nss);
PHI_C = log(PHI_Css);
PHI_Tt = log(PHI_Ttss);
PHI_Nt = log(PHI_Ntss);
PHI_Ct = log(PHI_Ctss);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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