clear all
close all

addpath '..\Data'; 
addpath '..\Programs'; 
% ------------------------
% FF exercise
% eigenvalue decomposition of ff 25 returns
% ------------------------

st = 247; % start in 1947
%st = 67; % start 1932

x = load('25_Portfolios_5x5_monthly.txt'); 
T = size(x,1);
x = x(st:end,:); % start after NAs end
r = x(:,2:end); 

x = load('F-F_Research_Data_Factors_monthly.txt'); 
x = x(st:T,:);
rmrf = x(:,2);  
smb = x(:,3); 
hml = x(:,4); 
rf  = x(:,5);
dates = x(:,1); 
dates = floor(dates/100)+(dates-100*floor(dates/100))*1/12; % turns 192603 in to 1926+3/12; 
% dates not used,  wrds date bug in 05 data

rx = r-rf*ones(1,25); % Form ff excess returns. 

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

Sigma = cov(rx); 
[Q,L] = eig(Sigma); 
loads = Q*L^0.5; 

disp('Eigenvalue decomposition of FF data '); 
disp('Standard deviation of factors'); 
disp(100*diag(L)'.^0.5);

figure; 
plot(diag(L).^0.5,'-v')
title('standard deviation of factors in FF data'); 
print -depsc2 ps10ffa.eps;

figure; 
bar3(reshape(Q(:,end),5,5)')
ylabel('size'); 
xlabel('value'); 
title('First factor loadings'); 
print -depsc2 ps10q.eps; 

figure; 
subplot(2,2,1); 
bar3(reshape(Q(:,end),5,5)')
ylabel('size'); 
xlabel('value'); 
title('First factor loadings'); 

subplot(2,2,2); 
bar3(reshape(Q(:,end-1),5,5)')
ylabel('size'); 
xlabel('value'); 
title('Second factor loadings'); 

subplot(2,2,3); 
bar3(reshape(Q(:,end-2),5,5)')
ylabel('size'); 
xlabel('value'); 
title('Third factor loadings'); 

subplot(2,2,4); 
bar3(reshape(Q(:,end-3),5,5)')
ylabel('size'); 
xlabel('value'); 
title('Fourth factor loadings'); 
print -depsc2 ps10ffb.eps;

% factor models


figure; 
subplot(2,2,1); 
bar3(reshape(mean(rx)',5,5)')
ylabel('size'); 
xlabel('value'); 
title('mean FF 25 portfolio xs returns'); 

% one factor model. 
factors = rx*Q(:,end); % produces T x 1 factors time series
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK,u] = tsregress2(rx,factors);

subplot(2,2,2); 
bar3(reshape(alpha,5,5)')
ylabel('size'); 
xlabel('value'); 
title('\alpha, 1 factor model'); 

% three factor model
factors = rx*Q(:,end-2:end); % produces T x 3 factors time series
[alpha, beta, R23, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK,u] = tsregress2(rx,factors);

disp('betas 25 portfolios on 3 factors'); 
disp(beta); 
disp('columns of Q -- note same as betas'); 
disp(Q(:,end-2:end)); 

subplot(2,2,3); 
bar3(reshape(alpha,5,5)')
ylabel('size'); 
xlabel('value'); 
title('\alpha, 3 factor model'); 

disp('3 factor alphas'); 
disp(reshape(alpha,5,5)')
disp('rmse alpha') 
disp(mean(alpha.^2)^0.5); 
% 4 factor model

factors = rx*Q(:,end-3:end); % produces T x 3 factors time series
[alpha, beta, R24, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK,u] = tsregress2(rx,factors);

subplot(2,2,4); 
bar3(reshape(alpha,5,5)')
ylabel('size'); 
xlabel('value'); 
title('\alpha, 4 factor model'); 
print -depsc2 ps10ffc.eps;

disp('4 factor alphas'); 
disp(reshape(alpha,5,5)')
disp('rmse alpha') 
disp(mean(alpha.^2)^0.5); 

disp('3 factor r2'); 
disp(reshape(R23,5,5)')

disp('4 factor R2'); 
disp(reshape(R24,5,5)')


% fama french 

factors = [rmrf hml smb]; % produces T x 3 factors time series
[alpha, beta, R2FF, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK,u] = tsregress2(rx,factors);

subplot(2,2,4); 
bar3(reshape(alpha,5,5)')
ylabel('size'); 
xlabel('value'); 
title('\alpha, fama french factor model'); 
print -depsc2 ps10ff_ff.eps;

disp('FF alphas'); 
disp(reshape(alpha,5,5)')
disp('rmse alpha') 
disp(mean(alpha.^2)^0.5); 

disp('FF R2'); 
disp(reshape(R2FF,5,5)')


