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_T = 1;
A_N = 1;
A_C = 1;

% Exogenous Prices:
P_C = 0.69;
R = 1.02;

% Debt stock in the steady state: 
D = 0.12;

% Subjective discount factor
%betaa =  1/R;

% Capital rate in the steady state:
mu_T = (1/betaa) - (1-delta);
mu_N = mu_T;
mu_C = mu_T;

% Total trade balance in the steady state:
TB = (R - 1)*D;

% Commdity sector: Computing Capital to Labor ratio: KL_C = K_C/L_C and W_C
KL_C = (alpha_C*P_C/mu_C)^(1/(1-alpha_C));
W_C = (1-alpha_C)*P_C*(KL_C^alpha_C);

% Tradable sector: Computing Capital to Labor ratio: KL_T = Kss_T/Lss_T and Wss_T
KL_T = ((alpha_T/mu_T)*(((gammaa_T*mu_T)/(P_C*alpha_T))^gammaa_T))^(1/(1 - alpha_T - gammaa_T));
W_T = (1 - gammaa_T - alpha_T)*((gammaa_T*mu_T/P_C*alpha_T)^gammaa_T)*KL_T^(alpha_T + gammaa_T); 

% Nontradable sector: {solve numerically for reer, Kss_N/Lss_N, Lss_N and Pss_N}
[FVAL y exitflag] = fsolve(@(t)(F(t, betaa, varphi, sigma, omega_T, omega_N, omega_C, chi, delta, alpha_T, alpha_N, alpha_C, gammaa_T, phi_T, phi_N, phi_C, v_C, v_D, rho_AC, rho_AT, rho_AN,  D, R, mu_T, mu_N, mu_C, TB, P_C, KL_C, W_C, KL_T, W_T)), [2;2;2;2], 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_N = FVAL(3);
P_N = FVAL(4);


% Using labor supply conditions:
L_C = (W_C/reer)^(1/(omega_C - 1));
L_T = (W_T/reer)^(1/(omega_T - 1));
W_N = (L_N^(omega_N - 1))*reer;

% Capital stock in the steady state is found by
K_C = KL_C*L_C;
K_T = KL_T*L_T;
K_N = KL_N*L_N;

% Using the definition of commodity demand:
CM_T = (gammaa_T/P_C)*(((gammaa_T*mu_T)/(P_C*alpha_T))^gammaa_T)*(KL_T^(alpha_T + gammaa_T))*L_T;

% Production in the steady state:
Y_T = (K_T^alpha_T)*(CM_T^gammaa_T)*(L_T^(1 - gammaa_T - alpha_T));
Y_C = (K_C^alpha_C)*(L_C^(1 - alpha_C));
Y_N = (K_N^alpha_N)*(L_N^(1 - alpha_N));

% Trade balance for Commodity and Tradable sector:
TB_C = P_C*(Y_C - CM_T);
TB_T = TB - TB_C;

% Investment in the steady state for each sector: 
I_T = delta*K_T;
I_C = delta*K_C;
I_N = delta*K_N;

% Consumption of Tradable and nontradable goods:
C_N = Y_N;
C_T = C_N*(((chi/(1 - chi))*P_N)^varphi);
C = (chi*(C_T^((varphi - 1)/varphi)) + (1-chi)*(C_N^((varphi - 1)/varphi)))^(varphi/(varphi-1));

% Output:
Y = Y_T +  P_C*Y_C + P_N*Y_N;

% Total Investment:
I = I_T + I_N + I_C;

% lagrange multiplier:
lambda = (C - (L_T/omega_T)^omega_T - (L_N/omega_N)^omega_N - (L_C/omega_C)^omega_C)*chi*((C/C_T)^(1/varphi));


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