%Steady state 
%Calibration to Mexico
%Time unit is a quarter

%%% utility params and taxes 
elast      = 8;           
EE         = 1 - 1/elast; %substitutability bw cf and ci (elasticity)
ETA        = 4;           %Frisch elasticity of labor supply
RBAR       = 0.027;       %SS country rate,EMBI(1998-2018)data=0.027
s          = 1;           %country spread
r          = RBAR;        %SS interest rate
BETA       = 1/(1+RBAR);  %HH's subjective discount factor
TAU        = 0.07;        %formal labor income tax rate
TAUNN      = 0.11;        %tax on formal firm's wage bill

%%% informal sector
urate      = 0.0384;  %SS unemployment rate
IEsize     = 0.15;     %size of informal economy
mismeasure = 0.6;     %what fraction of IE size is not included in national accounts
phiIM      = 1;       %fraction of ni (informal workers) substracted from unemployed(u=1-nf-phiIM*ni)
nni        = IEsize*(1-urate);
nnf        = 1-urate-nni;
empl       = nnf;
u          = 1 - nni - nnf;
ZETA       = 0.5;     %penalty when caught operating informal technology
NU         = 1;       %exponent on nni in piA function

%%% search and match params
SS         = 0.06;  %separation rate
ALPHAB     = 0.5;   %worker's bargaining power
ALPHAM     = 0.5;   %elasticity of (job) matching function
q          = 0.7;   %q=m/v :job finding rate/prob of fill vacancy
v          = SS*nnf/q;
theta      = v/(u+nni);  %vacancies to unemployment rate
OMEGA      = q*theta^(ALPHAM); %matching efficiency
BB         = 0;     %unemployment benefit
KAPPAv_yf  = 0.01;  

%%% production formal and informal params
ALPHAF     = 0.65;    %formal labor share
ALPHAI     = 0.7;     %informal labor share
d_yf       = 0.45;    %SS debt to GDP ratio Mexico data from FRED 1996-2019 45%
PSI        = 0.7;     %country interest rate debt elasticity

%%% shocks processes
RHOZF          = 0.78;  %persistence of TFP shock in formal sector
STD_EPS_ZF     = 1.168; %SD of TFP shock in percent in formal sector
RHOZI          = RHOZF;
STD_EPS_ZI     = STD_EPS_ZF;
RHOS           = 0.97;  %country spread persistence
STD_EPS_S      = 0.77;  %country spread volatility in percent

CORR_EPS_ZF_ZI = 0.0;
CORR_EPS_S_ZF = 0.0;
CORR_EPS_S_ZI = 0.0;

%%% remaining vars/parameters 
zf        = 1;    %tfp SS value
zi        = 1;    
tfpf      = zf;   %tfp SS value, transfer of state(z) to control var(tfp) for plotting IRFs/moments
tfpi      = zi;
yi        = zi*nni^ALPHAI; %informal output
yf        = zf*nnf^ALPHAF;   %formal output
KAPPA     = (KAPPAv_yf)/v;     %unit cost of posting a vacancy 
cf_yf     = 0.67;        %formal C/Y ratio, Mexico IMF data = 0.67 over 1990-2016
cf        = cf_yf*yf;    %formal consumption
DBAR      = yf*d_yf;     %SS level of debt
d         = DBAR;        


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% system of (implicit) equations solved using fsolve

% gov  = (TAU+TAUNN)*wf*nnf/(1-ZETA*(1-IOTA)/IOTA*(ci/cf)^(EE-1)*yi*nni^NU);
% ci   = (1-piA*ZETA)*yi;          
% IOTA = ( 1 + (1-TAU)/(1+TAUNN)*(cf/ci)^(EE-1)*nni/(ALPHAI*yi*(1-ALPHAB)*(1-piA*ZETA))*(ALPHAF*yf/nnf*(1-ALPHAB)+KAPPA/q*(1-SS-1/BETA-q*theta*ALPHAB)) )^(-1); 
% wf   = ALPHAB/(1+TAUNN)*(ALPHAF*yf/nnf + KAPPA*theta) + (1-ALPHAB)/(1-TAU)* ( (ETA*(nnf+nni)^(ETA-1))/(IOTA*cf^(EE-1)*(IOTA*cf^EE+(1-IOTA)*ci^EE)^(-1)) + BB ) ; 

F = @(X) [-X(1) + (TAU+TAUNN)*X(4)*nnf/(1-ZETA*(1-X(3))/X(3)*(X(2)/cf)^(EE-1)*yi*nni^NU); 
            - X(2) + (1-X(1)*nni^NU*ZETA)*yi;        
- X(3) + ( 1 + (1-TAU)/(1+TAUNN)*(cf/X(2))^(EE-1)*nni/(ALPHAI*yi*(1-ALPHAB)*(1-X(1)*nni^NU*ZETA))*(ALPHAF*yf/nnf*(1-ALPHAB)+KAPPA/q*(1-SS-1/BETA-q*theta*ALPHAB)) )^(-1); 
- X(4) + ALPHAB/(1+TAUNN)*(ALPHAF*yf/nnf + KAPPA*theta) + (1-ALPHAB)/(1-TAU)* ( (ETA*(nnf+nni)^(ETA-1))/(X(3)*cf^(EE-1)*(X(3)*cf^EE+(1-X(3))*X(2)^EE)^(-1)) + BB )];

X0 = [0.2; 0.4; 0.6; 1.8]; % initial values for gov(X1),ci(X2),IOTA(X3),wf(X4)
options = optimoptions('fsolve','Display','off');
[X,fval] = fsolve(F,X0,options);

gov   = X(1);
ci    = X(2);
IOTA  = X(3);
wf    = X(4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

piA = gov*nni^NU;
wi  = (1-IOTA)/IOTA*(ci/cf)^(EE-1) * ALPHAI*yi/nni * (1-piA*ZETA);
PHI = IOTA*cf^(EE-1)/( (IOTA*cf^EE + (1-IOTA)*ci^EE)*ETA*(nnf+nni)^(ETA-1) ) * (1-IOTA)/IOTA*(ci/cf)^(EE-1)*ALPHAI*yi/nni*(1-piA*ZETA);

tb        = yf-cf-KAPPA*v-gov + BB*(1-nnf); %trade balance in SS
tbyf      = tb/yf;
pc        = (1-IOTA)/IOTA*(ci/cf)^(EE-1);
SDF       = BETA;

%%% robustness on IE mismeasurement
omegaIM  = (1-mismeasure)*IEsize; %omegaIM=15% of IE (50% of IE=30% is captured) is included in Nat.accounts
thetaIM  = omegaIM*yf/( (1-IOTA)/IOTA*(ci/cf)^(EE-1) *(1-omegaIM)*yi);

ctotal   = cf + ci * (1-IOTA)/IOTA*(ci/cf)^(EE-1) * thetaIM;  
ytotal   = yf + yi * (1-IOTA)/IOTA*(ci/cf)^(EE-1) * thetaIM;
tbtotal  = ytotal - ctotal - KAPPA*v - gov + (1-nnf)*BB;
utotal   = 1 - nnf - phiIM*nni;
utotalf   = 1 - nnf - (1+1-phiIM)*nni;
tbtotalytotal = tbtotal/ytotal;
empltotal = nnf + nni;
ni_n = nni/empltotal;
irate = 1+r;