% Estimate SDFMs and use max FEV share approach to replicate Moench and 
% Soofi-Siavash (2022) results for yield news, implied and realized stock 
% market volatility and business cycle news shocks. 
% For our application, we adapted some programs from James H. Stock and 
% Mark W. Watson (2016), kindly provided on Mark W. Watson's website. 
% We also used some of the FRED-Databases Matlab codes kindly provided by Michael W. McCracken.
% Emanuel Moench, e.moench@fs.de
% Soroosh Soofi-Siavash, ssoofi@lb.lt
% March 2022

clear
close all
rng(63287545);

addpath(genpath('m_utilities'));
addpath(genpath('data'));

%%----------- Sample Period and Calendars
firstdate_m = [1962 , 7];       % first date
lastdate_m  = [2019 , 6];       % last  date
[T_m,calvec_m,calds_m] = calendar_make(firstdate_m,lastdate_m,12);

firstdate_q = [1962 , 3];       % first date
lastdate_q  = [2019 , 2];       % last date
[T_q,calvec_q,calds_q] = calendar_make(firstdate_q,lastdate_q,4);

% ----------- Load Data
m2q=0; % 0, use the monthly data; 1, use the quarterly data
load_data=0; % 0, use existing file; 1, load the data from the Excel file
load_all_data;


% *********************************************************************
% Parameters set here are used in all the models 
% *********************************************************************
%- Sampling parameters
if m2q==0
    est_par.smpl_par.nfirst = firstdate_m;      % first date
    est_par.smpl_par.nlast  = lastdate_m;       % last date
    est_par.smpl_par.calvec = calvec_m;         % calendar
    est_par.smpl_par.nper   = 12;               % number of periods a year
else
    est_par.smpl_par.nfirst = firstdate_q;      % first date
    est_par.smpl_par.nlast  = lastdate_q;       % last date
    est_par.smpl_par.calvec = calvec_q;         % calendar
    est_par.smpl_par.nper   = 4;                % number of periods a year
end


%- VAR parameters
if m2q==0 
    p=2; % 2 when using the monthly data
else 
    p=1; % 1 when using the quarterly data
end
est_par.var_par.nlag   = p;  % number of lags
est_par.var_par.iconst = 1;  % include constant
est_par.var_par.icomp  = 1;  % compute Companion form of model 

%- Univariate AR parameters for errors
est_par.n_uarlag = p;  % Number of arlags for uniqueness

%- Factor analysis parameters
est_par.fac_par.catcode = catcode_vec; % vector of codes indicating how variables are used in the factor estimation -- 1 include in the yield curve factors estimation;2 include in the macro factors estimation; 3 otherwise
est_par.lambda.nt_min   = 40;          % min number of obs for any series used to estimate lamba, irfs, etc.

%- Number of bootstap simulations for computing SEs for IRFs and FEVDs
n_rep  = 500;          


% *********************************************************************
% Estimate the only-yields model and the IRFs and FEVDs for the three 
% yield shocks
% *********************************************************************
disp('---------------------------------------------------------------------');
disp('Estimate the only-yields model and the IRFs and FEVDs for the three yield shocks');

only_yields_sdfm_yield_shocks;
% Save results for future charts
fac_est_out_1=fac_est_out;
irf_fevd_out_1=irf_fevd_out;
se_irf_vdecomp_out_1=se_irf_vdecomp_out;


% *********************************************************************
% Estimate the macro-yields model and the IRFs and FEVDs for 
% the three yield shocks
% *********************************************************************
disp('---------------------------------------------------------------------');
disp('Estimate the macro-yields model and the IRFs and FEVDs for the three yield shocks');

macro_yields_sdfm_yield_shocks;
% Save results for future charts
fac_est_out_2=fac_est_out;
irf_fevd_out_2=irf_fevd_out;
se_irf_vdecomp_out_2=se_irf_vdecomp_out;


