%Created: 19-Aug-2024 
%Author: Dennis Nsafoah

%Description: In this file, we calibrate parameters and solve for steady 
%state of the model.
%% Pre-determined parameters
%These parameters do not depend on the steady state of the model.
R_pol_para = 1/0.9925;
g_bar  = 0.2;
G_C = 1;
s_1 = 0.01; 
s_2 = 0.01; 
s_AC = 0.001;
s_AH = 0.001;
s_G = 0.01; 
s_J = 0.0416; 
s_R = 0.0025;
s_R_star = 0.01; 
s_Q_star = 0.01;
s_Y_F = 0.01; 
s_Z = 0.01; 
s_tau = 0.04; 
s_theta = 0.04; 
s_barphi = 0.04; 
Psi_1 = -7;
Psi_2 = -7;
alpha_C= 0.35;
alpha_L= 0.1;
alpha_H= 0.2;
alpha_N= 0.33; 
beta_ = 0.9925;
beta_prime = 0.97;
epsilon=0.32;
epsilon_prime=0.58;
Gamma_C= (G_C - epsilon)/(G_C - beta_*epsilon*G_C);
Gamma_C_prime = (G_C - epsilon_prime)/(G_C - beta_prime*epsilon_prime*G_C);
gamma_p = 0.469;
gamma_w = 0.763; 
delta_H = 0.01; 
delta0 = 0.025; 
delta2 = 0.01;
epsilon_p = 11;
epsilon_w = 11;
eta = 0.1; %0.52 is IN2010 estimate 
eta_prime =0.1; %0.51 is IN2010 estimate 
%xi = 0.66;
%xi_prime = 0.97;
eta1 = 1.5;
eta2 = 1.5;
theta_p = 0.83;
theta_wc = 0.79;
theta_wh = 0.91;
kappa = 39/40;
kappa_IC = 2;
kappa_IH = 2;
nu1 = 0.4;
rho_1 = 0.8; 
rho_2 = 0.8; 
rho_A = 0.0065;
rho_AC = 0.95;
rho_AH =  0.9997;
rho_G = 0.95; 
rho_J =  0.96;
rho_R =  0.8;
rho_R_star =  0.95;
rho_Q_star = 0.95;
rho_Y_F = 0.95; 
rho_Z =  0.96; 
rho_tau = 0.92; 
rho_theta =  0.95;
rho_barphi =  0.95;
barphi_s = 1;
phi_pi = 1.25;
phi_y = 0.25;
psi_ = 0.81;
m=0.85;

%% Targets
debt_to_GDP_ratio = 0.4054; 
ma_gov_bond_to_GDP_ratio  = 0.0334; 
ma_pvt_bond_to_GDP_ratio  = 0; %

LR = 4; % Leverage ratio = FI's Assets/Networth

