close all;

%format long

%%%%%%%%%%%%%%%%%%%Defining variables%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 var c iB piS pS mP d iD iV hd pH piH n omega y yD yF p pD pF yX yS e wpF wpX pX k iR rK
     piD mc I iL rhot sigmat iFB iW LFB LI q RFT BFP LW vE vR v LCB NFA thetaFB;


 varexo eps_W ; 


 parameters
beta THETA_K delta eta_D mu_R kappa phi_1 phi_2 phi_3 eta_F rho eta_x nu sigma THETA_V 
eta_H eta_N eta LAMBDA_D alpha theta_D phi_I gamma_vv gamma_v psi phi_F mu_F chi eps_1 eps_2 
phiR1 phiR2 phi_q theta_FB theta_FP GAMMA_Y veta phi_E rho_W theta_C

yG_ss iR_ss iB_ss iL_ss iD_ss iV_ss iW_ss r_ss piS_ss K_ss I_ss rK_ss q_ss PH_ss LI_ss V_ss 
L_ss mP_ss C_ss d_ss Hd_ss PS_ss piH_ss BFP_ss N_ss omega_ss Y_ss YD_ss YF_ss P_ss PD_ss
PF_ss YX_ss E_ss PX_ss LW_ss piDG_ss mc_ss piD_ss VE_ss VR_ss LCB_ss LFB_ss RF_ss
RFT_ss T_ss thetaFP_ss thetaFB_ss NFA_ss
Hbar WPF WPX Bbar BCbar BPbar
YS_ss;

%%%%%%%%Calibration%%%%%%%
sigma=0.6;
beta=0.93;
eta_x=1.5;
nu=0.2;
mu_R=0.1;
eta_N=1.5;
LAMBDA_D=0.7;
eta=0.6;
mu_F=0.5;
GAMMA_Y=0.5 ;
veta=0.6;
alpha=0.35;
theta_D=10;
phi_I=74.5;
delta=0.01;
THETA_K=14;
rho=0.08;
gamma_v=0.08;
theta_FB = 1; 
theta_FP = 1;
phi_1=0.03;
phi_2=0;
phi_3=1.5;
chi=0;
eps_1=2.5; 
eps_2=0.5;
phi_F=1;
phi_E=0.5;
phi_q=0.05;
theta_C=0.5; 
psi=0.5; 

eta_H=0.02;
THETA_V=0.5;
kappa=0.2;
gamma_vv=0.004;
eta_D=0;
eta_F=0; 
phiR1=1; 
phiR2=1;
rho_W=0.8;
theta_C=0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Steady state%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iB_ss = beta^(-1) - 1; %(4)
r_ss = beta^(-1) - 1; %(4)
iR_ss = beta^(-1) - 1; %(4)

%(THETA_K/2)*((K_ss/K_ss)-1)^2*K_ss = 0; %(5)
I_ss = delta*K_ss; %(6)
rK_ss = (1+iL_ss)*(1-beta*(1-delta)); %(7)

iD_ss = (1+(1/eta_D))^(-1)*(1-mu_R)*iR_ss; %(8)
q_ss = ((kappa*PH_ss*Hbar)/LI_ss)^phi_1*(V_ss*L_ss)^phi_2; %(9)

iL_ss = (1/(1+eta_F^(-1)*q_ss))*((1-rho)*beta^(-1)+rho*((1+iV_ss)+gamma_v))-1; %(10)

%mP_ss = (eta_x*nu*C_ss^(1/sigma)*(1+iB_ss))/iB_ss; %(11)
%d_ss = (eta_x*(1-nu)*C_ss^(1/sigma)*(1+iB_ss))/(iB_ss-iD_ss); %(12)   
%(V_ss/PS_ss) = (1/THETA_V)*((iV_ss-iB_ss)/(1+iB_ss)); %(13)
%(PH_ss*Hd_ss)/PS_ss = (1-((1+piH_ss)/(1+iB_ss)))^(-1)*(eta_H/C_ss^(-1/sigma)); %(14)
%BFP_ss = ((1+iW_ss)-(1+iB_ss))/(theta_FP*(1+iW_ss)); %(15)  
                                         
mP_ss = (eta_x*nu*C_ss^(1/sigma))/(1-beta); %(16)
d_ss = (eta_x*(1-nu)*C_ss^(1/sigma))/((1-beta)*mu_R); %(17)   
V_ss = ((beta/THETA_V)*(1+iV_ss-beta^(-1)))*PS_ss; %(18)
Hd_ss = (((1-beta)^(-1))*(eta_H/C_ss^(-1/sigma)))*(PS_ss/PH_ss); %(19)
BFP_ss = (iW_ss-iB_ss)/(theta_FP*(1+iW_ss)); %(20)

iV_ss = ((THETA_V*V_ss)/(beta*PS_ss))+beta^(-1)-1; %(21)  

N_ss = 1-((eta_N*C_ss^(1/sigma))/omega_ss); %(22)

