var NF NI U S C CF PI LAMBDA R PAI ZF ZI Z_ZF Z_ZI WF THETA V MC G Z_G YF TLS Y SU CI;

varexo EPS_ZF EPS_ZI EPS_G;

parameters RHO ETA OMEGA TAU_C CHI VARPHI BETA BIG_M KAPPA EPSILON PHI PSI GAMMA RHO_ZF RHO_ZI b XI RHO_G TAU_W TLS_SS
           NF_SS NI_SS S_SS C_SS CF_SS PI_SS LAMBDA_SS R_SS PAI_SS ZF_SS ZI_SS Z_ZF_SS Z_ZI_SS WF_SS THETA_SS V_SS MC_SS  
           G_SS Z_G_SS U_SS YF_SS Y_SS SU_SS CI_SS; 

RHO=0.05;
ETA=0.75;
OMEGA=0.7;
TAU_C=0.21;
TAU_W=0.15;
CHI=0.789156338187430;
VARPHI=1;
BETA=0.99;
BIG_M=0.451276094896353;
KAPPA=0.024933241246562;
EPSILON=6;
PHI=60;
PSI=0.5;
GAMMA=0.836734923023498;
XI=0.5;
G_SS=0.166989086733278;

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;

%DECENTRALIZED STEADY STATE

PAI_SS=1.0026;
R_SS=PAI_SS/BETA; 
MC_SS=(EPSILON-1+PHI*(PAI_SS-1)*PAI_SS*(1-BETA))/EPSILON;
//ss0=[inival(3) inival(4) inival(2) inival(1) inival(6) inival(5) inival(11) inival(9) inival(7) inival(8)];
ss0=[0.7917 0.1183 0.623230439734604 0.452285754661247 0.777863134476648 1.180567500555204 0.458212756106031 0.0593775 0.831108041552078 0];
options = optimset('Display','iter', 'MaxFunEvals', 10000, 'MaxIter', 10000,'TolX',1e-20);
[steady,fval,exitflag]=fsolve(@ss_sol_mod,ss0,options,OMEGA, ETA, TAU_C, R_SS, VARPHI, CHI,BIG_M, XI, GAMMA, G_SS, PAI_SS, PHI, KAPPA, RHO, BETA, MC_SS, PSI, TAU_W);

f=ss;
NF_SS=steady(1);      
NI_SS=steady(2);
CF_SS=steady(3);
C_SS=steady(4);
PI_SS=steady(5);
LAMBDA_SS=steady(6);
THETA_SS=steady(7);
V_SS=steady(8);
WF_SS=steady(9);
TLS_SS=steady(10);

S_SS=1-(1-RHO)*NF_SS-NI_SS;
U_SS=1-NF_SS-NI_SS;
b=GAMMA*WF_SS;
MATCH_SS=BIG_M*(S_SS^XI)*(V_SS^(1-XI));
YF_SS=ZF_SS*NF_SS;
Y_SS=YF_SS+PI_SS*ZI_SS*NI_SS;
SU_SS=PI_SS*ZI_SS*NI_SS/Y_SS;
CI_SS=ZI_SS*NI_SS;

model;

1=NF+NI+U; //1
S=1-(1-RHO)*NF(-1)-NI; //2
C=(OMEGA*CF^ETA+(1-OMEGA)*(ZI*NI)^ETA)^(1/ETA); //3
OMEGA*(C^(-ETA))*(CF^(ETA-1))=(1+TAU_C)*LAMBDA; //4
(1-OMEGA)*(C^(-ETA))*((ZI*NI)^(ETA-1))=PI*LAMBDA*(1+(R-1)/R); //5
PI*ZI=CHI*(NI^VARPHI)/LAMBDA+BIG_M*(THETA^(1-XI))*(1-TAU_W)*WF+(1-BIG_M*(THETA^(1-XI)))*b; //6
1/R=BETA*(LAMBDA(+1)/LAMBDA)/PAI(+1); //7
THETA=V/S; //8
NF=(1-RHO)*NF(-1)+V*BIG_M*(THETA^-XI); //9
KAPPA/(BIG_M*(THETA^-XI))=MC*ZF-WF+BETA*(LAMBDA(+1)/LAMBDA)*(1-RHO)*KAPPA/(BIG_M*(THETA(+1)^-XI)); //10
EPSILON-1-EPSILON*MC+PHI*(PAI-1)*PAI-BETA*(LAMBDA(+1)/LAMBDA)*PHI*(PAI(+1)-1)*PAI(+1)*ZF(+1)*NF(+1)/(ZF*NF); //11
WF=PSI*MC*ZF+b*(1-PSI)/(1-TAU_W)+PSI*BETA*(LAMBDA(+1)/LAMBDA)*(1-RHO)*KAPPA*THETA(+1); //12
G+b*U=TAU_W*NF*WF+TAU_C*CF+PI*ZI*NI-PI(-1)*ZI(-1)*NI(-1)/PAI+TLS; //13
ZF*NF*(1-(PHI/2)*(PAI-1)^2)=CF+G+KAPPA*V; //14
YF=ZF*NF; //15


