% problem set 4; follow up on Fama French lecture
% fama french simulation is at end
% remember to post tsregress for the students!

clear all; 
close all; 

format bank; % two decimal points
format compact; 

addpath('../Data');  % modify according to your directory structure. This is where to find regress_jc on my computer
addpath('../Programs');

sample = 2; 

if sample == 1;
    st = 67; % start date. 247 = 1947:01, 67 = 193201, 439 = 196301
    en = 438; % 0 : end of sample. Otherwise end date
end;

if sample == 2;
    st = 439; % start date. 247 = 1947:01, 67 = 193201, 439 = 196301
    en = 810; % 0 : end of sample. Otherwise end date
end

if sample == 3;
    st = 811; % start date.811 = 199401
    en = 0; % 0 : end of sample. Otherwise end date
end

if sample == 4;
    st = 439; % start date. 247 = 1947:01, 67 = 193201, 439 = 196301
    en = 0; % 0 : end of sample. Otherwise end date
end



% -----------------
% load data
% -----------------

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

T1 = size(x,1); 
if en == 0; 
    en = T1; 
end; 

x = x(st:en,:); 
r = x(:,2:end); % date is first column 

        % 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-1+size(r,1),:); % factors typically goes further than portfolios, so cut off  

en2 = size(r,1); 

rmrf = x(1:en2,2); 
smb = x(1:en2,3); 
hml = x(1:en2,4); 
rf  = x(1:en2,5);
dates = x(1:en2,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(any(isnan(r))); % 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))|any(any(isnan(x(2:end)))); 
    disp('you have missing values in the factors data'); 
end; 

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

rx = r-rf*ones(1,N);

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

f = [rmrf];
[alpha, beta, R2, siga, sigb] = tsregress(rx,f);


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

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 == 2;
    print -depsc2 ps3a.eps; 
end
if sample == 1;
    print -depsc2 ps3a_early.eps; 
end
if sample == 3;
    print -depsc2 ps3a_late.eps; 
