% Master thesis - EL MAHI Gabriel
% Toulouse School of Economics
% 2023 - 2024

close all
clearvars
clc

% pour eviter des fonctions du meme nom
% it is the path where all functions are available
% restoredefaultpath
warning('off')
% addpath(genpath('routines'))
% addpath(genpath('/Applications/Dynare/4.5.7'))
%addpath('/Applications/Dynare/4.5.3/matlab')


%% %%%%%%%%%%%%%%%%%%%%%%%%% CALIBRATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Steady states
% Table 2
% par.y = 0.9231;
% par.y_F_star = 0.21*par.y;
% par.y_F = 0.25*par.y;
% par.y_H = 0.6983;
% par.y_F = 0.2308;
% par.y_i = 0.8922;
% par.x_i = 1.5463;
% par.c_h_H = 1.3570;
% par.c_l_H = 0.3445;
% par.c_h_F = 1.5975;
% par.c_l_F = 0.4055;
% par.i = 0.18*par.y;
% par.k = 3.95*par.y;
% par.d = 0.10*par.y;
% par.n_h_H = 0.62;
% par.n_l_H = 0.49;
% par.n_h_l_H = 0.19;
% par.n_h_F = 0.0048;
% par.n_l_F = 0.0025;
% par.l_h = 0.1152;
% par.l_l = 0.3875;
% par.u_h = 0.07;
% par.u_l = 0.12;
% par.r_k = 0.0986;
% par.w_h_H_H = 1.1619;
% par.w_l_H_H = 0.7746;
% par.w_hl_H_H = 0.8133;
% par.m_h_H = 0.0115;
% par.m_l_H = 0.0270;
% par.m_hl_H = 0.0159;
% par.v_h = 0.0165;
% par.v_l = 0.1116;
% par.psi_l_H_H = 0.3277;
% par.psi_hl_H_H = 0.90;
% par.psi_h_H_F = 0.70;
% par.psi_l_H_F = 0.45;
% par.psi_hl_H_F = 0.3090;
% par.p_H = 1.0417;
% par.p_F = 0.8481;
% par.e = 0.8481;
% 
% 
% 
% 
% par.r_d = 1.0530;
% par.r_star = 1.0530;
% par.z = 1.5498;
% par.b_z = 0.3001;
% par.phi_z = 1.8988;
% par.X_h = 0.6992;
% par.X_l = 1.4809;
% par.tau = 0.1874;
% par.x = 0.4612;
% par.s = 1 - 0.8353;
% par.O_h = 0.0256;
% par.O_l = 0.0033;
% 
% % table 1: parametrization
% % 1- Data / Simulation targets
% par.t_l = 0.69;
% par.t_h = 0.31;
% par.epsilon = 1.46;
% par.s_c = 0.71;
% par.s_g_c = 0.05;
% par.s_g_t = 0.08;
% par.kappa_h = 0.1;
% par.x_h_2 = 0.1;
% par.x_l_2 = 0.1;
% 
% % 2- Steady-state equations
% par.beta = 0.9497;
% par.delta = 0.0456;
% par.gamma = 1.8004;
% par.zeta = 0.7927;
% par.omega = 0.8142;
% par.q_h = 1.3316;
% par.Phi_h = 0.0166;
% par.Phi_l = 0.4327;
% par.kappa_l = 0.3163;
% par.omega_bar = 0.6444;
% par.mu_1 = 0.4129;
% par.mu_2 = 0.7288;
% par.b_1 = 0.1249;
% par.phi_2 = 1.4636;
% par.x_h_1 = 1.3161;
% par.x_l_1 = 3.2438;
% par.theta_h = 0.0555;
% par.theta_hl = 0.7841;
% par.theta_l = 0.7009;
% 
% % 3- Literature
% par.rho = 0.77;
% par.alpha = 0.47;
% par.phi_h = 1.7;
% par.phi_l = 1.7;
% par.eta = 1.01;
% par.sigma_h = 0.06;
% par.sigma_h_F = 0.06;
% par.sigma_l = 0.08;
% par.sigma_l_F = 0.08;
% par.Theta = 4;
% 
% % 4- Other
% par.phi_1 = 1;
% par.b_2 = 2;
% 
% 
% par.psi_h_H_H = 0.10;
% par.psi_h_F = 1.6 * par.psi_h_H_H;
% par.m_h_H = 0.0115;
% par.m_l_H = 0.0270;
% par.psi_l_F = 0.5;
% par.O_l = 0.0033;
% par.psi_l_H_H = 0.3277;
% par.psi_h_F = 0.70;
% par.O_h = 0.0256;
% 
% 
% 
% par.e_bar = 0.8481;
% par.gdp_bar = 1;



% par.w_h_F = 1;
% par.tau_h = 0.3;
% par.Lambda_c_h = 1;
% par.w_l_F = 1;
% par.tau_l = 0.3;
% par.psi_l_F = 1;
% par.omega_star = 1;
% par.gamma_star = 1;
% par.y_star = 2;
% par.tau = 1;
% par.psi = 1;
% par.gdp = 1;
% par.p_F = 1;
% par.g_t_h = 0.2;
% par.g_t_l = 0.2;
% par.P_star = 0.5;




