% problem set 3 CAPM; 
% ps3.m has the full FF replication

clear all; 
close all; 

format bank; % two decimal points
format compact; 
% MAKE SURE you run the later subperiod last or it will overwrite graphs. 

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

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

if 0;
    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 st == 439;
    print -depsc2 ps3a.eps; 
end
if en == 438;
    print -depsc2 ps3a_early.eps; 
end
disp('CAPM betas'); 
disp(reshape(beta,5,5)'); 
figure; 
bar3(reshape(beta,5,5)')
ylabel('size'); 
xlabel('value'); 
zlabel('Capm beta'); 
if st == 439;
    print -depsc2 ps3b.eps;
end
if en == 438;
    print -depsc2 ps3b_early.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 st == 439;
    print -depsc2 capm_bar_plots_together.eps
end
if en == 438;
    print -depsc2 capm_bar_plots_together_early.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 st == 439;
    axis([ 0 1.5 0 1.3]);
    print -depsc2 capm_1.eps;
end;
if en == 438;
    axis([ 0 2 0 2.5]); 
    print -depsc2 capm_1_early.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 st == 439;
    axis([ 0 1 0 1.3]); 
    print -depsc2 capm_2.eps;
end;
if en == 438;
    axis([ 0 2.5 0 2.5]); 
    print -depsc2 capm_2_early.eps;
end;


