% Function to calculate the steady state in DYNARE.
% This replaces the command steady since such steady state must be
% calculated bearing in mind that there is an endogenous parameter B0 which
% needs to be calibrated simultaneously with the steady state.

function [ys,check] = Underground_economy_v1_steadystate(ys,exe)

% % -------------------------------------- % %
% % ------ Extraccion de parametros ------ % %
% % -------------------------------------- % %

global M_

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

% % --------------------------------------------- % %
% % ------ C?lculo del estado estacionario ------ % %
% % --------------------------------------------- % %
%% EXOGENOUS PROCESSES
tauc = tauc_bar;
tauh = tauh_bar;
taus = taus_bar;
p    = p_bar;
a    = a_bar;
b    = b_bar;
xii  = xii_bar; 
xih  = xih_bar;
%% ENDOGENOUS VARS
% % Impositions:
% Labor:
h    = log(h_bar);                                         %Calibrated from ISTAT.See Orsi                                                     %et al.(2014), p.424.
hu   = log(exp(h)*hu_h_bar);                             %0.13 calibrated from ISTAT in paper.
hm   = log(exp(h)-exp(hu));

% Euler:
%r    = [1-bet*(1-deltak)]/[bet*(1-tauh)];
% r    = log((exp(gamma_Q/100)-bet*(1-deltak))/(bet*(1-exp(tauh))));    %Taken from sis_matrix(143).
% r    = log((1-bet*(1-deltak))/(bet*(1-exp(tauh))));    %Taken from sis_matrix(143).
r = log((1/bet -  (1-deltak))/(1-exp(tauh))) ;
ru   = log((1-exp(tauh))*exp(r));                                   %Eq(5) in notes, Eq(14) in paper. see Orsi et al. (2014) p.424.

% Firm Demands for PF:
khm  = log(((1-exp(tauc))*(1-alp)/exp(r))^(1/alp));             % AUX for SS         %same as sis_matrix(70)
khu  = log(((1-exp(tauc)*exp(p)*s)*(1-alpu)/exp(ru))^(1/alpu)); % AUX for SS         %same as sis_matrix(71)
wm   = log(alp*(1-exp(tauc))/(1-exp(tauc)+exp(taus))*exp(khm)^(1-alp));       %same as sis_matrix(162)
wu   = log(alpu*exp(khu)^(1-alpu));                            %same as sis_matrix(164)
km   = log(exp(khm)*exp(hm));                                       %same as sis_matrix(185) 
ku   = log(exp(khu)*exp(hu));                                       %same as sis_matrix(186)
k    = log(exp(ku)+exp(km));                                        %same as sis_matrix(187)


% Capital LOM and Resource constraint:
% iv   = log((exp(gamma_Q/100)-1+deltak)*exp(k));                         %comes from sis_matrix(189)
iv   = log(deltak*exp(k)); 
c    = log((1-exp(tauh))*(exp(wm)*exp(hm)+exp(r)*exp(km))+ exp(wu)*exp(hu)+exp(ru)*exp(ku)-exp(iv));         %sis_matrix(191-192)
                                                     %the difference is "iv".
% Production functions:
ym   = log(exp(km)^(1-alp)*exp(hm)^alp);                            %same as sis_matrix(194)
yu   = log(exp(ku)^(1-alpu)*exp(hu)^alpu);                          %same as sis_matrix(195)
y    = log(exp(ym)+exp(yu));                                        %missing from sis_matrix
% Labor Supplies:
B0   = ( (1-exp(tauh))*exp(wm)*exp(c)^-sig )*exp(h)^-xi ;                   %same as sis_matrix(203).                           
B1   = (exp(wu)*exp(c)^-sig - B0*exp(h)^xi)*exp(hu)^-phi;          %same as sis_matrix(84). 

% Government Variables:
gc   = log(exp(tauc)*(exp(ym)-exp(hm)*exp(wm)+exp(p)*s*(exp(yu)-exp(wu)*exp(hu))));               %same as sis_matrix(198-199).
gh   = log(exp(tauh)*(exp(wm)*exp(hm)+exp(r)*exp(km)));                            %same as sis_matrix(202).
gs   = log(exp(taus)*exp(wm)*exp(hm));                                   %same as sis_matrix(200).
frev = log(exp(gc)+exp(gh)+exp(gs));
ev   = log(exp(wu)*exp(hu)*exp(taus)+(1-exp(p))*exp(tauc)*(exp(yu)-exp(wu)*exp(hu))+...
        exp(tauh)*(exp(wu)*exp(hu)+exp(ru)*exp(ku)));                         %missing from sis_matrix

rat  = log(exp(yu)/exp(y));                                         %missing from sis_matrix.

%----------------------
% Measurement Equations
%----------------------
    % Consumption growth
    C_obs = gamma_Q;
    
    % Investment growth
    X_obs = gamma_Q;
    
    % Corporate taxation earnings growth
    G_c_obs = gamma_Q;
    
    % Social security taxation earnings growth
    G_s_obs = gamma_Q;
    
    % Personal income taxation earnings growth
    G_h_obs = gamma_Q;
    
    % Gross earnings in the regular market growth
    W_h_obs = gamma_Q;
    
    % Average probability of being inspected
    l_p_obs = log(p_mean) + p;
    p_obs   = p_mean + p;

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