if m2q==0
    % *********************************************************************
    % Estimate the IRFs and FEVDs for the implied stock market volatility shock
    % *********************************************************************
    disp('---------------------------------------------------------------------');
    disp('Estimate the IRFs and FEVDs for the implied stock market volatility shock');
    
    macro_yields_sdfm_vxo_shock;
    % Save results for future charts
    irf_fevd_out_3=irf_fevd_out;
    se_irf_vdecomp_out_3=se_irf_vdecomp_out;  
    
    
    % *********************************************************************
    % Estimate the IRFs and FEVDs for the realized stock market volatility shock
    % *********************************************************************   
    disp('---------------------------------------------------------------------');
    disp('Estimate the IRFs and FEVDs for the realized stock market volatility shock');
    
    macro_yields_sdfm_rvol_shock;
    % Save results for future charts
    irf_fevd_out_4=irf_fevd_out;
    se_irf_vdecomp_out_4=se_irf_vdecomp_out;      
    
    
    % *********************************************************************
    % Estimate the IRFs and FEVDs for the business cycle news shock
    % *********************************************************************  
    disp('---------------------------------------------------------------------');
    disp('Estimate the IRFs and FEVDs for the business cycle news shock');   
    
    macro_yields_sdfm_business_cycle_news_shock;
    % Save results for future charts
    irf_fevd_out_5=irf_fevd_out;
    se_irf_vdecomp_out_5=se_irf_vdecomp_out;    
    
    
    % *********************************************************************
    % Estimate the IRFs and FEVDs for the residual yield news shock that
    % is orthogonal to the realized stock market volatility and business-
    % cycle news shocks
    % *********************************************************************  
    disp('---------------------------------------------------------------------');
    disp('Estimate the IRFs and FEVDs for the residual yield news shock'); 
    
    macro_yields_sdfm_residual_yield_news_shock;
    % Save results for future charts
    irf_fevd_out_6=irf_fevd_out;
    se_irf_vdecomp_out_6=se_irf_vdecomp_out; 
    
else 
    % *********************************************************************
    % Estimate the quarterly macro-yields model and the IRFs and FEVDs
    % for a André Kurmann and Christopher Otrok (2013) type slope shock
    % *********************************************************************   
    disp('---------------------------------------------------------------------');
    disp('Estimate the IRFs and FEVDs for a André Kurmann and Christopher Otrok (2013) type slope shock');
    
    quarterly_macro_yields_sdfm_ko_slope_shock;
    % Save results for future charts
    irf_fevd_out_7=irf_fevd_out;
    se_irf_vdecomp_out_7=se_irf_vdecomp_out;    
end


if m2q==0
    % *********************************************************************
    % Tabulate results from OLS regressions of one-year change in Treasury yields
    % on lagged factors and the correlations among the estimated shocks
    % *********************************************************************
    tabulate_OLS_regressions_correlations;
    
end

if m2q==0
	% ********************************************************************* 
    % Figures 1, 2, 3 and 4 -- IRFs and FEVDs for the three yield shocks
    % *********************************************************************     
    figures_yield_shocks;
    
    
	% ********************************************************************* 
    % Figures 6, 7 and 8 -- IRFs and FEVDs for the implied and realized 
    % volatility shocks
    % *********************************************************************     
    figures_rvol_vxo_shocks;
    
    
	% ********************************************************************* 
    % Figures 9, 10 and 11 -- IRFs and FEVDs for the business cycle news
    % shock
    % *********************************************************************      
    figures_business_cycle_news_shock;
    
    
	% ********************************************************************* 
    % Figures 5 and 12 -- Plots of the estimated shocks
    % *********************************************************************    
    figures_filtered_shock_series; 
    
    
    % *********************************************************************
    % Figure 13 -- IRFs and FEVDs of the ten-year yield and components for 
    % the US, Germany, the UK and Japan to a yield news shock that is 
    % orthogonal to the realized volatility and business cycle news shocks
    % ********************************************************************
    figures_purged_yield_news_shock;
else
    % *********************************************************************
    % Figure A.8 -- IRFs and FEVDs from quarterly macro-yields DFM for 
    % the Kurmann and Otrok slope shock
    % *********************************************************************
    figures_kurmann_otrok_slope_shock;
end