par.t_h = 0.31;
par.t_l = 0.69;
par.mu_1 = 0.4129;
par.mu_2 = 0.7288;
% Unemployment benefits
par.omega_bar = 0.6444;
par.w_h_F = 1.1619;
par.w_l_F = 0.7746;
par.eta = 1.01;
par.delta = 0.0456;
par.Theta = 4;
par.beta = 0.9497;
par.sigma_h = 0.06;
par.sigma_l = 0.08;
par.sigma_h_F = 0.06;
par.sigma_l_F = 0.08;
par.psi_h_F = 0.70;
par.psi_l_F = 0.45;
par.Phi_h = 0.0166;
par.Phi_l = 0.4327;
par.phi_h = 1.7;
par.phi_l = 1.7;
par.X_l = 1.4809;
par.X_h = 0.6992;
par.theta_h = 0.0555;
par.theta_l = 0.7009;
par.theta_hl = 0.7841;
par.kappa_h = 0.1;
par.kappa_l = 0.3163;
par.alpha = 0.47;
par.epsilon = 1.46;
par.rho = 0.77;
par.zeta = 0.7927;
par.q_h = 1.3316;
par.e_bar = 0.8481;
par.gdp_bar = 0.9231;
par.d_bar = 0.10*0.9231;
par.omega = 0.8142;
par.gamma = 1.8004;
par.omega_star = par.omega;
par.gamma_star = par.gamma;
par.psi_rp = 1;

par.y_star = 1.5107;
par.tau = 0.1874;
par.phi_1 = 1;
par.phi_2 = 1.4636;
par.b_1 = 0.1249;
par.b_2 = 2;
par.p_F = 0.8481;

% pas calculé
par.P_star = 0.8481;


par.A = 1;
par.epsilon_rp = 0;
par.g_t_h_bar = 0.04;
par.g_t_l_bar = 0.04;


     
eval(['dynare model'])

% Nombre de variables endogènes
N = length(M_.endo_names);
% Liste des libellés des variables endogènes (Vecteur: Nx1)
pos = M_.endo_names;

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IRFs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loop over parameters
    par_name = 'beta';
    
    if exist('par_name') & strcmp(par_name,'beta')
        PAR = [0.94 0.95];
    else
        PAR = 1; 
    end


for iExp = 1:length(PAR)
       if exist('par_name')
           par.(par_name) = PAR(iExp);
       end

       
% fill M_.params with chosen parameter values
for i = 1:length(M_.param_names)
    M_.params(i) = par.(M_.param_names{i});
end

% Get policy functions
options_.order = 1; options_.nograph = 1;
[~,oo_,options_] = stoch_simul(M_,options_,oo_,M_.endo_names);

% Save IRFs in a structure
% for i=1:N
%     irf.(pos{i})(iExp,:) = oo_.irfs.([pos{i},'_shock_a']);
% end
% 
% for i=1:N
%     irf.(pos{i})(iExp,:) = oo_.irfs.([pos{i},'_shock_R']);
% end

for i=1:N
    irf.(pos{i})(iExp,:) = 5*oo_.irfs.([pos{i},'_shock_a'])/oo_.irfs.G_shock_g(1);
end

end
% 
% % plotting
% hFig = figure(1);
% set(gcf,'Color','w');
% set(hFig,'Units','normalized','Position',[0 0 0.8 0.8]);
% 
% % choose which variables to plot
% var = {'A', 'y', 'Inv', 'c_I', 'K', 'rr', 'G','Pi','unemp_W','unemp_B','w_W','w_B','theta_W','theta_B','R','wprem'};
% 
% varName.A = 'A';
% varName.y = 'output';
% varName.Inv = 'I';
% varName.c_I = 'c_I';
% varName.K = 'K';
% varName.rr = 'real interest rate';
% varName.G = 'G';
% varName.Pi = 'Inflation';
% varName.unemp_W = 'Unemployment W';
% varName.unemp_B = 'Unemployment B';
% varName.w_W = 'w_W';
% varName.w_B = 'w_B';
% varName.theta_W = 'theta_W';
% varName.theta_B = 'theta_B';
% varName.R = 'Nominal IR';
% varName.wprem = 'Wage premium';
% 
% 
% 
% 
% for j=1:length(var)
%     v=var{j}; %variable to be plotted
%     T = length(irf.(v)); % Lenght of irfs
% 
%     subplot(4,4,j)
%     hold on
%     plot(1:T,irf.(v),'LineWidth',3.5);
%     plot(1:T,zeros(1,T),'k','LineWidth',3);
%     hold off
%     set(gca,'Fontsize',12);
%     xlabel('quarter','Fontsize',18);
%     ylabel('percent','Fontsize',18);
%     title(varName.(v),'Fontsize',22);
%     grid minor;
%     axis tight;
% 
% end
% legend
% 