% I define the next two targets as parameters for ease of notation                                             
sigma_gov = 0.0025;   %SS excess return on government bonds                                               
sigma_pvt = 0.0075;   %SS excess return on private bonds
sigma_B = 0.008;  %SS excess return on loans to impatient households
barN = 0.8;
barN_prime = 0.8;
barI = 0.232;
%% Steady state 
ACSS = 1;
AHSS = 1;
ZSS = 1;
%JSS = 1;
tauSS= 1;
uCSS = 1;
uHSS = 1;
pi_CSS = 1;
upsilon_wCSS =1;
upsilon_wHSS =1;
upsilon_pSS =1;
Q_starSS=1;
R_starSS=1.01;
barphiSS=1;
PhiSS=exp(1);
thetaSS=1;
R_polSS = R_pol_para;
R_DSS = 1/beta_*G_C;
R_govSS = R_DSS + sigma_gov;                                     
R_pvtSS = R_DSS + sigma_pvt; 
R_BSS = R_DSS + sigma_B; 
LambdaSS = beta_*G_C;
Lambda_primeSS = beta_prime*G_C;
Q_govSS = 1/(R_govSS-kappa);
Q_pvtSS = 1/(R_pvtSS-kappa);
R_SSS = R_DSS;
p_FSS = 1;
p_DSS = 1;
p_D_starSS = p_DSS;
p2SS = p_D_starSS*(epsilon_p -1)/epsilon_p;
p_KCSS = p_DSS;
M2SS = LambdaSS/( Q_pvtSS*(pi_CSS - kappa*LambdaSS) );
M1SS = 1 + psi_*(M2SS-1);
N_HSS = 1;
N_CSS = barN*N_HSS;
N_H_primeSS = 1;
N_C_primeSS = barN_prime*N_H_primeSS;
N2_CSS = (N_CSS + N_C_primeSS)/(upsilon_wCSS);
N2_HSS = (N_HSS + N_H_primeSS)/(upsilon_wHSS);
delta1 = 1/beta_*G_C - (1-delta0); % Target: uSS = 1 %This is a parameter
K_CSS = ( alpha_C*p2SS*N2_CSS^(1-alpha_C)/(delta1*p_KCSS*M1SS) )^(1/(1-alpha_C));
K_HSS = ( alpha_H*N2_HSS^(1-alpha_H-alpha_L)/(delta1*M1SS) )^(1/(1-alpha_H));%because qSS = p_KHSS
K_LSS = 1;
I_C_hatSS = delta0*K_CSS;
I_H_hatSS = delta0*K_HSS;
I_CSS = I_C_hatSS; 
I_HSS = I_H_hatSS;
GDPSS = (I_CSS+I_HSS)/barI; 
%qSS=(GDPSS-p_DSS*YSS)/H2SS;
qSS=(GDPSS-p_DSS*ACSS*(K_CSS)^alpha_C*N2_CSS^(1-alpha_C))/(AHSS*(K_HSS)^alpha_H*N2_HSS^(1-alpha_H-alpha_L)*K_LSS^alpha_L);
p_KHSS = qSS;
w2_CSS = (1-alpha_C)*(N2_CSS^(-alpha_C))*p2SS*(K_CSS^alpha_C); 
w2_HSS = (1-alpha_H-alpha_L)*(N2_HSS^(-alpha_H-alpha_L))*qSS*(K_HSS^alpha_H);  
w2_C_starSS = w2_CSS;
w2_H_starSS = w2_HSS;
w_CSS = w2_C_starSS*((epsilon_w -1)/epsilon_w);
w_HSS = w2_H_starSS*((epsilon_w -1)/epsilon_w);
f1_CSS = (w_CSS*(w2_CSS^(epsilon_w))*N2_CSS)/(1-theta_wc*LambdaSS); 
f1_HSS = (w_HSS*(w2_HSS^(epsilon_w))*N2_HSS)/(1-theta_wh*LambdaSS); 
f2_CSS = (w2_CSS^(epsilon_w)*N2_CSS)/(1-theta_wc*LambdaSS); 
f2_HSS = (w2_HSS^(epsilon_w)*N2_HSS)/(1-theta_wh*LambdaSS); 
Y2SS = ACSS*(K_CSS)^alpha_C*N2_CSS^(1-alpha_C);
H2SS = AHSS*(K_HSS)^alpha_H*N2_HSS^(1-alpha_H-alpha_L)*K_LSS^alpha_L;
YSS = Y2SS/upsilon_pSS;
G_bar_para = g_bar*GDPSS;
GSS = G_bar_para;
p_LSS = LambdaSS*(alpha_L*qSS*AHSS)*(K_HSS)^alpha_H*K_LSS^(alpha_L-1)*N2_HSS^(1-alpha_H-alpha_L);
xi = log(w_CSS/w_HSS)/log(N_CSS/N_HSS); %This is a parameter
xi_prime = log(w_CSS/w_HSS)/log(N_C_primeSS/N_H_primeSS); %This is a parameter
MU_NCSS = ZSS*tauSS*(1+eta)*(N_CSS^(xi))*(N_CSS^(1+xi)+N_HSS^(1+xi))^((eta-xi)/(1+xi));
MU_NHSS = ZSS*tauSS*(1+eta)*(N_HSS^(xi))*(N_CSS^(1+xi)+N_HSS^(1+xi))^((eta-xi)/(1+xi));
MU_CSS = MU_NCSS/w_CSS;
MU_NC_primeSS = ZSS*tauSS*(1+eta_prime)*(N_C_primeSS^(xi_prime))*(N_C_primeSS^(1+xi_prime)+N_H_primeSS^(1+xi_prime))^((eta_prime-xi_prime)/(1+xi_prime));
MU_NH_primeSS = ZSS*tauSS*(1+eta_prime)*(N_H_primeSS^(xi_prime))*(N_C_primeSS^(1+xi_prime)+N_H_primeSS^(1+xi_prime))^((eta_prime-xi_prime)/(1+xi_prime));
MU_C_primeSS = MU_NC_primeSS/w_CSS;
CSS = 1/MU_CSS;
C_FSS = nu1*(p_FSS^(-eta1))*CSS;
C_DSS = (1-nu1)*(p_DSS^(-eta1))*CSS;
C_primeSS = 1/MU_C_primeSS;
C_F_primeSS = nu1*(p_FSS^(-eta1))*C_primeSS;
C_D_primeSS = (1-nu1)*(p_DSS^(-eta1))*C_primeSS;
MU_HSS = MU_CSS*qSS - beta_*MU_CSS*qSS*(1-delta_H); 
XSS = YSS-C_DSS-C_D_primeSS-I_CSS-I_HSS -GSS;
lambda3SS =(MU_C_primeSS - MU_C_primeSS*R_BSS)/(R_BSS) ; 
MU_H_primeSS = MU_C_primeSS*qSS*delta_H - lambda3SS*m*qSS; 
Y_FSS=10*GDPSS;
nu2 = XSS/Y_FSS;%This is a parameter
b_pvtSS = psi_*(p_KCSS*I_CSS+p_KHSS*I_HSS)/((1-kappa)*Q_pvtSS); 
b_gov_ma_para = (ma_gov_bond_to_GDP_ratio*4*GDPSS)/Q_govSS; % This is a parameter
b_pvt_ma_para = (ma_pvt_bond_to_GDP_ratio*4*GDPSS)/Q_pvtSS; % This is a parameter
b_gov_maSS = b_gov_ma_para;
b_pvt_maSS = b_pvt_ma_para;
sSS = Q_pvtSS*b_pvt_maSS + Q_govSS*b_gov_maSS;
b_pvt_fiSS = b_pvtSS - b_pvt_maSS;
b_gov_bar = (debt_to_GDP_ratio*4*GDPSS)/Q_govSS; % This is a parameter
b_gov_fiSS = b_gov_bar - b_gov_maSS;
Pi_maSS = R_pvtSS*Q_pvtSS*b_pvt_maSS + ...
          R_govSS*Q_govSS*b_gov_maSS - (R_SSS)*sSS;
