clc; clear all; close all
%% Parameters that affects the steady state 
tau_k_star = 0.0;
tau_w_star = 0.0;
tau_d_star = 0.0;

bbeta_star = 0.99;
sigma_s_star = 0.8;
sigma_r_star = 0.5;
varphi_s_star =  1.5;
varphi_r_star =  2.5;
aalpha_star  =  0.33;  
lambda_star = 0.5;  
varepsilon_star = 6; 

chi_star =  0.7;  
zzeta_star = 52; % to ensure the steady state Nr_star = 1/3
phi_I_star  = 3;   
ddelta_star  =  0.025;  

% Monetary policy feedbacks
rho_r_star = 0.7;
theta_r_y_star  = 0.1;
theta_r_pi_star =  1.5;
% Fiscal policy feedbacks
rho_g_star =  0.75;
theta_g_bg_star = 0.1;
theta_g_y_star = 0.1;
Pi_star_bar  = (1.02)^(1/4); 
Pi_star = Pi_star_bar;
PiF_star = Pi_star_bar;
% Other Shock precess
rho_ms_star = 0.75;
rho_g_star  = 0.75;
rho_z_star  = 0.75;
rho_eta_star  = 0.75;
rho_m_star  = 0.75;
bg_star = 0.7;

calvo_star=0.6;               
gamma_p_star=(varepsilon_star-1)*calvo_star/(Pi_star^2*(1-calvo_star)*(1-bbeta_star*calvo_star)); 


%% Steady state
q_star = 1;
Y_star = 1;
eta_star = 1;
M_star = 1;
Pi_star = Pi_star_bar;
PiF_star = Pi_star_bar;
LAM_star = bbeta_star;
SX_star = 0;
X_star = 0;
R_star = Pi_star/bbeta_star;
r_star = R_star/Pi_star;
rk_star = (r_star + ddelta_star -1)/(1-tau_k_star);
mc_star = (varepsilon_star -1)/varepsilon_star + (1-bbeta_star)/varepsilon_star*gamma_p_star*(PiF_star -1)*PiF_star;
K_star = (aalpha_star*mc_star*Y_star)/rk_star;
N_star = 1/3;
Nr_star = (1/(zzeta_star*(1-chi_star)))^(1/(1+varphi_s_star));
w_star = ((1-aalpha_star)*Y_star*mc_star)/N_star;
Cr_star = w_star*(1-tau_w_star)*Nr_star;
Ns_star = N_star/(1-lambda_star) - (lambda_star)/(1-lambda_star)*Nr_star;
I_star = ddelta_star*K_star;
Gamma_star = (gamma_p_star/2)*(PiF_star -1)^2*Y_star;
d_star = (1-mc_star)*Y_star - Gamma_star;
d_s_star = lambda_star/(1-lambda_star)*d_star;
BGF_star = bg_star*Y_star;
BGF_s_star = lambda_star/(1-lambda_star)*BGF_star;
Cs_star = w_star*(1-tau_w_star)*Ns_star + (R_star/Pi_star -1)*BGF_s_star + (1-tau_d_star)*d_s_star;
C_star = lambda_star*Cr_star + (1-lambda_star)*Cs_star;
G_star = Y_star - C_star - I_star - Gamma_star;
J_star = G_star - w_star*N_star*tau_w_star - rk_star*K_star*tau_k_star -tau_d_star*d_star;
Z_star = (Y_star/(K_star^aalpha_star*N_star^(1-aalpha_star)))^(1/(1-aalpha_star));
Us_star = (((1-chi_star)*Cs_star)^(1-sigma_s_star))/(1-sigma_s_star)*exp(zzeta_star*((sigma_s_star -1)/(1+varphi_s_star))*Ns_star^(1+varphi_s_star));
Ur_star = (((1-chi_star)*Cr_star)^(1-sigma_r_star))/(1-sigma_r_star)*exp(zzeta_star*((sigma_r_star -1)/(1+varphi_s_star))*Nr_star^(1+varphi_r_star));
U_star = lambda_star*Ur_star + (1-lambda_star)*Us_star;
UCs_star = ((1-chi_star)*Cs_star)^(-sigma_s_star)*exp(zzeta_star*((sigma_s_star -1)/(1+varphi_s_star))*Ns_star^(1+varphi_s_star));
Vs_star = Us_star + bbeta_star*Us_star;
Vr_star = Ur_star + bbeta_star*Ur_star;
V_star = lambda_star*Vr_star + (1-lambda_star)*Vs_star;

%% Steady-state values (for policy part)
R_star_ss=R_star;
Pi_star_ss=Pi_star;
Y_star_ss=Y_star;
C_star_ss=C_star;
I_star_ss=I_star;
G_star_ss=G_star;
Z_star_ss=Z_star;  
BGF_star_ss=BGF_star;   
w_star_ss =w_star;
N_star_ss=N_star;
Nr_star_ss=Nr_star;
Ns_star_ss=Ns_star;


%% Save parameters
save param tau_k_star tau_w_star tau_d_star bbeta_star sigma_s_star sigma_r_star varphi_s_star varphi_r_star...
     aalpha_star lambda_star varepsilon_star chi_star zzeta_star phi_I_star  ddelta_star rho_r_star theta_r_y_star...
     theta_r_pi_star rho_g_star theta_g_bg_star theta_g_y_star Pi_star_bar Pi_star rho_ms_star rho_g_star...
     rho_z_star rho_eta_star rho_m_star bg_star calvo_star gamma_p_star G_star_ss Z_star_ss...
     Y_star_ss  BGF_star_ss  R_star_ss Pi_star_ss I_star_ss C_star_ss   w_star_ss  N_star_ss  Nr_star_ss  Ns_star_ss