%*********************************************************************
%  Program        : 3.FMRegression_residual.m
%  Description    : Fama-MacBeth two-pass regressions comparing total
%                   uncertainty and macro uncertainty with the
%                   residual uncertainty measures (Table 3, Panel B
%                   of the paper). The Hou et al. (2015) q-factor
%                   model is the base model, augmented with two
%                   uncertainty series at the same time.
%                   First stage: full-sample betas for each asset.
%                   Second stage: cross-sectional regressions each
%                   year. 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 : Table3\ q_annual_disp6_negativeprc.csv
%          (annualized q-factors and [d_tfpres, d_pcahat]: the
%          residual of total uncertainty and macro uncertainty)
%          -> Column (1) of Table 3, Panel B (default)
%          Table3\ q_annual_disp7_negativeprc.csv
%          (annualized q-factors and [d_tfpstd_star, d_hatres]:
%          total uncertainty and the residual of macro uncertainty)
%          -> Column (2) of Table 3, Panel B
%          Both files place the two uncertainty series in the same
%          columns, so only the file name needs to be changed.
%          Table3\ 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 3.FMRegression_residual.txt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Read the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
projpath = 'D:\SAS\TFP\Dispersion_replication';
% Factor returns (T by K): T---# of obs, K---# of factors
% Column (1): q_annual_disp6_negativeprc.csv (residual of total + macro)
% Column (2): q_annual_disp7_negativeprc.csv (total + residual of macro)
%data = csvread(fullfile(projpath,'Table3','q_annual_disp6_negativeprc.csv'),1,0);
data = csvread(fullfile(projpath,'Table3','q_annual_disp7_negativeprc.csv'),1,0);
FF6_disp = data(1:45,2:8);
% Columns 5:6 of FF6_disp are the two uncertainty series.
UNC = FF6_disp(1:45,5:6);
HXZ = FF6_disp(1:45,1:4);
rf = data(1:45,8);
Factors = [HXZ UNC];
Factors1 = [HXZ];
%
% Portfolio returns (T by N): T---# of obs, N---# of assets
data2 = csvread(fullfile(projpath,'Table3','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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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-s7 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)));

s = [s1 s2 s3 s4 s5 s6 s7 ];
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