TSS = GSS + b_gov_bar*(1-(1-kappa)*Q_govSS) - Pi_maSS;
b_starSS = (p_DSS*XSS - p_FSS*(C_FSS+C_F_primeSS))/(Q_starSS-R_starSS*PhiSS*Q_starSS);
H_primeSS = (w_CSS*N_C_primeSS + w_HSS*N_H_primeSS-C_primeSS)/(qSS -(m*qSS/R_BSS)+m*qSS-qSS*(1-delta_H));
b_primeSS = m*((qSS*H_primeSS)/R_BSS);
HSS = (H2SS -delta_H*H_primeSS)/delta_H; 
Assets_fi = Q_pvtSS*b_pvt_fiSS + Q_govSS*b_gov_fiSS+ Q_starSS*b_starSS + sSS+b_primeSS;
nSS = Assets_fi/LR;
dSS = (LR-1)*nSS;
Delta = sigma_gov/sigma_pvt; %This is a parameter
bar_phi_a = (R_starSS*PhiSS - R_pvtSS)/(R_starSS*PhiSS*b_starSS);
phiSS = (Q_pvtSS*b_pvt_fiSS + Delta*Q_govSS*b_gov_fiSS +Q_starSS*b_starSS)/nSS;
OmegaSS = (thetaSS*phiSS*nSS*thetaSS -thetaSS*(R_BSS-R_DSS)*b_primeSS)/...
           (thetaSS*LambdaSS*(R_DSS*nSS-(R_starSS*bar_phi_a*PhiSS*Q_starSS*b_starSS^2))+thetaSS*phiSS*nSS*LambdaSS*sigma_pvt);
lambda1SS = OmegaSS*LambdaSS*sigma_pvt/(thetaSS- OmegaSS*LambdaSS*sigma_pvt);
vartheta = (1-OmegaSS)/(1-thetaSS*phiSS); %This is a parameter
chi = nSS*(1-vartheta*R_DSS) - ... %This is a parameter
    vartheta*( sigma_pvt*Q_pvtSS*b_pvt_fiSS + sigma_gov*Q_govSS*b_gov_fiSS +(R_starSS*PhiSS-R_DSS)*Q_starSS*b_starSS+sigma_B*b_primeSS);
divSS = CSS+qSS*HSS-w_CSS*N_CSS-w_HSS*N_HSS-qSS*(1-delta_H)*HSS+TSS - (R_DSS-1)*dSS +chi; 
x1SS = (p2SS*YSS)/(1-theta_p*LambdaSS);
x2SS = YSS/(1-theta_p*LambdaSS);
JSS = MU_H_primeSS*H_primeSS/ZSS;


