% problem set 5; follow up on 2nd Fama French lecture; dissecting
% anomalies, cross sectional regressions. 

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 = 439; % start date. 247 = 1947:01, 67 = 193201, 439 = 196301
en = 0; % 0 : end of sample. Otherwise end date
end

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

if 0;
st = 67; % start date. 247 = 1947:01, 67 = 193201, 439 = 196301
en = 438; % 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); 

        % 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 goes further than portfolios  

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

% ----------------------
% be/me data, cross sectional regressions
% -------------------------

x = load('25_portfolios_5x5_BEME.txt'); 
if st == 67; 
    beme = x(6:end-1,2:end); % be/me of the 25 portfolios, starting in 1926. Start in 1932 to avoid -99 This timing uses year t-1 beme for all returns in year t
end;
if st == 247; 
    beme = x(21:end-1,2:end); % be/me of the 25 portfolios, starting in 1926. Start in 1932 to avoid -99 This timing uses year t-1 beme for all returns in year t
end;
if st == 439; 
    beme = x(37:end-1,2:end); % be/me of the 25 portfolios, starting in 1926. Start in 1932 to avoid -99 This timing uses year t-1 beme for all returns in year t
end;

x = load('25_portfolios_5x5_size.txt');
mc = x(st-1:end-1,2:end); % shifted one back from returns, so you can forecast

% expand beme to annual; note two extra months of other variables. 
beme = kron(beme,ones(12,1)); 

T = min([size(r,1) size(mc,1) size(beme,1)]');


beme = beme(1:T-1,:);% ) 
r = r(2:T,:);
rx = rx(2:T,:);
mc = mc(1:T-1,:);
T = T-1; % now represents sample size

% r and rx [1:T] are now next month's return; me and beme [1:T] are this
% month's values, ready for forecasting regressions

% plots at the end of the solutions to explain units choices

