%*********************************************************************
%  Program        : 4.FMRegression_extending.m
%  Description    : Extending-window Fama-MacBeth two-pass regressions
%                   (Table 4, Panels B and D of the paper).
%                   The mimicking-factor weights are re-estimated every
%                   year in SAS. This program estimates the
%                   extending-window betas (first stage) and the
%                   cross-sectional prices of risk (second stage).
%                   For each window, only the last twelve months of
%                   cross-sectional estimates are kept, so the gammas
%                   do not overlap across windows.
%                   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 : output\ FMR1996 .csv - FMR2016 .csv
%          (annually re-estimated mimicking factors, computed in SAS;
%          the exported file names contain a trailing blank before
%          ".csv")
%          Table4\extendingwindow\ testasset.csv
%          (45 monthly test asset returns in percent)
%  Helper: nw.m (Newey-West covariance matrix)
%          r2_extending.m (R-squared for the bootstrap; MATLAB
%          function names cannot start with a digit, so the helper
%          file has no table-number prefix)
%*********************************************************************
diary 4.FMRegression_extending.txt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Read the data: Projection method is extended every year.
% The weights of projection method is estimated every year, and they are computed in SAS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
projpath = 'D:\SAS\TFP\Dispersion_replication';
% Factor returns (T by K): T---# of obs, K---# of factors
% One file per estimation year; FMR{t} corresponds to the year 1995+t.
FMR = cell(1,21);
for t = 1:21
    FMR{t} = csvread(fullfile(projpath,'output',sprintf('FMR%d .csv',1995+t)),1,0);
end
% Matched extending-window return
data = csvread(fullfile(projpath,'Table4','extendingwindow','testasset.csv'),1,0);
Rawreturn2016=data(:,2:end)/100;
rf2016= FMR{21}(1:end,end);
R2016=Rawreturn2016-rf2016(:,ones(1,size(Rawreturn2016,2)))  ;

Factors2016=FMR{21}(1:end,2:6); %Change-E.g.,in HXZ (4 factors)+UNC, Factors2016=FMR{21}(1:end,2:6);
                                %       E.g.,in FF6 (6 factors)+UNC, Factors2016=FMR{21}(1:end,2:8);

% Split the test returns and factors for each extending-window.
    for t=1:20
    R_ext{t}=R2016(1:300+12*(t-1),1:end);
    Ret_ext{t}=R2016(1:300+12*(t-1),1:end)';
    Factors1 = Factors2016(1:300+12*(t-1),1:end); 
    R_ext{t}=[R_ext{t} Factors1];
    %returns in t+1 is used in the second stage in Fama-MacBeth
    F_ext{t}=[R2016(2:301+12*(t-1),1:end) Factors2016(2:301+12*(t-1),1:end)];
    Fut_ext{t}=[R2016(2:301+12*(t-1),1:end) Factors2016(2:301+12*(t-1),1:end)]';
    end
    % t=21, the last month observations are used in the second stage.
    R_ext{21}=R2016(1:299+12*(21-1),1:end);
    Ret_ext{21}=R2016(1:299+12*(21-1),1:end)';
    Factors1 = Factors2016(1:299+12*(21-1),1:end); 
    R_ext{21}=[R_ext{21} Factors1];
    F_ext{21}=[R2016(2:300+12*(21-1),1:end) Factors2016(2:300+12*(21-1),1:end)];
    Fut_ext{21}=[R2016(2:300+12*(21-1),1:end) Factors2016(2:300+12*(21-1),1:end)]';
% FM regressions: First stage-Estimate the extending window beta
for t=1:20
    
Factors = FMR{t}(1:300+12*(t-1),2:6); %Change-E.g.,in HXZ (4 factors)+UNC, Factors2016=FMR{21}(1:end,2:6);
                                      %       E.g.,in FF6 (6 factors)+UNC, Factors2016=FMR{21}(1:end,2:8);

rf = FMR{t}(1:300+12*(t-1),end);
R=R_ext{t};
pfr{t}=R(289+12*(t-1):300+12*(t-1),1:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %2. Estimate the extending 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);            % Extending-window 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))),6)*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);

%Save estimations
Ax{t}=X;
Sensitivity{t}=A;
Berror{t}=error;
clear error
end
% t=21
Factors = FMR{21}(1:299+12*(21-1),2:6); %Change-E.g.,in HXZ (4 factors)+UNC, Factors2016=FMR{21}(1:end,2:6);
                                        %       E.g.,in FF6 (6 factors)+UNC, Factors2016=FMR{21}(1:end,2:8);
rf = FMR{21}(1:299+12*(21-1),end);
R=R_ext{21};
pfr{21}=R(289+12*(21-1):299+12*(21-1),1:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %2. Estimate the extending 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);            % Extending-window 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))),6)*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);

%Save estimations
Ax{21}=X;
Sensitivity{21}=A;
Berror{21}=error;
clear error
% FM regressions: Second stage-Run Cross-sectonal Regressions at each time period
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for t = 1:20
c  = zeros(size(Ax{t},2),size(Fut_ext{t},2));
e  = zeros(size(pfr{t},2),size(Fut_ext{t},2));
for i      = 1:size(Fut_ext{t},2);
    c(:,i) = inv([ones(size(pfr{t},2),1) Sensitivity{t} ]'*...
                 [ones(size(pfr{t},2),1) Sensitivity{t} ])*...
                 [ones(size(pfr{t},2),1) Sensitivity{t} ]'*Fut_ext{t}(:,i);

    e(:,i) = Fut_ext{t}(:,i) - [ones(size(pfr{t},2),1) Sensitivity{t} ]*c(:,i);
    gamma{t}=c(1:end, end-11:end);
    epsilon{t}=e(1:end, end-11:end);
