var NF NI U C CF CI M PI LAMBDA R PAI ZF ZI Z_ZF Z_ZI F WF THETA V Q MATCH MC YF YI G Z_G TLS; 

varexo EPS_ZF EPS_ZI EPS_G EPS_R;

parameters RHO ETA OMEGA TAU_C CHI VARPHI BETA BIG_M KAPPA EPSILON PHI PSI GAMMA ALPHA_R ALPHA_Y ALPHA_PAI RHO_ZF RHO_ZI RHO_G b XI ALPHA_U
           NF_SS NI_SS C_SS CF_SS CI_SS M_SS PI_SS LAMBDA_SS R_SS PAI_SS ZF_SS ZI_SS Z_ZF_SS Z_ZI_SS F_SS WF_SS THETA_SS V_SS Q_SS MATCH_SS MC_SS YF_SS 
           YI_SS G_SS Z_G_SS U_SS NI_TO_N RHO_R TAU_W TLS_SS;

ALPHA_R=0.8;
ALPHA_PAI=1.5;
ALPHA_Y=0.5/4;
ALPHA_U=0.5/4;
RHO_ZF=0.9;
RHO_ZI=0.9;
RHO_G=0.7;
ZF_SS=1;
ZI_SS=1;
Z_ZF_SS=0;
Z_ZI_SS=0;
Z_G_SS=0;
RHO_R=0.7;

%DECENTRALIZED STEADY STATE

BETA=0.99;
PAI_SS=1; %1
EPSILON=6;
PHI=60;
ZF_SS=1;
MC_SS=(EPSILON-1+PHI*(PAI_SS-1)*PAI_SS*(1-BETA))/EPSILON;
VBAR=0.045; %2
RHO=0.025;
Q_SS=2/3; %3
U_SS=0.09; %4
NI_TO_N=0.13; %5
NI_SS=NI_TO_N*(1-U_SS);
NF_SS=1-NI_SS-U_SS;
XI=0.5;
V_SS=RHO*NF_SS/Q_SS;
THETA_SS=V_SS/U_SS;
BIG_M=Q_SS*THETA_SS^XI;
KAPPA=BETA*MC_SS*ZF_SS*Q_SS/(1-BETA*(1-RHO)+BETA/VBAR);
WF_SS=(KAPPA/Q_SS)/VBAR;
PSI=0.5;
TAU_W=0.15;
GAMMA=(WF_SS-PSI*MC_SS*ZF_SS-PSI*KAPPA*THETA_SS)/(WF_SS*(1-PSI)/(1-TAU_W));
F_SS=BIG_M*THETA_SS^(1-XI);
OMEGA=0.7;
ETA=0.75;
TAU_C=0.21;
VARPHI=1;
YI_SS=NI_SS;
CI_SS=YI_SS;
R_SS=PAI_SS/BETA;
YF_SS=ZF_SS*NF_SS;
ZI_SS=1;

ss0=[0.65 0.7 2 0.15 0.6 0.1 1.5];

options = optimset('Display','iter', 'MaxFunEvals', 10000, 'MaxIter', 10000,'TolX',1e-20);
[ss,fval,exitflag]=fsolve(@ss_sol,ss0,options,BETA, OMEGA, ETA, NI_SS, NF_SS, TAU_C, TAU_W, R_SS, ZI_SS, VARPHI, F_SS, GAMMA, U_SS, PAI_SS, YF_SS, PHI, V_SS, CI_SS, WF_SS, KAPPA);

f=ss;
C_SS=ss(1);
CF_SS=ss(2); 
CHI=ss(3);
G_SS=ss(4);
PI_SS=ss(5);
M_SS=ss(6);
LAMBDA_SS=ss(7);

b=GAMMA*WF_SS;
MATCH_SS=BIG_M*(U_SS^XI)*(V_SS^(1-XI));
TLS_SS=0;

model; 
1=NF+NI+U; //1
C=(OMEGA*CF^ETA+(1-OMEGA)*CI^ETA)^(1/ETA); //2
M=PI*CI; //3
OMEGA*(C^(-ETA))*(CF^(ETA-1))=(1+TAU_C)*LAMBDA; //4
(1-OMEGA)*(C^(-ETA))*(CI^(ETA-1))=PI*LAMBDA*(1+(R-1)/R); //5
PI*ZI=CHI*(NI^VARPHI)/LAMBDA+BETA*(LAMBDA(+1)/LAMBDA)*(F*(1-TAU_W)*WF(+1)+(1-F)*b); //6
1/R=BETA*(LAMBDA(+1)/LAMBDA)/PAI(+1); //7
THETA=V/U; //8
MATCH=BIG_M*(U^XI)*(V^(1-XI)); //9
Q=MATCH/V; //10
F=MATCH/U; //11
NF=(1-RHO)*NF(-1)+MATCH(-1); //12
KAPPA/Q=BETA*(LAMBDA(+1)/LAMBDA)*(MC(+1)*ZF(+1)-WF(+1)+(1-RHO)*KAPPA/Q(+1)); //13
EPSILON-1-EPSILON*MC+PHI*(PAI-1)*PAI-BETA*(LAMBDA(+1)/LAMBDA)*PHI*(PAI(+1)-1)*PAI(+1)*YF(+1)/YF=0; //14
WF=PSI*MC*ZF+(1-PSI)*b/(1-TAU_W)+PSI*KAPPA*THETA; //15
G+b*U=TAU_W*NF*WF+TAU_C*CF+M-M(-1)/PAI+TLS; //16
YF=ZF*NF; //17
YI=ZI*NI; //18
YI=CI; //19
YF*(1-(PHI/2)*(PAI-1)^2)=CF+G+KAPPA*V; //20
log(R/R_SS)=ALPHA_R*log(R(-1)/R_SS)+(1-ALPHA_R)*(ALPHA_PAI*log(PAI/PAI_SS)+ALPHA_Y*log(YF/YF(-1))-ALPHA_U*log(U/U_SS))-EPS_R; //21

//EXOGENOUS PROCESSES
ZF=ZF_SS*exp(Z_ZF); //22
ZI=ZI_SS*exp(Z_ZI); //23
Z_ZF=RHO_ZF*Z_ZF(-1)+EPS_ZF; //24
Z_ZI=RHO_ZI*Z_ZI(-1)+EPS_ZI; //25
G=G_SS*exp(Z_G); //26
Z_G=RHO_G*Z_G(-1)+EPS_G; //27
end;

steady_state_model; 
NF=NF_SS; 
NI=NI_SS; 
U=U_SS; 
C=C_SS; 
CF=CF_SS; 
CI=CI_SS; 
M=M_SS; 
PI=PI_SS; 
LAMBDA=LAMBDA_SS; 
R=R_SS; 
PAI=PAI_SS; 
ZF=ZF_SS; 
ZI=ZI_SS; 
Z_ZF=Z_ZF_SS; 
Z_ZI=Z_ZI_SS; 
F=F_SS; 
WF=WF_SS; 
THETA=THETA_SS; 
V=V_SS; 
Q=Q_SS; 
MATCH=MATCH_SS; 
MC=MC_SS; 
YF=YF_SS; 
YI=YI_SS; 
G=G_SS; 
Z_G=Z_G_SS; 
TLS=TLS_SS;
end;

resid; 
steady; 
check;

shocks; 
var EPS_ZF; 
stderr 0.01; 
var EPS_ZI; 
stderr 0.01;
var EPS_G; 
stderr 0.01;
var EPS_R; 
stderr 0.01;
end;

stoch_simul(irf=20, nodisplay, graph_format=pdf, order=1, pruning);