% follow up on Fama French 

clear all; 
close all; 

addpath('..\data')
sample = 3; % 1:193201-196212 2:all  3:63-end

x = load('25_Portfolios_5x5_monthly.txt');

% -----------------
% load data
% -----------------
addpath('..\data'); 

x = load('25_Portfolios_5x5_monthly.txt'); 

if sample == 1; 
    st = 67; en = 438; % This used to be to 1967, adding 5*12. Why, was that better?
    % The FF model is really weird in the early sample. Small growth has a
    % huge beta! 
end;
if sample == 2; 
    st = 67; en = size(x,1)
end;
if sample == 3; 
    %st = 523; %197001 -- negative cross sectional slope is even stronger. 
    st = 439; 
    %st = 247; % choose 67 for beginning of sample, 247 for 1947 on
    en = size(x,1); 
end;

x = x(st:en,:); % start after NAs end

r = x(:,2:end); 

        % data organization 
        %                      Small                                 2                 ...
        %          Low     2      3      4    High    Low     2      3      4       
        %         
        
labels = { 's1v5' 's1v4' 's1v3' 's1v2' 's1v1' ... 
           's2v5' 's2v4' 's2v3' 's2v2' 's2v1' ... 
           's3v5' 's3v4' 's3v3' 's3v2' 's3v1' ...
           's4v5' 's4v4' 's4v3' 's4v2' 's4v1' ... 
           's5v5' 's5v4' 's5v3' 's5v2' 's5v1'};

x = load('F-F_Research_Data_Factors_monthly.txt'); 
x = x(st:st+size(r,1)-1,:); % allows factors to go further than portfolios  

rmrf = x(1:end,2); 
smb = x(1:end,3); 
hml = x(1:end,4); 
rf  = x(1:end,5);
dates = x(1:end,1); 
dates = floor(dates/100)+(dates-100*floor(dates/100))*1/12; % turns 192603 in to 1926+3/12; 

disp('Data sample'); 
disp([x(1,1) x(end,1)]);

if any(any(r < -98)); % any operates along one dimension only; twice tells you about a matrix
    disp('you have missing values in the return data'); 
end; 
if any(any(x(:,2:end)<-98)); 
    disp('you have missing values in the factors data'); 
end; 

T = length(rmrf);
N = size(r,2); 

% ----------------------------------------------
% run time series regression on the market alone
% ----------------------------------------------

rx = r-rf*ones(1,N);
f = [rmrf];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK, e] = tsregress2(rx,f);
sigma = cov(e); 

% ---------------------------
% display results
% ---------------------------