%% Save the entire workspace in 'calibration.mat'
save calibration.mat
%% Save parameters as a separate file called 'params.mat'
save params.mat ...
R_pol_para g_bar G_C s_1 s_2 s_AC s_AH s_G s_J s_R s_R_star ... 
s_Q_star s_Y_F s_Z s_tau s_theta s_barphi Psi_1 Psi_2 alpha_C alpha_L alpha_H ... 
alpha_N beta_ beta_prime epsilon epsilon_prime G_bar_para Gamma_C Gamma_C_prime gamma_p ... 
gamma_w delta_H delta0 delta2 epsilon_p epsilon_w eta eta1 eta_prime xi xi_prime ... 
eta2 theta_p theta_wc theta_wh kappa kappa_IC kappa_IH nu1 rho_1 rho_2 rho_A rho_AC ... 
rho_AH rho_G rho_J rho_R rho_R_star rho_Q_star rho_Y_F rho_Z rho_tau rho_theta ... 
barphi_s phi_pi phi_y psi_ m b_gov_ma_para  b_pvt_ma_para b_gov_bar chi ...
Delta vartheta delta1 sigma_pvt sigma_gov sigma_B barN barI nu2 bar_phi_a qSS p_DSS b_starSS R_starSS...
R_DSS thetaSS Y_FSS Q_starSS R_starSS barphiSS rho_barphi JSS

   
%% Save steady-state values as a separate file called 'ss_values.mat'
save ss_values.mat ...
    ACSS AHSS ZSS tauSS uCSS uHSS pi_CSS upsilon_wCSS upsilon_wHSS upsilon_pSS ...
    Y_FSS Q_starSS R_starSS barphiSS PhiSS thetaSS R_polSS R_DSS R_govSS      ...                            
    R_pvtSS R_BSS LambdaSS Lambda_primeSS Q_govSS Q_pvtSS R_SSS p_FSS p_DSS ...
p_D_starSS p2SS p_KCSS M2SS M1SS MU_CSS MU_HSS qSS p_KHSS CSS C_FSS C_DSS ...
N_HSS N_CSS MU_NCSS MU_NHSS w_CSS w_HSS MU_NC_primeSS MU_NH_primeSS MU_C_primeSS ...
N_H_primeSS N_C_primeSS C_primeSS C_F_primeSS C_D_primeSS N2_CSS N2_HSS K_CSS ...
K_HSS K_LSS w2_CSS w2_HSS  w2_C_starSS w2_H_starSS f1_CSS f1_HSS f2_CSS f2_HSS ...
Y2SS H2SS YSS p_LSS I_C_hatSS I_H_hatSS I_CSS I_HSS GSS XSS lambda3SS MU_H_primeSS ...
JSS HSS H_primeSS b_primeSS GDPSS b_pvtSS b_gov_maSS b_pvt_maSS sSS  ...
b_pvt_fiSS b_gov_fiSS Pi_maSS TSS b_starSS nSS dSS ...
phiSS OmegaSS lambda1SS divSS x1SS x2SS
%% Create a vector called 'ys' that contains all steady-state values
ys = [ACSS 
AHSS 
ZSS 
tauSS
uCSS 
uHSS 
pi_CSS 
upsilon_wCSS 
upsilon_wHSS 
upsilon_pSS 
Y_FSS
Q_starSS
R_starSS
barphiSS
PhiSS
thetaSS
R_polSS 
R_DSS 
R_govSS                                  
R_pvtSS 
R_BSS 
LambdaSS 
Lambda_primeSS 
Q_govSS 
Q_pvtSS
R_SSS 
p_FSS 
p_DSS 
p_D_starSS 
p2SS
p_KCSS 
M2SS 
M1SS 
MU_CSS 
MU_HSS 
qSS 
p_KHSS 
CSS 
C_FSS 
C_DSS 
N_HSS 
N_CSS 
MU_NCSS 
MU_NHSS 
w_CSS 
w_HSS 
MU_NC_primeSS 
MU_NH_primeSS 
MU_C_primeSS 
N_H_primeSS 
N_C_primeSS 
C_primeSS 
C_F_primeSS
C_D_primeSS
N2_CSS 
N2_HSS 
K_CSS
K_HSS
K_LSS
w2_CSS 
w2_HSS  
w2_C_starSS 
w2_H_starSS 
f1_CSS 
f1_HSS
f2_CSS 
f2_HSS
Y2SS 
H2SS 
YSS 
p_LSS 
I_C_hatSS 
I_H_hatSS
I_CSS 
I_HSS 
GSS 
XSS 
lambda3SS
MU_H_primeSS
JSS 
HSS 
H_primeSS
b_primeSS 
GDPSS 
b_pvtSS 
b_gov_maSS
b_pvt_maSS 
sSS 
b_pvt_fiSS 
b_gov_fiSS 
Pi_maSS 
TSS 
b_starSS 
nSS 
dSS 
phiSS 
OmegaSS 
lambda1SS 
divSS
x1SS
x2SS
];
%% These commands count the number of parameters and variables in .mat files
clear
load calibration.mat
everything_vec = who;
disp(['There are ' num2str(length(everything_vec)) ' variables and parameters in calibration.mat.'])

clear
load params.mat
param_name_vec = who;
disp(['There are ' num2str(length(param_name_vec)) ' parameters in params.mat.'])

clear
load ss_values.mat
ssvar_name_vec = who;
disp(['There are ' num2str(length(ssvar_name_vec)) ' variables in ss_values.mat.'])

%% Some important Statistics
clear
load calibration.mat
Consumption = (CSS+C_primeSS)/GDPSS;
disp ([Consumption])