end
end
% t=21
c  = zeros(size(Ax{21},2),size(Fut_ext{21},2));
e  = zeros(size(pfr{21},2),size(Fut_ext{21},2));
for i      = 1:size(Fut_ext{21},2);
    c(:,i) = inv([ones(size(pfr{21},2),1) Sensitivity{21} ]'*...
                 [ones(size(pfr{21},2),1) Sensitivity{21} ])*...
                 [ones(size(pfr{21},2),1) Sensitivity{21} ]'*Fut_ext{21}(:,i);

    e(:,i) = Fut_ext{21}(:,i) - [ones(size(pfr{21},2),1) Sensitivity{21} ]*c(:,i);
    gamma{21}=c(1:end, end-10:end);
    epsilon{21}=e(1:end, end-10:end);
end
z  = [gamma{1} gamma{2} gamma{3} gamma{4} gamma{5} gamma{6} gamma{7} gamma{8}...
      gamma{9} gamma{10} gamma{11} gamma{12} gamma{13} gamma{14} gamma{15} gamma{16}...
      gamma{17} gamma{18} gamma{19} gamma{20} gamma{21} ]';	% Time-series of cross-sectional coefficients
m  = mean(z);	                                                                            % Time-series average of the cross-sectional coefficients
x  = [epsilon{1} epsilon{2} epsilon{3} epsilon{4} epsilon{5} epsilon{6} epsilon{7} epsilon{8}...
      epsilon{9} epsilon{10} epsilon{11} epsilon{12} epsilon{13} epsilon{14} epsilon{15} epsilon{16}...
      epsilon{17} epsilon{18} epsilon{19} epsilon{20} epsilon{21} ]';      % Time-series of cross-sectional errors

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
FF=F_ext{21};
R1996 = FF(288:539,1:end);                                   % Test asset and factor returns over the evaluation period
M1996 = mean(R1996);
V1996 = var(M1996);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4 R-squared as in Lettau/Ludvigson (1999), Jagannathan/Wang (1996)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gf = (V1996-U)/V1996;
gf = 1-(1-gf)*((size(R,2)-1)/((size(R,2))-size(z,2)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5 Shanken's EIV correction (formula taken from Cochrane, p.239), Shanken (1992) Theorem 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Factors2=FMR{21}(289:540,2:6);  %Change
Q = [ones(size(Factors2,1),1) Factors2];
F  = cov(Q(:,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-s6 are adjusted stds of the intercept and the factors
s2 = sqrt(sh*(((m(2)/t(2))^2)-(var(Q(:,2))/size(Q,1)))+(var(Q(:,2))/size(Q,1)));
s3 = sqrt(sh*(((m(3)/t(3))^2)-(var(Q(:,3))/size(Q,1)))+(var(Q(:,3))/size(Q,1)));
s4 = sqrt(sh*(((m(4)/t(4))^2)-(var(Q(:,4))/size(Q,1)))+(var(Q(:,4))/size(Q,1)));
s5 = sqrt(sh*(((m(5)/t(5))^2)-(var(Q(:,5))/size(Q,1)))+(var(Q(:,5))/size(Q,1)));
s6 = sqrt(sh*(((m(6)/t(6))^2)-(var(Q(:,6))/size(Q,1)))+(var(Q(:,6))/size(Q,1)));
%s7 = sqrt(sh*(((m(7)/t(7))^2)-(var(Q(:,7))/size(Q,1)))+(var(Q(:,7))/size(Q,1)));
%s8 = sqrt(sh*(((m(8)/t(8))^2)-(var(Q(:,8))/size(Q,1)))+(var(Q(:,8))/size(Q,1)));
%s9 = sqrt(sh*(((m(9)/t(9))^2)-(var(Q(:,9))/size(Q,1)))+(var(Q(:,9))/size(Q,1)));
s = [s1 s2 s3 s4 s5 s6 ];
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Factors_b{1}=FMR{1}(1:end, 2:end-1);
for t = 2:21
Factors_b{t}=FMR{t}(end-11:end, 2:end-1);
end
Factors_b2=[Factors_b{1};Factors_b{2};Factors_b{3};Factors_b{4};Factors_b{5};
            Factors_b{6};Factors_b{7};Factors_b{8};Factors_b{9};Factors_b{10};
            Factors_b{11};Factors_b{12};Factors_b{13};Factors_b{14};Factors_b{15};
            Factors_b{16};Factors_b{17};Factors_b{18};Factors_b{19};Factors_b{20};Factors_b{21}];
%
rng(0,'twister')  % For reproducibility
% The bootstrap is repeated 10,000 times as described in the paper.
bootstrapr2 = bootstrp(10000,@r2_extending,Factors_b2,R2016);
disp('The [5% 95%] interval of R-squared as in Lettau/Ludvigson (1999), Jagannathan/Wang (1996)');
[prctile(bootstrapr2,5) prctile(bootstrapr2,95)]
%%
 diary off