figure; 
subplot(2,2,1)
bar3(reshape(mean(beme),5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('BE/ME');
title('BE/ME')

subplot(2,2,2)
bar3(reshape(mean(log(beme)),5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('log(BE/ME)'); 
title('log(BE/ME)')

subplot(2,2,3)
bar3(reshape(mean(mc),5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('size'); 
title('Size')

subplot(2,2,4)
bar3(reshape(mean(log(mc)),5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('log(size)'); 
title('log(Size)');

print -depsc2 ps3_units.eps; 

if st == 67; 
    beg = 1932;
end

if st == 247; 
    beg = 1947;
end
if st == 439; 
    beg = 1963;
end

last = beg+1/12*size(mc,1);

%logmc2 = log(mc) - (mean(log(mc)')'*ones(1,25)); %-- old, take out cross sectional average
 logmc2 = log(mc./(sum(mc,2)*ones(1,25))); 
% A better way, size relative to
% total market cap. The diference is whether 1/25 goes in or out of the
% log. It seems to make no difference at all
% except intercept. 

figure;
subplot(2,2,1)
plot((beg+1/12:1/12:last)',log(mc));
title('log size over time'); 

 subplot(2,2,2);
plot((beg+1/12:1/12:last)',logmc2);
title('log size / total market value at each point in  time'); 

lastbeme = beg+1/12*size(beme,1);

subplot(2,2,3)
plot((beg+1/12:1/12:lastbeme)',log(beme));
title('log beme over time'); 
print -depsc2 ps3_size_beme_over_time.eps

% -----------------------------
% ols regression of average return on average beme, size
% -----------------------------

y = mean(rx)'; 
x = [ones(25,1) mean(logmc2)' mean(log(beme))']; 
[b,bse,bt,r2,F] = regress_jc(y,x,0);
fitmeans = x*b; 

fprintf('regression of mean xs return on mean log beme and mean log size\n'); 
fprintf('       %8.8s %8.8s %8.8s %8.8s \n','const','size','beme','R2'); 
fprintf('    b  %8.3f %8.3f %8.3f  %8.3f \n',b,var(fitmeans)/var(y)); 


disp('') ;
disp('actual mean return'); 
disp(reshape(mean(rx)',5,5)'); 
disp('regression on average log size and average log beme'); 
disp(reshape(fitmeans,5,5)'); 

figure; 
subplot(1,2,1)
bar3(reshape(fitmeans,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %'); 
title('mean returns from cs regression'); 
%axis([-inf inf -inf inf -0.2 2.5]);

subplot(1,2,2)
bar3(reshape(y',5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %'); 
title('mean returns from mean beme,size'); 
%axis([-inf inf -inf inf -0.2 2.5]);
print -depsc2 ps3_xs.eps; 

figure; 
plot(fitmeans,(y)','.','MarkerSize',10);
hold on
plot([0 2]',[0 2]','-r');
xlabel('Fitted mean a + b*E(size) + c*E(bm)')
ylabel('Actual mean')
if st == 247; 
    axis([0 1.2 0 1.2]); 
end
if st == 67; 
    axis([0 1.6 0 1.6]); 
end
if st == 439; 
    axis([0 1.3 0 1.3]); 
end
text(fitmeans,(y)',labels);
print -depsc2 fitmean2.eps

% ----------------------------
% cross term to improve fit
% -----------------------------

y = mean(rx)'; 
x = [ones(25,1) mean(logmc2)' mean(log(beme))' mean(logmc2)'.*mean(log(beme))' ]; 
[b,bse,bt,r2,F] = regress_jc(y,x,0);
fitmeans = x*b; 

fprintf('regression of mean xs return on mean log beme and mean log size and cross term \n'); 
fprintf('       %8.8s %8.8s %8.8s %8.8s %8.8s \n','const','size','beme','s x bm','R2'); 
fprintf('    b  %8.3f %8.3f %8.3f %8.3f %8.3f \n',b,var(fitmeans)/var(y)); 

figure; 
subplot(1,2,1)
bar3(reshape(fitmeans,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %'); 
title('mean returns from cs regression'); 
%axis([-inf inf -inf inf -0.2 2.5]);

subplot(1,2,2)
bar3(reshape(y',5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %'); 
title('mean returns from mean beme,size, and cross term'); 
%axis([-inf inf -inf inf -0.2 2.5]);
print -depsc2 ps3_xs_crossterm.eps; 

figure; 
plot(fitmeans,(y)','.','MarkerSize',10);
hold on
plot([0 2]',[0 2]','-r');
xlabel('Fitted mean a + b*E(size) + c*E(bm) + d*E(s)*E(bm)')
ylabel('Actual mean')
if st == 247; 
    axis([0 1.2 0 1.2]); 
end
if st == 67; 
    axis([0 1.6 0 1.6]); 
end
if st == 439; 
    axis([0 1.3 0 1.3]); 
end
text(fitmeans,(y)',labels);
print -depsc2 fitmean2_crossterm.eps


% -----------------------------
% pooled time-series cross section forecast
% -----------------------------

y = reshape(rx,T*25,1); 
x = [ones(T*25,1) reshape(logmc2,25*T,1) reshape(log(beme),T*25,1)]; 
[b_p,bse_p,bt_p,r2_p,F_p] = regress_jc(y,x,0);
fit_p = x*b_p; 
r2 = var(fit_p)/var(y);
Erx_p = reshape(fit_p,T,25); 

fprintf('pooled forecasting regression of return on log beme and log size\n');
fprintf('    b  %8.3f %8.3f %8.3f %8.3f \n',b_p,r2); 

% ------------------------------------------------------------
% forecast using average and differences from averages
% ---------------------------------------------------------

%  the portfolio means themselves

logmc4 = ones(T,1)*mean(logmc2);
logbeme4 = ones(T,1)*mean(log(beme));
x = [ones(T*25,1)  reshape(logmc4,25*T,1) reshape(logbeme4,T*25,1)]; 
[b_ppd,bpe_ppd,bt_ppd,r2_ppd,F_ppd] = regress_jc(y,x,0);
fprintf('pooled forecasting regression of sample mean x for each portfolio\n');
fprintf('     b %8.3f %8.3f %8.3f \n',b_ppd); 

logmc4 = logmc2 - ones(T,1)*mean(logmc2);
logbeme4 = log(beme) - ones(T,1)*mean(log(beme));
x = [ones(T*25,1)  reshape(logmc4,25*T,1) reshape(logbeme4,T*25,1)]; 
[b_ppd,bpe_ppd,bt_ppd,r2_ppd,F_ppd] = regress_jc(y,x,0);
fprintf('pooled forecasting regression of x - sample mean x for each portfolio\n');
fprintf('     b %8.3f %8.3f %8.3f \n',b_ppd); 

% taking out time-means

logmc3 = mean(logmc2,2)*ones(1,25);
logbeme3 = mean(log(beme),2)*ones(1,25);
x = [ ones(T*25,1) reshape(logmc3,25*T,1) reshape(logbeme3,T*25,1)]; 
[b_ptd,bse_ptd,bt_ptd,r2_ptd,F_ptd] = regress_jc(y,x,0);
fprintf('pooled forecasting regression, time means \n');
fprintf('    b  %8.3f %8.3f %8.3f \n',b_ptd); 

logmc3 = logmc2 - mean(logmc2,2)*ones(1,25);
logbeme3 = log(beme) - mean(log(beme),2)*ones(1,25);
x = [ ones(T*25,1) reshape(logmc3,25*T,1) reshape(logbeme3,T*25,1)]; 
[b_ptd,bse_ptd,bt_ptd,r2_ptd,F_ptd] = regress_jc(y,x,0);
fprintf('pooled forecasting regression, removing time means \n');
fprintf('    b  %8.3f %8.3f %8.3f \n',b_ptd); 

% ------------------------
% Fama MacBeth 
% ------------------------

% pure cross sectional regression

y = mean(rx)'; 
x = [ones(25,1) mean(logmc2)' mean(log(beme))']; 
[b,bse,bt,r2,F] = regress_jc(y,x,0);
fitmeans = x*b; 

fprintf('regression of mean xs return on mean log beme and mean log size\n'); 
fprintf('       %8.8s %8.8s %8.8s  \n','const','size','beme'); 
fprintf('    b  %8.3f %8.3f %8.3f  \n',b); 
fprintf('   se  %8.3f %8.3f %8.3f  \n',bse); 
fprintf('    t  %8.3f %8.3f %8.3f  \n',bt); 

% conventional calculation
e = y-x*b;
s2 = var(e);
sigb = diag(inv(x'*x)*s2).^0.5; 

err_t = rx-ones(T,1)*(x*b)'; 
sig = cov(err_t)/T;
sigb_corr = inv(x'*x)*x'*sig*x*inv(x'*x); 
sigb_corr = diag(sigb_corr).^ 0.5; 

fprintf('conv se%8.3f %8.3f %8.3f  \n',sigb); 
fprintf('    t  %8.3f %8.3f %8.3f  \n',b./sigb);

fprintf('corr se%8.3f %8.3f %8.3f  \n',sigb_corr); 
fprintf('    t  %8.3f %8.3f %8.3f  \n',b./sigb_corr);

% fama macbeth

btv = zeros(T,3);
for t = 1:T;
    y = rx(t,:)';
    x = [ones(25,1) logmc2(t,:)' log(beme(t,:))'];
    [bt,bse,btstat,r2,F] = regress_jc(y,x,0);
    btv(t,:)=bt';
end;
    
bfb = mean(btv)';
sefb = std(btv)'/T^0.5; 

fprintf(' FMB b %8.3f %8.3f %8.3f  \n',bfb); 
fprintf('    se %8.3f %8.3f %8.3f  \n',sefb); 
fprintf('    t  %8.3f %8.3f %8.3f  \n',bfb./sefb);



stop

% ---------------------
% the remaining stuff preserves my previous attempt to do both time and
% firm dummies with kronecker products
% ------------------------------------

% -----------------
% log beme alone
% ------------------

y = reshape(rx,T*25,1); 
x = [ones(T*25,1) reshape(log(beme),T*25,1)]; 
[b_b,bse_b,bt_b,r2_b,F_b] = regress_jc(y,x,0);
fit_b = x*b_b; 
r2 = var(fit_b)/var(y);

fprintf('forecasting regression of return on log beme alone\n');
fprintf('    b  %8.3f %8.3f %8.3s %8.3f  \n',b_b,' ',r2); 




% ------------------------------------------------------------
% pooled time-series cross section forecast with time dummies
% ---------------------------------------------------------

x = [ kron(ones(25,1), eye(T)) reshape(logmc2,25*T,1) reshape(log(beme),T*25,1)];  
[b_ptd,bse_ptd,bt_ptd,r2_ptd,F_ptd] = regress_jc(y,x,0);
fit_ptd = x*b_ptd; 
r2 = var(fit_ptd)/var(y);  % including dummies in r2 calculation 
Erx_ptd = reshape(fit_ptd,T,25); 

fprintf('forecasting regression of return on log beme and log size with time dummies\n');
fprintf('    b  %8.3s %8.3f %8.3f  %8.3f \n',' ',b_ptd(end-1:end),r2); 

% check with taking out time-means

logmc3 = logmc2 - mean(logmc2,2)*ones(1,25);
logbeme3 = log(beme) - mean(log(beme),2)*ones(1,25);
x = [ ones(T*25,1) reshape(logmc3,25*T,1) reshape(logbeme3,T*25,1)]; 
[b_ptd,bse_ptd,bt_ptd,r2_ptd,F_ptd] = regress_jc(y,x,0);
fprintf('    b  %8.3f %8.3f %8.3f \n',b_ptd); 

% ---------------------------------------------------------
% pooled time-series cross section forecast with portfolio dummies
% ---------------------------------------------------------

x = [kron(eye(25), ones(T,1)) reshape(logmc2,25*T,1) reshape(log(beme),T*25,1)]; 
[b_ppd,bse_ppd,bt_ppd,r2_ppd,F_ppd] = regress_jc(y,x,0);
fit_ppd = x*b_ppd; 
Erx_ppd = reshape(fit_ppd,T,25); 
r2 = var(fit_ppd)/var(y);

fprintf('forecasting regression of return on log beme and log size with portfolio dummies\n');
fprintf('    b  %8.3s %8.3f %8.3f %8.3f \n',' ',b_ppd(end-1:end),r2); 

% check taking out portfolio means
logmc3 = logmc2 - ones(T,1)*mean(logmc2);
logbeme3 = log(beme) - ones(T,1)*mean(log(beme));
x = [ones(T*25,1)  reshape(logmc3,25*T,1) reshape(logbeme3,T*25,1)]; 
[b_ppd,bpe_ppd,bt_ppd,r2_ppd,F_ppd] = regress_jc(y,x,0);
fprintf('    b  %8.3f %8.3f %8.3f \n',b_ppd); 




disp('') ;
disp('actual mean return'); 
disp(reshape(mean(rx)',5,5)'); 
disp('regression on average log size and average log beme'); 
disp(reshape(fitmeans,5,5)'); 
disp('forecasting regression using log size and log beme'); 
disp(reshape(mean(Erx_p)',5,5)'); 
disp('forecasting regression using log size and log beme'); 
disp(reshape(mean(Erx_p)',5,5)'); 


% ---------------------------------------------------------
% some graphs I did not ask for. They show that even the pooled regression
% gives a decent fit. 
% ---------------------------------------------------------

figure; 
subplot(1,2,1)
bar3(reshape(fitmeans,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %'); 
title('mean returns from mean beme,size'); 
axis([-inf inf -inf inf -0.2 2.5]);

subplot(1,2,2)
bar3(reshape(mean(Erx_p)',5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('mean mo return %'); 
title('mean returns from forecast regression'); 
axis([-inf inf -inf inf -0.2 2.5]);
print -depsc2 ps3_xs_2.eps; 

figure; 
plot(mean(Erx_p)',mean(rx)','.','MarkerSize',10);
xlabel('Fitted mean -- pooled regression')
ylabel('Actual mean') 
axis([0 1.6 0 1.6]); 
print -depsc2 fitmean3.eps