disp(' Sample');
disp([x(1,1) x(end,1)]);
disp(' '); 
disp('Time series regression results'); 
disp('As in FF all results are in boxes with size and book to market'); 
disp('mean return'); 
meanr = mean(rx)'; 
disp(reshape(meanr,5,5)'); 

figure; 
bar3(reshape(meanr,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %');

if sample == 3; 
    print -depsc2 ps3a.eps; 
end;
if sample == 1; 
    print -depsc2 ps3a_early.eps; 
end;



 
figure; 
bar3(reshape(beta,5,5)')
ylabel('size'); 
xlabel('value'); 
title('Capm beta'); 
if sample == 3; 
    print -depsc2 ps3b.eps;
end
if sample == 1; 
    print -depsc2 ps3b_early.eps; 
end;

figure
bar3(reshape(alpha,5,5)');
ylabel('size'); 
xlabel('value'); 
title('capm alphas %');
if sample == 3; 
    print -depsc2 ps3b_alpha.eps;
end
if sample == 1; 
    print -depsc2 ps3b_alpha_early.eps; 
end;

%set(gcf,'paperposition',[0.25 2.5 8 3]) use for two graphs side to side,
%subplot(1,2,1)

disp('CAPM betas'); 
disp(reshape(beta,5,5)');
disp('CAPM alphas'); 
disp(reshape(alpha,5,5)'); 
disp('T on CAPM alphas'); 
disp(reshape(alpha./siga,5,5)'); 
disp('CAPM R2'); 
disp(reshape(R2,5,5)'); 


fprintf('CAPM chi2 statistic, N, prob value (%%)   %7.2f %7.2f %7.2e \n' , [chi2stat N2 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)            %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, N, T-N-K, prob value (%%)    %7.2f %7.2f %7.2f %7.2e \n', [Fstat N2 TNK 100*Fpv]);
fprintf('10, 5, 1%% p values of F(N,T-N-K)         %7.2f %7.2f %7.2f \n', Fpvals);
fprintf('rms alpha, mean abs alpha,               %7.2f %7.2f  \n',mean(alpha.^2)^0.5,mean(abs(alpha))); 

% ------------------------------------------------
% ols cross sectional regression with free constant
% ------------------------------------------------

Erx = mean(rx)'; 
X = [ones(N,1) beta]; 
lambda2 = X\Erx;
alpha2 = Erx-X*lambda2; 
gamma2 = lambda2(1);
lambda2 = lambda2(2:end);

% not in problem set -- measures of cross-sectional R2 -- not all the same,
% Beware cross sectional R2!
%xsr21 = var(X*lambda2,1)/var(Erx,1); % var(beta*lambda)/var(Ere)
%xsr22 = 1 - var(alpha2,1)/var(Erx,1); % 1 - var(alpha)/var(E(re)
%xsr23 = 1 - (alpha2'*alpha2/size(alpha2,1))/var(mean(rx),1); % 1 - mean(alpha^2)/var(E(Re))

sigmaf = cov(rmrf);
K = 1; 
sigma_lambda2 = 1/T*(inv(X'*X)*X'*sigma*X*inv(X'*X)*(1+lambda2'*inv(sigmaf)*lambda2) + [0 zeros(1,K); zeros(K,1) sigmaf]);
sigma_lambda2 = diag(sigma_lambda2).^0.5; 
covalpha = 1/T*((eye(N)-X*inv(X'*X)*X')*sigma*(eye(N)-X*inv(X'*X)*X')'*(1+lambda2'*inv(sigmaf)*lambda2));
siga = diag(covalpha).^0.5; 
chi2stat2 = alpha2'*pinv(covalpha)*alpha2; 
        % This is a little tricky. cov(alpha) is a singular matrix 
        % The mean of the residual in any regression is always zero, so
        % there is a linear combination of alpha which does not vary from
        % sample to sample. Estimating lambda loses another degree of
        % freedom. For this reason, I use pinv the "pseudoinverse" and the
        % chi2 distribution has N-2 rather than N degrees of freedom. This
        % is a bit beyond the class and not in the problem set, but I might
        % as well do it right. 
disp('rank and size of covalpha. DOF should be rank of alpha'); 
disp([rank(covalpha) size(covalpha,1)]);
chi2vals2 = chi2inv([0.9 0.95 0.99],rank(covalpha)); 
chi2pv2 = 1-chi2cdf(chi2stat2,rank(covalpha)); 

% ------------------------------------------------
% ols cross sectional regression with free constant and factors included as
% portfolios 
% ------------------------------------------------

Erx = [mean(rx)'; 0; mean(rmrf)]; 
X = [ones(N+2,1) [beta;0;1]]; 
lambda2_a = X\Erx;
alpha2_a = Erx-X*lambda2_a; 
gamma_a = lambda2_a(1);
lambda2_a = lambda2_a(2:end);
sigma_a = [sigma zeros(N,2); zeros(2,N) zeros(2,2)]; 

sigmaf = cov(rmrf);
sigma_lambda_a = 1/T*(inv(X'*X)*X'*sigma_a*X*inv(X'*X)*(1+lambda2_a'*inv(sigmaf)*lambda2_a) + [0 zeros(1,K); zeros(K,1) sigmaf]);
sigma_lambda_a = diag(sigma_lambda_a).^0.5; 
covalpha_a = 1/T*((eye(N+2)-X*inv(X'*X)*X')*sigma_a*(eye(N+2)-X*inv(X'*X)*X')'*(1+lambda2_a'*inv(sigmaf)*lambda2_a));
siga_a = diag(covalpha_a).^0.5; 
chi2stat_a = alpha2_a'*pinv(covalpha_a)*alpha2_a; 
disp('rank and size of covalpha. DOF should be rank of alpha'); 
disp([rank(covalpha_a) size(covalpha_a,1)]);
chi2vals_a = chi2inv([0.9 0.95 0.99],rank(covalpha_a)); 
chi2pv_a = 1-chi2cdf(chi2stat_a,rank(covalpha_a)); 



% -------------------
% plot of Er vs beta and predicted
% -------------

figure; 
plot([0;1.8],[0;1.8]*mean(rmrf),'-r','LineWidth',2); 
hold on
plot([0;1.8], [gamma2+[0;1.8]*lambda2],'-vg','linewidth',2); % cross sectional slope 
plot([0;1.8], [gamma_a+[0;1.8]*lambda2_a],'-ob','linewidth',2); % cross sectional slope with factors too
plot(beta,mean(rx),'.','MarkerSize',10);
errorbar(beta,mean(rx),siga,siga,'.');
plot(0,0,'sr','markerfacecolor','r');
plot(1,mean(rmrf),'sr','markerfacecolor','r');
xlabel('\beta_i');
ylabel('E(R^{ei})')
legend('TS','CS, only 25','CS, 25+rmrf+rf'); 
text(beta,mean(rx),labels)
text(0,0,'Rf');
text(1,mean(rmrf),'rmrf');
if sample == 3; 
    axis([0 1.5 0.0 1.3]); 
end;
title('CAPM'); 
if sample == 3; 
    print -depsc2 xscapm.eps;
end;
if sample == 1; 
    print -depsc2 xscapm_early.eps;
end;


% actual vs. predicted, cross section
figure; 
errorbar(gamma2 + beta*lambda2,mean(rx)',siga,'o','MarkerFaceColor','b'); % plot actual means on right, and model predicted rhs. Error bars are std errors of time series alpha
hold on; 
plot([gamma2;gamma2 + lambda2], [0; mean(rmrf)],'rs', 'MarkerFaceColor','r'); % Add points for  the market and the risk free rate. 
hold on; 
text(gamma2,0,' Rf'); % identify Rm and Rf
text(gamma2 + lambda2,mean(rmrf),' Rm'); 
plot([0;2*mean(rmrf)], [0;2*mean(rmrf)]);  % add 45 degree line
text(gamma2+lambda2*beta,mean(rx),labels); 
ylabel('actual mean return E(R^{ei})'); 
xlabel('predicted mean return \gamma + beta_i x \lambda') 
title('Actual vs. predicted, cross-section');
if sample == 3; 
    print -depsc2 actual_predicted_capm.eps; 
end; 

disp(' '); 
disp(' comparison table'); 
fprintf('CAPM, Sample %8.0i to %8.0i \n', x(1,1), x(end,1));
fprintf('            %7.7s %7.7s %7.7s %7.7s %7.7s %7.7s %7.7s %7.7s %7.7s \n', 'gamma','lambda','s(gamma)','s(lambda)','rmse(a)','E|a|','chi2','5%','%p');
fprintf('TS         %7.7s %7.2f %7.7s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',' ',mean(rmrf),' ', std(rmrf)/T^0.5,mean(alpha.^2)^0.5,mean(abs(alpha)),chi2stat,chi2vals(2),100*chi2pv);
fprintf('CS, free g %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',gamma2,lambda2,sigma_lambda2,mean(alpha2.^2)^0.5,mean(abs(alpha2)),chi2stat2,chi2vals2(2),100*chi2pv2);
fprintf('CS, with f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',gamma_a,lambda2_a,sigma_lambda_a,mean(alpha2_a.^2)^0.5,mean(abs(alpha2_a)),chi2stat_a,chi2vals_a(2),100*chi2pv_a);


% ------------------------------------------------------
% Run time series regressions of the three factor model
% ------------------------------------------------------

f = [rmrf hml smb];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK] = tsregress2(rx,f);

disp('3F model b'); 
disp(reshape(beta(:,1),5,5)'); 
figure; 
subplot(2,2,1)
bar3(reshape(beta(:,1),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('3F b'); 

disp('3F model h'); 
disp(reshape(beta(:,2),5,5)'); 
subplot(2,2,2)
bar3(reshape(beta(:,2),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('3F h'); 
    
disp('3F model s'); 
disp(reshape(beta(:,3),5,5)'); 
subplot(2,2,3)
bar3(reshape(beta(:,3),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('3F s'); 

disp('3F model alphas'); 
disp(reshape(alpha,5,5)'); 
subplot(2,2,4)
bar3(reshape(alpha,5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('3F \alpha'); 
if sample == 3; 
    print -depsc2 ps3c.eps;
end

disp('T on 3F alphas'); 
disp(reshape(alpha./siga,5,5)'); 
disp('3F R2'); 
disp(reshape(R2,5,5)'); 


fprintf('chi2 statistic, N, prob value (%%)    %7.2f %7.2f %7.2e \n' , [chi2stat N2 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)        %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, N, T-N-K, prob value(%%) %7.2f %7.2f %7.2f %7.2e \n', [Fstat N2 TNK 100*Fpv]);
fprintf('10, 5, 1%% p values of F(N,T-N-K)     %7.2f %7.2f %7.2f \n', Fpvals);
fprintf('rms alpha,mean abs alpha,            %7.2f %7.2f  \n',mean(alpha.^2)^0.5,mean(abs(alpha))); 


figure; 
errorbar(beta*mean(f)',mean(rx),siga,'v');
hold on
%plot(beta*mean(f)',mean(rx),'v');
plot(mean(f)', mean(f)','sr','MarkerFaceColor','r'); 
plot([0 2],[0 2],'-r','linewidth',2)
xlabel('b_i*E(rmrf) + h_i*E(hml) + s_i*E(smb)');
text(mean(rmrf)', mean(rmrf)','rmrf'); 
text(mean(hml)', mean(hml)','hml'); 
text(mean(smb)', mean(smb)','smb'); 
text(beta*mean(f)',mean(rx),labels);

ylabel('E(R^{ei})')
axis([0 1.2 0 1.2]); 
title('FF3F model'); 
print -depsc2 xsff3f.eps;

% ----------------
% drop smb? 
% ---------------

f = [rmrf hml smb];
disp(' '); 
disp('testing for dropping factors'); 
disp('sample ' )
disp([x(1,1) x(end,1)]);
fprintf('         %8.8s %8.8s %8.8s \n','rmrf','hml','smb'); 
fprintf('E(f)     %8.3f %8.3f %8.3f \n',mean(f))
fprintf('t        %8.3f %8.3f %8.3f \n',mean(f)./std(f)*T^0.5);
fprintf('sharpe   %8.3f %8.3f %8.3f \n',(mean(f)./std(f))*12^0.5);



f = [rmrf smb ];
[alphah, beta, R2, sigah, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK,eh] = tsregress2(hml,f);
f = [rmrf hml ];
[alphas, beta, R2, sigas, sigh, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK,es] = tsregress2(smb,f);
f = [hml smb ];
[alpham, beta, R2, sigam, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK,em] = tsregress2(rmrf,f);

fprintf('E(f*)    %8.3f %8.3f %8.3f \n',[alpham alphah alphas])
fprintf('t        %8.3f %8.3f %8.3f \n',[alpham alphah alphas]./[sigam sigah sigas]);
fprintf('sharpe   %8.3f %8.3f %8.3f \n',[alpham/std(em) alphah/std(eh) alphas/std(es)]*12^0.5);
f = [rmrf hml smb];


% ------------------------
% Eigenvalues
% ------------------------

[Q,L]= eig(cov(rx));
disp('eigenvalues')

figure;
plot((diag(L).^0.5),'-v');
title('eigenvalues')
ylabel('\lambda^{0.5}'); 
if sample ==3; 
    print -depsc2 lambdas.eps
end

figure; 
subplot(2,2,1);
bar3(reshape(Q(:,end),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('Loading 1'); 

subplot(2,2,2);
bar3(reshape(Q(:,end-1),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('Loading 2'); 


subplot(2,2,3);
bar3(reshape(Q(:,end-2),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('Loading 3'); 

subplot(2,2,4);
bar3(reshape(Q(:,end-3),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('Loading 4'); 

if sample == 3; 
    print -depsc2 eig_1.eps;
end;

fe = rx*Q;
figure;
plot(abs(mean(fe)),'-v'); 
hold on
plot(abs(mean(fe))./std(fe)*size(fe,1)^0.5,'-or');
legend('abs(means)','t stats',2)
title('means of factors')
if sample == 3; 
    print -depsc2 mean_t_eigs.eps;
end;


disp('compare FF3F and eigenvalue models') 

fprintf('            %8.7s %8.7s %8.7s %8.7s %8.7s %8.7s %8.7s \n','chi2','N','%pv','5%','rmsa','|a|','R2')
[alphaff, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK, e] = tsregress2(rx,f);
fprintf(' FF3F       %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f  \n' , chi2stat, N2, 100*chi2pv,chi2vals(2),mean(alphaff.^2)^0.5,mean(abs(alphaff)),mean(R2));
disp('3 FF alphas'); 
disp(reshape(alphaff,5,5)')

figure; 
subplot(2,2,1); 
errorbar(beta*mean(f)',mean(rx)',siga,'v');
hold on
plot([0 1.2],[0 1.2],'-r');
plot(mean(f)',mean(f),'s','markerfacecolor','r')
title('FF3F')
axis([0 1.2 0 1.2]); 
xlabel('predicted')
ylabel('actual'); 

[alpha2, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK, e] = tsregress2(rx,fe(:,end:-1:end-1));
fprintf(' eig F2     %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n' , chi2stat, N2, 100*chi2pv,chi2vals(2),mean(alpha2.^2)^0.5,mean(abs(alpha2)),mean(R2));
disp('2 factor alphas'); 
disp(reshape(alpha2,5,5)')

subplot(2,2,2); 
errorbar(beta*mean(fe(:,end:-1:end-1))',mean(rx)',siga,'v');
hold on
plot([0 1.2],[0 1.2],'-r');
plot(abs(mean(fe(:,end:-1:end-1))')/4,abs(mean(fe(:,end:-1:end-1))')/4,'s','markerfacecolor','r')
title('2 eigenvalue') 
axis([0 1.2 0 1.2]); 
xlabel('predicted')
ylabel('actual'); 


[alpha3, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK, e] = tsregress2(rx,fe(:,end:-1:end-2));
fprintf(' eig F3     %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n' , chi2stat, N2, 100*chi2pv,chi2vals(2),mean(alpha3.^2)^0.5,mean(abs(alpha3)),mean(R2));
disp('3 factor alphas'); 
disp(reshape(alpha3,5,5)')

subplot(2,2,3); 
errorbar(beta*mean(fe(:,end:-1:end-2))',mean(rx)',siga,'v');
hold on
plot([0 1.2],[0 1.2],'-r');
plot(abs(mean(fe(:,end:-1:end-2))')/4,abs(mean(fe(:,end:-1:end-2))')/4,'s','markerfacecolor','r')
title('3 eigenvalue') 
axis([0 1.2 0 1.2]); 
xlabel('predicted')
ylabel('actual'); 

[alpha4, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK, e] = tsregress2(rx,fe(:,end:-1:end-3));
fprintf(' eig F4     %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n' , chi2stat, N2, 100*chi2pv,chi2vals(2),mean(alpha4.^2)^0.5,mean(abs(alpha4)),mean(R2));
disp('4 factor alphas'); 
disp(reshape(alpha4,5,5)')

subplot(2,2,4); 
errorbar(beta*mean(fe(:,end:-1:end-3))',mean(rx)',siga,'v');
hold on
plot([0 1.2],[0 1.2],'-r');
plot(abs(mean(fe(:,end:-1:end-3))')/4,abs(mean(fe(:,end:-1:end-3))')/4,'s','markerfacecolor','r')
title('4 eigenvalue') 
axis([0 1.2 0 1.2]); 
xlabel('predicted')
ylabel('actual'); 


if sample == 3; 
    print -depsc2 er_vs_pred_eigs.eps
end;

figure
subplot(2,2,1)
bar3(reshape(alphaff,5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('\alpha');
title('Fama-French')

subplot(2,2,2)
bar3(reshape(alpha2,5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('\alpha');
title('2 eigenvalue')

subplot(2,2,3)
bar3(reshape(alpha3,5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('\alpha');
title('3 eigenvalue')

subplot(2,2,4)
bar3(reshape(alpha4,5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('\alpha');
title('4 eigenvalue')

print -depsc2 eigalphas.eps