Y_ss = (LAMBDA_D*YD_ss^((eta-1)/eta)+(1-LAMBDA_D)*YF_ss^((eta-1)/eta))^(eta/(eta-1)); %(23)
YD_ss = LAMBDA_D^eta*(PD_ss/P_ss)^(-eta)*Y_ss; %(24)
YF_ss = (1-LAMBDA_D)^eta*(PF_ss/P_ss)^(-eta)*Y_ss; %(24)
P_ss = (LAMBDA_D^eta*(PD_ss)^(1-eta)+(1-LAMBDA_D)^eta*(PF_ss)^(1-eta))^(1/(1-eta)); %(25)
PF_ss = E_ss*WPF; %(26)
Y_ss = (GAMMA_Y*YS_ss^(1+(1/veta))+(1-GAMMA_Y)*YX_ss^(1+(1/veta)))^(veta+(1+veta)); %(27)
PX_ss = E_ss*WPX; %(28)
Y_ss = (((1-GAMMA_Y)/GAMMA_Y)^veta*(PS_ss/PX_ss)^veta)*YX_ss; %(29)

YD_ss = N_ss^(1-alpha)*K_ss^alpha; %(30)
%K_ss = ((alpha/(1-alpha))*(((1+iR_ss)*omega_ss)/rK_ss))*N_ss; %(32)
omega_ss = ((1-alpha)/alpha)*((K_ss*(beta^(-1)-1+delta))/(N_ss*(1+iR_ss))); %(32)

LW_ss = PS_ss*omega_ss*N_ss; %(33)

%piDG_ss = 1;
mc_ss=((1-theta_D)+(phi_I*((piD_ss/piDG_ss)-1)*(piD_ss/piDG_ss))-beta*phi_I*((piD_ss/piDG_ss)-1)*(piD_ss/piDG_ss))/theta_D;
mc_ss = (theta_D-1)/theta_D; %(34)

LI_ss = PS_ss*I_ss; %(35)

VE_ss = ((gamma_vv/(iV_ss+gamma_v-beta^(-1)+1))^(1/(1-phi_E)))*PS_ss; %(36)
VR_ss = rho*LI_ss; %(37)
V_ss = VR_ss+VE_ss; %(38)

LCB_ss = LW_ss+LI_ss-E_ss*LFB_ss-(1-mu_R)*d_ss*PS_ss-V_ss; %(39)
LFB_ss = (beta^(-1)-(1+iW_ss))/(theta_FB*(1+iW_ss)); %(40)

RF_ss = (phiR1*WPF*YF_ss)^phi_F*(phiR2*(LFB_ss-BFP_ss))^(1-phi_F);
RFT_ss = (phiR1*WPF*YF_ss)^phi_F*(phiR2*(LFB_ss-BFP_ss))^(1-phi_F);

YS_ss = (C_ss+delta*K_ss)/(1-psi); %(41)
PS_ss = (P_ss*Y_ss-PX_ss*YX_ss)/YS_ss; %(42)
%BCbar/PS_ss = (V_ss/PS_ss)+eta_x*C_ss^(1/sigma)*(1+iB_ss)*((nu/iB_ss)+((1-nu)/(iB_ss-iD_ss))); 
BCbar = ((V_ss/PS_ss)+((eta_x*C_ss^(1/sigma))/(1-beta))*(nu+((1-nu)/mu_R)))*PS_ss; %(43)
BPbar = Bbar - BCbar; %(44)
T_ss = psi*PS_ss*Y_ss+iB_ss*BPbar-iR_ss*LCB_ss-iW_ss*E_ss*RF_ss; %(45)

LFB_ss=((WPX*YX_ss)-WPF*YF_ss+iW_ss*NFA_ss+thetaFP_ss*BFP_ss)/thetaFB_ss; %(46)
NFA_ss = RF_ss+BFP_ss-LFB_ss; %(47)





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model(linear);
c(+1) = c+sigma*(iB-piS(+1)); %(1)
piS(+1) = pS(+1)-pS; %(2)

%mP*mP_ss = ((eta_x*nu*C_ss^(1/sigma))/(1-beta))*((c/sigma)-(beta/(1-beta))*iB); %(3)-----
mP = (c/sigma)-(beta/(1-beta))*iB; %(4)
d = (c/sigma)+((1-mu_R+mu_R*beta)/(mu_R*(1-beta)))*(iD-iB); %(5)
v-pS = iV_ss/(iV_ss-iB_ss)*(iV-iB); %(6)
hd = pS-pH+((piH_ss*(piH(+1)-iB))/(iB_ss-piH_ss))+(c/sigma); %(7)

%n*N_ss = ((eta_N*C_ss^(1/sigma))/(omega_ss))*omega - ((eta_N*C_ss^(1/sigma)*c)/(sigma*omega_ss));% ()
n = ((eta_N*C_ss^(1/sigma))/(omega_ss-eta_N*C_ss^(1/sigma)))*(omega-(c/sigma)); %(8)
y = (LAMBDA_D*(YD_ss/Y_ss)^((eta-1)/eta))*yD+((1-LAMBDA_D)*(YF_ss/Y_ss)^((eta-1)/eta))*yF; %(9)

yD = eta*(p-pD)+y; %(9bis)
yF = eta*(p-pF)+y; %(9ter)