//EXOGENOUS PROCESSES
ZF=ZF_SS*exp(Z_ZF); //16
ZI=ZI_SS*exp(Z_ZI); //17
Z_ZF=RHO_ZF*Z_ZF(-1)+EPS_ZF; //18
Z_ZI=RHO_ZI*Z_ZI(-1)+EPS_ZI; //19
G=G_SS*exp(Z_G); //20
Z_G=RHO_G*Z_G(-1)+EPS_G; //21

//DEFINE Y, SU AND CI
Y=YF+PI*ZI*NI;
SU=PI*ZI*NI/Y;
CI=ZI*NI;
end;

initval; 
NF=NF_SS; 
NI=NI_SS; 
U=U_SS; 
S=S_SS;    
C=C_SS; 
CF=CF_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; 
WF=WF_SS; 
THETA=THETA_SS; 
V=V_SS; 
MC=MC_SS; 
G=G_SS; 
Z_G=Z_G_SS; 
TLS=TLS_SS;
YF=YF_SS;
Y=Y_SS;
SU=SU_SS;
CI=CI_SS;
end;

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

planner_objective(log(C)-CHI*(NI^(1+VARPHI))/(1+VARPHI));
ramsey_policy(planner_discount=0.99, irf=41, nodisplay, graph_format=pdf, order=1, instruments=(PAI)) Y NF NI U CF CI C PAI R WF V THETA SU;
//stoch_simul(irf=20, nodisplay, graph_format=pdf, order=1) NF NI U CF CI WF V PAI R;

Y_ZF=oo_.irfs.Y_EPS_ZF';
NF_ZF=oo_.irfs.NF_EPS_ZF';
NI_ZF=oo_.irfs.NI_EPS_ZF';
U_ZF=oo_.irfs.U_EPS_ZF';
CF_ZF=oo_.irfs.CF_EPS_ZF';
CI_ZF=oo_.irfs.CI_EPS_ZF';
C_ZF=oo_.irfs.C_EPS_ZF';
PAI_ZF=oo_.irfs.PAI_EPS_ZF';
R_ZF=oo_.irfs.R_EPS_ZF';
WF_ZF=oo_.irfs.WF_EPS_ZF';
V_ZF=oo_.irfs.V_EPS_ZF';
THETA_ZF=oo_.irfs.THETA_EPS_ZF';
SU_ZF=oo_.irfs.SU_EPS_ZF';

save('IRFs_ZF_RAMSEY', 'Y_ZF', 'NF_ZF', 'NI_ZF', 'U_ZF', 'CF_ZF', 'CI_ZF', 'C_ZF', 'PAI_ZF', 'R_ZF', 'WF_ZF', 'V_ZF', 'THETA_ZF', 'SU_ZF');

Y_G=oo_.irfs.Y_EPS_G';
NF_G=oo_.irfs.NF_EPS_G';
NI_G=oo_.irfs.NI_EPS_G';
U_G=oo_.irfs.U_EPS_G';
CF_G=oo_.irfs.CF_EPS_G';
CI_G=oo_.irfs.CI_EPS_G';
C_G=oo_.irfs.C_EPS_G';
PAI_G=oo_.irfs.PAI_EPS_G';
R_G=oo_.irfs.R_EPS_G';
WF_G=oo_.irfs.WF_EPS_G';
V_G=oo_.irfs.V_EPS_G';
THETA_G=oo_.irfs.THETA_EPS_G';
SU_G=oo_.irfs.SU_EPS_G';

save('IRFs_G_RAMSEY', 'Y_G', 'NF_G', 'NI_G', 'U_G', 'CF_G', 'CI_G', 'C_G', 'PAI_G', 'R_G', 'WF_G', 'V_G', 'THETA_G', 'SU_G');

for i=1:length(Y_ZF)
WELF_RAMSEY_ZF(i,1)=(BETA^(i-1))*(log(C_ZF(i))-CHI*(NI_ZF(i)^(1+VARPHI))/(1+VARPHI));
end

WELF_RAMSEY_ZF_1=sum(WELF_RAMSEY_ZF(1:4));
WELF_RAMSEY_ZF_5=sum(WELF_RAMSEY_ZF(1:20));
WELF_RAMSEY_ZF_10=sum(WELF_RAMSEY_ZF(1:40));

save('WELF_RAMSEY_ZF', 'WELF_RAMSEY_ZF_1', 'WELF_RAMSEY_ZF_5', 'WELF_RAMSEY_ZF_10');

for i=1:length(Y_G)
WELF_RAMSEY_G(i,1)=(BETA^(i-1))*(log(C_G(i))-CHI*(NI_G(i)^(1+VARPHI))/(1+VARPHI));
end

WELF_RAMSEY_G_1=sum(WELF_RAMSEY_G(1:4));
WELF_RAMSEY_G_5=sum(WELF_RAMSEY_G(1:20));
WELF_RAMSEY_G_10=sum(WELF_RAMSEY_G(1:40));

save('WELF_RAMSEY_G', 'WELF_RAMSEY_G_1', 'WELF_RAMSEY_G_5', 'WELF_RAMSEY_G_10');