%*********************************************************************
%  Program        : 2.FMRegression_FF.m
%  Description    : Fama-MacBeth two-pass regressions of the
%                   Fama-French factor models augmented with the
%                   annual uncertainty shocks (Table 2 of the paper).
%                   First stage: full-sample betas for each asset.
%                   Second stage: cross-sectional regressions each
%                   year; the prices of risk are the time-series
%                   averages of the cross-sectional coefficients.
%                   Reported statistics: Fama-MacBeth t-statistics,
%                   Shanken (1992) EIV-corrected t-statistics, the
%                   R-squared as in Lettau/Ludvigson (1999) and
%                   Jagannathan/Wang (1996), and its bootstrapped
%                   [5% 95%] interval (10,000 simulations).
%
%  Input : Table2\ ff6_annual_disp_negativeprc.csv
%          (annualized Fama-French six factors and the uncertainty
%          shock; the last factor column is the uncertainty shock)
%          - ff6_annual_disp_negativeprc.csv  : total uncertainty
%            (UNC, d_tfpstd_star)
%          - ff6_annual_disp2_negativeprc.csv : macro uncertainty
%            (UNCma, d_pcahat)
%          - ff6_annual_disp3_negativeprc.csv : micro uncertainty
%            (UNCmi, d_pcares)
%          Table2\ testasset45y_negativeprc.csv
%          (45 annual test asset returns: 6 size/BM, 6 size/OP,
%          6 size/INV, 6 size/MOM, 6 size/EG, 10 operating accrual,
%          and 5 industry portfolios)
%  Helper: nw.m (Newey-West covariance matrix)
%          r2.m (R-squared for the bootstrap)
%*********************************************************************
diary 2.FMRegression_FF.txt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Read the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
projpath = 'D:\SAS\TFP\Dispersion_replication';
% Factor returns (T by K): T---# of obs, K---# of factors
% The default input below uses the total uncertainty shock (UNC).
% To reproduce the other panels of Table 2, change the file name:
%   ff6_annual_disp_negativeprc.csv  -> Panel A: total uncertainty (UNC, d_tfpstd_star)
%   ff6_annual_disp2_negativeprc.csv -> Panel B: macro uncertainty (UNCma, d_pcahat)
%   ff6_annual_disp3_negativeprc.csv -> Panel C: micro uncertainty (UNCmi, d_pcares)
%data = csvread(fullfile(projpath,'Table2','ff6_annual_disp_negativeprc.csv'),1,0);
%data = csvread(fullfile(projpath,'Table2','ff6_annual_disp2_negativeprc.csv'),1,0);
data = csvread(fullfile(projpath,'Table2','ff6_annual_disp3_negativeprc.csv'),1,0);
FF6_disp = data(1:45,2:8);
FF6 = FF6_disp(1:45,1:6);
FF3 = FF6_disp(1:45,1:3);
FF4 = FF6_disp(1:45,[1,2,3,6]);
FF5 = FF6_disp(1:45,1:5);
UNC = FF6_disp(1:45,7);
% Choose the base factor model: FF3, FF4, FF5, or FF6.
% When switching the model, also adjust the Shanken correction block in
% Section 5 (the number of s terms is the number of factors plus one).
FF = FF6;
rf = data(1:45,9);
Factors = [FF UNC]; %Delete UNC for FF3-FF6
Factors1 = [FF];
%
% Portfolio returns (T by N): T---# of obs, N---# of assets
data2 = csvread(fullfile(projpath,'Table2','testasset45y_negativeprc.csv'),1,0);
Rawreturn=data2(1:45,2:end);
R=Rawreturn-rf(1:45,ones(1,size(Rawreturn,2)))  ;
% The tested pricing factors are included as test assets.
R=[R Factors1];
%
fprintf('Average portfolio returns\n');
mean_port_ret=mean(Rawreturn)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Estimate full sample betas for each asset
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X  = [ones(size(Factors,1),1) Factors];

ST = [];