end
disp('CAPM betas'); 
disp(reshape(beta,5,5)'); 
figure; 
bar3(reshape(beta,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('Capm beta'); 
if sample == 2;
    print -depsc2 ps3b.eps;
end
if sample == 1;
    print -depsc2 ps3b_early.eps;
end
if sample == 3;
    print -depsc2 ps3b_late.eps;
end


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)'); 


% bar graphs together
figure; 
subplot(1,2,1)
bar3(reshape(meanr,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %');
title('E(R) %') 

subplot(1,2,2)
bar3(reshape(beta,5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('beta'); 
title('Market Beta')

set(gcf,'paperposition',[0.25 2.5 8 3])
if sample == 2;
    print -depsc2 capm_bar_plots_together.eps
end
if sample == 1;
    print -depsc2 capm_bar_plots_together_early.eps
end
if sample == 3;
    print -depsc2 capm_bar_plots_together_late.eps
end

% er beta graphs
figure; 
plot(beta,mean(rx),'.','MarkerSize',10);
hold on
errorbar(beta,mean(rx),siga,siga,'.');
plot([0;2.5],[0;2.5]*mean(rmrf),'-r','LineWidth',2); 
plot(0,0,'sr','markerfacecolor','r');
plot(1,mean(rmrf),'sr','markerfacecolor','r');
text(0,0,'Rf');
text(1,mean(rmrf),'rmrf');
xlabel('\beta_i');
ylabel('E(R^{ei})')
text(beta,mean(rx),labels)
title('CAPM'); 
if st == 247; 
    axis([0 1.5 0.0 1.3]); 
end;
if sample == 2;
    axis([ 0 1.5 0 1.3]);
    print -depsc2 capm_1.eps;
end;
if sample == 1;
    axis([ 0 2 0 2.5]); 
    print -depsc2 capm_1_early.eps;
end;
if sample == 3;
    axis([ 0 2 0 2.5]); 
    print -depsc2 capm_1_late.eps;
end;

figure; 
plot(beta*mean(rmrf),mean(rx),'.','MarkerSize',10);
hold on
errorbar(beta*mean(rmrf),mean(rx),siga,siga,'.');
plot([0;2.5],[0;2.5],'-r','LineWidth',2); 
text(beta*mean(rmrf),mean(rx),labels)
xlabel('Predicted E(R^{ei}) = \beta_i*E(rmrf)');
ylabel('E(R^{ei})')
title('CAPM'); 
if st == 67; 
    axis([0.6 1.8 0.6 1.8]); 
end;
if st == 247; 
    axis([0 0.9 0 1.3]); 
end;
if sample == 2;
    axis([ 0 1 0 1.3]); 
    print -depsc2 capm_2.eps;
end;
if sample == 1;
    axis([ 0 2.5 0 2.5]); 
    print -depsc2 capm_2_early.eps;
end;
if sample == 3;
    axis([ 0 2.5 0 2.5]); 
    print -depsc2 capm_2_late.eps;
end;


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

f = [rmrf hml smb];
[alpha, beta, R2, siga, sigb] = tsregress(rx,f);

disp('3F model b'); 
disp(reshape(beta(:,1),5,5)'); 
figure; 
bar3(reshape(beta(:,1),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('3F b'); 
if sample == 2;
    print -depsc2 ps3c.eps;
end
if sample == 3;
    print -depsc2 ps3c_late.eps;
end


disp('3F model h'); 
disp(reshape(beta(:,2),5,5)'); 
figure; 
bar3(reshape(beta(:,2),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('3F h'); 
if sample == 2;
    print -depsc2 ps3d.eps;
end
if sample == 3;
    print -depsc2 ps3d_late.eps;
end


disp('3F model s'); 
disp(reshape(beta(:,3),5,5)'); 
figure; 
bar3(reshape(beta(:,3),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('3F s'); 
if sample == 2;
    print -depsc2 ps3e.eps;
end;
if sample == 3;
    print -depsc2 ps3e_late.eps;
end;

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


% bar graphs together

figure
subplot(2,2,1)
bar3(reshape(meanr,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %'); 
title('E(R) %'); 


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

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

subplot(2,2,4)
bar3(reshape(beta(:,3),5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('3F s'); 
title('SMB s');
if sample == 2; 
    print -depsc2 FF3F_bar_together.eps
end
if sample == 2; 
    print -depsc2 FF3F_bar_together_late.eps
end

% continue

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

figure; 
plot(beta*mean(f)',mean(rx),'.','MarkerSize',5);
hold on
errorbar(beta*mean(f)',mean(rx),siga,siga,'.');
plot([0;1.8],[0;1.8],'-r','LineWidth',2); 
text(beta*mean(f)',mean(rx),labels)
xlabel('b_i*E(rmrf) + h_i*E(hml) + s_i*E(smb)');
ylabel('E(R^{ei})')
title('FF3F model'); 
if st == 67;
    axis([0 1.8 0 1.8]); 
end;
if st == 247;
    axis([ 0 1.2 0 1.3]); 
end;
if sample == 2;
    axis([ 0 1.2 0 1.2]); 
    print -depsc2 ps3_ff3f_predEr.eps
end;
if sample == 3;
    axis([ 0 1.2 0 1.2]); 
    print -depsc2 ps3_ff3f_predEr_late.eps
end;

% -----------------------
% Premiums and dropping s?
% ----------------------

f = [rmrf hml smb];
fprintf('                  %8.8s %8.8s %8.8s \n', 'rmrf','hml','smb');
fprintf('mean (pct/mo)     %8.2f %8.2f %8.2f \n', mean(f));
fprintf('Sharpe            %8.2f %8.2f %8.2f \n', mean(f)./std(f));
fprintf('t                 %8.2f %8.2f %8.2f \n', size(f,1)^0.5*mean(f)./std(f));

[alpha, beta, R2, siga, sigb] = tsregress([hml smb],rmrf);
b = [ones(size(hml,1),1) rmrf]\ [hml smb];
e = [hml smb] - [ones(size(hml,1),1) rmrf]*b;   
    
disp('alpha and beta of on rmrf'); 
fprintf('alpha             %8.8s %8.2f %8.2f \n', '', alpha');
fprintf('alpha sharpe      %8.8s %8.2f %8.2f \n', '', alpha'./std(e));
fprintf('alpha t           %8.8s %8.2f %8.2f \n', '', alpha'./siga');
fprintf('beta              %8.8s %8.2f %8.2f \n', '', beta');
fprintf('beta t            %8.8s %8.2f %8.2f \n', '', beta'./sigb');

fprintf('               %8.8s %8.8s %8.8s \n','rmse(a)','mean|a|','R2');
f = [rmrf hml smb];
[alpha, beta, R2, siga, sigb] = tsregress(rx,f);
fprintf(' FF3F model    %8.3f %8.3f %8.3f \n',[mean(alpha.^2).^0.5 mean(abs(alpha))],mean(R2));

f = [rmrf hml ];
[alpha, beta, R2, siga, sigb] = tsregress(rx,f);
fprintf('  rmrf hml    %8.3f %8.3f %8.3f \n',[mean(alpha.^2).^0.5 mean(abs(alpha))],mean(R2));

f = [rmrf ];
[alpha, beta, R2, siga, sigb] = tsregress(rx,f);
fprintf(' rmrf only    %8.3f %8.3f %8.3f \n',[mean(alpha.^2).^0.5 mean(abs(alpha))],mean(R2));

w = inv(cov([rmrf hml smb]))*mean([rmrf hml smb])'*100;
f = [rmrf hml smb]*w; 
[alpha, beta, R2, siga, sigb] = tsregress(rx,f);
fprintf('1 factor ff   %8.3f %8.3f %8.3f \n',[mean(alpha.^2).^0.5 mean(abs(alpha))],mean(R2));
disp('weights'); 
disp(w');

w = inv(cov([rx]))*mean([rx])'*100;
f = rx*w; 
[alpha, beta, R2, siga, sigb] = tsregress(rx,f);
fprintf('1 factor mv  %8.3f %8.3f %8.3f \n',[mean(alpha.^2).^0.5 mean(abs(alpha))],mean(R2));
disp('weights'); 
disp(reshape(w,5,5));





% ******************************
% fama french simulation. What if there is no phenomenon?
% ******************************

stop

N = 25; % number of portfolios
T = 10000*12; % time span in months;

means = 1*(1:N)'/N; % mean excess return  per month.  
sd = 20^(1/12);

rets = sd*randn(T,N) + ones(T,1)*means'; 

disp('sample mean and sd returns')
disp(mean(rets))
disp(std(rets));

rm = mean(rets,2); % makes mean along horizontal dimension, so it's T x 1 average return
hml = (mean(rets(:,floor(2*N/3)+1:N),2) - mean(rets(:,1:floor(N/3)),2)); % makes 6-10 - 1:5 portfolio

% regressions on the market, hml, and both

[bv,sebv,R2v,R2vadj,v,F] = olsgmm(rets,[ones(T,1) rm],0,0);
disp('coefficients on constant and market; R2') 
fprintf('     ');
for i = 1:N; 
    fprintf(' %8i ',i); 
end;
fprintf('\n');
fprintf('const');
for i = 1:N; 
    fprintf(' %8.2f ',bv(1,i)); 
end;
fprintf('\n');
fprintf('mkt  ');
for i = 1:N; 
    fprintf(' %8.2f ',bv(2,i)); 
end;
fprintf('\n');
fprintf('R2   ');
for i = 1:N; 
    fprintf(' %8.2f ',R2v(i)); 
end;
fprintf('\n');

[bv,sebv,R2v,R2vadj,v,F] = olsgmm(rets,[ones(T,1) hml],0,0);
disp('coefficients on constant and hml; R2') 
for i = 1:N; 
    fprintf(' %8i ',i); 
end;
fprintf('\n');
fprintf('const');
for i = 1:N; 
    fprintf(' %8.2f ',bv(1,i)); 
end;
fprintf('\n');
fprintf('hml  ');
for i = 1:N; 
    fprintf(' %8.2f ',bv(2,i)); 
end;
fprintf('\n');
fprintf('R2   ');
for i = 1:N; 
    fprintf(' %8.2f ',R2v(i)); 
end;
fprintf('\n');

figure; 
plot((1:N)',means,'--b','linewidth',2);
hold on
plot((1:N)',bv(1,:)','-bv','linewidth',2);
plot((1:N)',bv(2,:)','--r*','linewidth',2);
plot((1:N)',R2v(:)','-k^','linewidth',2);
legend('mean','constant','hml','R2',4)
xlabel('portfolio');
%print -depsc2 ps3_simulation_hml.eps



[bv,sebv,R2v,R2vadj,v,F] = olsgmm(rets,[ones(T,1) rm hml],0,0);
disp('coefficients on constant market and hml; R2') 
for i = 1:N; 
    fprintf(' %8i ',i); 
end;
fprintf('\n');
fprintf('const');
for i = 1:N; 
    fprintf(' %8.2f ',bv(1,i)); 
end;
fprintf('\n');
fprintf('mkt  ');
for i = 1:N; 
    fprintf(' %8.2f ',bv(2,i)); 
end;
fprintf('\n');
fprintf('hml  ');
for i = 1:N; 
    fprintf(' %8.2f ',bv(3,i)); 
end;
fprintf('\n');

fprintf('R2   ');
for i = 1:N; 
    fprintf(' %8.2f ',R2v(i)); 
end;
fprintf('\n');

figure; 
plot((1:N)',means,'--b','linewidth',2);
hold on
plot((1:N)',bv(1,:)','-bv','linewidth',2);
plot((1:N)',bv(2,:)','-ro','linewidth',2);
plot((1:N)',bv(3,:)','--r*','linewidth',2);
plot((1:N)',R2v(:)','-k^','linewidth',2);
legend('mean','constant','market','hml','R2',4)
xlabel('portfolio');
print -depsc2 ps3_simulation.eps