p = (LAMBDA_D^eta*(PD_ss/P_ss)^(1-eta))*pD+((1-LAMBDA_D)^eta*(PF_ss/P_ss)^(1-eta))*pF; %(10)
pF = mu_F*(e+wpF)+(1-mu_F)*pF(-1); %(11)
y = GAMMA_Y*(YS_ss/Y_ss)^((1+veta)/veta)*yS+(1-GAMMA_Y)*(YX_ss/Y_ss)^((1+veta)/veta)*yX; %(12)
yS = yX+veta*(pS-pX); %(13)
pX = e+wpX; %(14)
yD = (1-alpha)*n+alpha*k; %(15)
n = k-iR-omega+((1+rK_ss)/rK_ss)*rK; %(16)
piD = ((theta_D-1)/phi_I)*mc+beta*piD(+1); %(17)
mc = (1-alpha)*(iR+omega)+alpha*(rK_ss/(rK_ss-1))*rK; %(18)
I = (1/delta)*(k(+1)-(1-delta)*k); %(19)
rK(+1) = ((1+iL_ss)/(1+rK_ss))*((iL-piS(+1))+THETA_K*(k(+1)-k))-((beta*(1+iL_ss))/(1+rK_ss))*((1-delta)*(iL(+1)-iB(+1))+THETA_K*(k(+2)-k(+1))); %(20)
iD = ((1-mu_R)/(1-(1-beta)*mu_R))*iR; %(21)
iL = (1/((1+iL_ss)*q_ss))*(((1-rho)/beta)*iR+rho*(1+iV_ss)*iV-(1+rho)*(1+iV_ss+gamma_v-(1/beta))*rhot+rho*((1+iV_ss)+gamma_v-(1/beta))*sigmat)-q;%(22)

iFB = iW+e(+1)-e+(((1+iW_ss)*(theta_FB/2)*LFB_ss)/((1+iW_ss)*(1+(theta_FB/2)*LFB_ss)))*LFB;
q = phi_1*(pH-LI)+phi_2*sigmat+phi_3*y; %(23)
iR = chi*iR(-1)+(1-chi)*(eps_1*piS+eps_2*(y-Y_ss)); %(24) -----
RFT = phi_F*(wpF+yF)+(1-phi_F)*((LFB_ss*LFB-BFP_ss*BFP)/(LFB_ss-BFP_ss)); %(25)
LW = n+omega+pS; %(26)
LI = pS+I; %(27)

%LW = (1/LFB_ss)*LW_ss*LW+LI_ss*LI;%-----
LW_ss*LW+LI_ss*LI = 0;
vE-pS = (1/(1-phi_E))*((beta^(-1)*iR-(1+iV_ss)*iV)/(1+iV_ss+gamma_v-beta^(-1))); %(28)
vR = ((1+rho)/rho)*rhot+sigmat+LI; %(29)
sigmat = -phi_q*q; %(30)

v = (VR_ss/V_ss)*vR+(VE_ss/V_ss)*vE; 

LCB = (1/LCB_ss)*(LW_ss*LW+LI_ss*LI-E_ss*LFB_ss*(e+LFB)-(1-mu_R)*PS_ss*d_ss*(d+pS)-V_ss*v); %(31)
LFB = ((1+iR_ss)/(iR_ss-iW_ss))*(iR-iW-e(+1)+e); %(32)

%(mP_ss*PS_ss)*(mP+pS)+(d_ss*PS_ss)*(d+pS) = E_ss*RF_ss*(RF+e)-E_ss*LFB_ss*(LFB+e)-V_ss*v; %(33)
(mP_ss*PS_ss)*(mP+pS)+(d_ss*PS_ss)*(d+pS) = E_ss*RF_ss*(RFT+e)-E_ss*LFB_ss*(LFB+e)-V_ss*v; %(33)

(1-psi)*(Y_ss/C_ss)*yS = c+(K_ss/C_ss)*(k(+1)-k)+delta*(K_ss/C_ss)*k; %(34)
hd = 0; %(35)
WPX*YX_ss*(wpX+yX)-WPF*YF_ss*(wpF+yF)+NFA_ss*iW_ss*iW(-1)+NFA_ss*(iW_ss-1)*NFA(-1) = thetaFB_ss*LFB_ss*thetaFB(-1)+(thetaFB_ss-1)*LFB_ss*LFB(-1)+NFA_ss*(NFA-NFA(-1)); %(36)

piH = pH(+1)-pH;
piD = pD(+1)-pD;
iB = iW+e(+1)-e;
NFA*NFA_ss = RFT*RFT_ss + BFP*BFP_ss - LFB*LFB_ss;
thetaFB = (theta_FB/2)*LFB_ss*LFB;
BFP = (iW_ss/(theta_FP*iW_ss))*iW - (iB_ss/(theta_FP*iW_ss))*iB;


%%%%%%%%%%var exo%%%%%%%%%%%%%
iW=rho_W*iW(-1)+eps_W; %Risk free world (37)

rhot = theta_C*(LI-LI(-1));


end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%Computation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
check;
steady;

shocks;
var eps_W;  stderr .01;
end;

stoch_simul(order=1,irf=50);