for i          = 1:size(R,2)
   beta(:,i)   = inv(X'*X)*X'*R(:,i);            % full sample betas
   error(:,i)  = R(:,i)-X*beta(:,i);
   covest      = size(X,1)*inv(X'*X)*nw(X.*kron(error(:,i),ones(1,size(X,2))),5)*inv(X'*X);    % adjust for autocorrelations
   stde        = sqrt(diag(covest));
   ST          = [ST stde];
   tstats(:,i) = beta(:,i)./stde;
   chr(i)      = var(X*beta(:,i))/var(R(:,i));
end

A = beta';
A = A(:,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. FM regressions: Run Cross-sectonal Regressions at each time period
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ret = R';

M  = mean(R);	% Time-series average for each portfolio return
V  = var(M);	% Cross-section variance of time-series averages of returns
CV = cov(R);

c  = zeros(size(X,2),size(Ret,2));
e  = zeros(size(R,2),size(Ret,2));

for i      = 1:size(Ret,2);
    c(:,i) = inv([ones(size(R,2),1) A ]'*...
                 [ones(size(R,2),1) A ])*...
                 [ones(size(R,2),1) A ]'*Ret(:,i);

    e(:,i) = Ret(:,i) - [ones(size(R,2),1) A ]*c(:,i);
end

z  = c';		                                             % Time-series of cross-sectional coefficients
x  = e';		                                             % Time-series of cross-sectional errors
m  = mean(z);	                                             % Time-series average of the cross-sectional coefficients
t  = (mean(z)./sqrt(var(z)))*sqrt(size(X,1));                % Fama-McBeth t-statistics
y  = mean(x);		                                         % Time-series average of the residuals for each portfolio
U  = var(y);                                                 % Cross-section variance of time-series averages of residuals

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4 R-squared as in Lettau/Ludvigson (1999), Jagannathan/Wang (1996)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gf = (V-U)/V;
gf = 1-(1-gf)*((size(R,2)-1)/((size(R,2))-size(c,1)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5 Shanken's EIV correction (formula taken from Cochrane, p.239), Shanken (1992) Theorem 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The default (FF6+UNC) uses s1-s8. For FF3+UNC use s1-s5; for FF4+UNC
% use s1-s6; for FF5+UNC use s1-s7.
F  = cov(X(:,2:end));
sh = (1+m(:,2:end)*inv(F)*m(:,2:end)'); % the multiplier for the adjustment
s1 = sqrt(sh)*(m(1)/t(1));   % s1-s8 are adjusted stds of the intercept and the factors
s2 = sqrt(sh*(((m(2)/t(2))^2)-(var(X(:,2))/size(X,1)))+(var(X(:,2))/size(X,1)));
s3 = sqrt(sh*(((m(3)/t(3))^2)-(var(X(:,3))/size(X,1)))+(var(X(:,3))/size(X,1)));
s4 = sqrt(sh*(((m(4)/t(4))^2)-(var(X(:,4))/size(X,1)))+(var(X(:,4))/size(X,1)));
s5 = sqrt(sh*(((m(5)/t(5))^2)-(var(X(:,5))/size(X,1)))+(var(X(:,5))/size(X,1)));
s6 = sqrt(sh*(((m(6)/t(6))^2)-(var(X(:,6))/size(X,1)))+(var(X(:,6))/size(X,1)));
s7 = sqrt(sh*(((m(7)/t(7))^2)-(var(X(:,7))/size(X,1)))+(var(X(:,7))/size(X,1)));
s8 = sqrt(sh*(((m(8)/t(8))^2)-(var(X(:,8))/size(X,1)))+(var(X(:,8))/size(X,1)));
%s9 = sqrt(sh*(((m(9)/t(9))^2)-(var(X(:,9))/size(X,1)))+(var(X(:,9))/size(X,1)));

s = [s1 s2 s3 s4 s5 s6 s7 s8];
newt = m./s;
disp('FM coefficients, t-stat, and Shanken adjusted t-stat: Intercept and factor risk premia');
[m;t;newt]
disp('R-squared as in Lettau/Ludvigson (1999), Jagannathan/Wang (1996)');
gf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6 compute the R2 interval, using bootstrap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rng(0,'twister')  % For reproducibility
bootstrapr2 = bootstrp(10000,@r2,Factors,R);
disp('The [5% 95%] interval of R-squared as in Lettau/Ludvigson (1999), Jagannathan/Wang (1996)');
[prctile(bootstrapr2,5) prctile(bootstrapr2,95)]
%%
diary off