% problem set 7; follow up on Fama French and characteristics
% rewrite to use an m file for all the cases where we do things over and
% over

clear all; 
close all; 
addpath('..\data'); 

sample = 3; % 1:32-67 2:all  3:47-end; 4: 63-end

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

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

if sample == 1; 
    st = 67; en = 438+5*12; 
end;
if sample == 2; 
    st = 67; en = size(x,1)
end;
if sample == 3; 
    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);
rx = r - rf*ones(1,size(r,2));
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); 

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

disp(' '); 

x = load('25_portfolios_5x5_avb_avme.txt'); 
figure; 
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
    plot(x(6:end,1),log(x(6:end,2:end))); 
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
    plot(x(21:end,1),log(x(21:end,2:end))); 
end;
title('log sum be/sum me'); 


x = load('25_portfolios_5x5_avbeme.txt'); 
figure; 
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
    plot(x(6:end,1),log(x(6:end,2:end))); 
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
    plot(x(21:end,1),log(x(21:end,2:end))); 
end;
title('log average beme'); 
print -depsc2 log_beme_time.eps;

% based on plots, avbeme looks better. I chose not to remove trend

x = load('25_portfolios_5x5_size.txt');
if st == 67;
    mc = x(66:end-1,2:end); % market cap 1929:12 to end; using 31:12 (forecaster for 32:01) on 
    plotsizedate = floor(x(66:end-1,1)/100)+(x(66:end-1,1)-100*floor(x(66:end-1,1)/100))/12;
end;
if st == 247;
    mc = x(246:end-1,2:end); % market cap 1946:12 to end 
    plotsizedate = floor(x(246:end-1,1)/100)+(x(246:end-1,1)-100*floor(x(246:end-1,1)/100))/12;
end;

figure;
plot(plotsizedate,log(mc));
title('log size'); 
text(plotsizedate(end)*ones(25,1),log(mc(end,:))',labels);
print -depsc2 log_size_time.eps;

% obviously log size can't be compared over time. 
% Make a new variable, size = sizei/total market cap

mc2 = mc./(sum(mc,2)*ones(1,25));
logmc2 = log(mc2);
figure;
plot(plotsizedate,logmc2);
title('log size/totalmarket cap'); 
text(plotsizedate(end)*ones(25,1),logmc2(end,:)',labels);
print -depsc2 log_relative_size.eps

% expand beme to annual. Beme in year t is in june of t, thus good for
% forecasting july (07) of t
beme = kron(beme,ones(12,1)); % now june BEME is spread Jan-Dec
%beme = beme(7:end-6,:);  though it's available in july, don't start using it until jan when the portfolio is reformed
beme = beme(1:T,:);% if return does not finish out the year (e.g. stop in oct) then beme is too long (stop in dec) 

% TIMING: BEME (t) is the BEME measured in the july preceding year t.
% size(t) is the size for the month preceding t. Thus regressions of rx on
% beme and size are forecsating regressions. 


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(2+mean(log(beme)),5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('log(BE/ME)'); 
title('log(BE/ME)')

subplot(2,2,3)
bar3(reshape(mean(mc2),5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('size'); 
title('Size -- fraction of mkt')

subplot(2,2,4)
bar3(reshape(7+mean(logmc2),5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('log(size)'); 
title('log(Size) -- fraction of mkt')

print -depsc2 ps3_units.eps; 


% -------------------------------------------------------------
% 2. cross sectional regressions.  
% -------------------------------------------------------------
% puzzle: how to construct standard errors for "alphas"? FMB gives a
% different answer than cov(alphat). Dropped it from the problem set. 

% only bm

X = [ones(N,1)  mean(log(beme))'];
lambda = X \ mean(rx)';
alphat = rx-ones(T,1)*(X*lambda)';
alphas = mean(rx)'-(X*lambda);
rmsea = mean(alphas.^2)^0.5;
meanabsa = mean(abs(alphas)); 
selambda = inv(X'*X)*X'*cov(alphat)/T*X*inv(X'*X);
selambda = diag(selambda).^0.5;

fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s','rmse_a','|a|');
fprintf('cross sectional regression on -- characteristics \n'); 
fprintf('        b %8.3f %8.3f %8.3s %8.3s %8.3f %8.3f \n',lambda(1),lambda(2),' ',' ',rmsea,meanabsa); 
fprintf('        t %8.3f %8.3f                         \n', lambda(1)./selambda(1),lambda(2:end)./selambda(2:end)); 


% 2a size and bm,

X = [ones(N,1)  mean(log(beme))' mean(logmc2)'];
lambda = X \ mean(rx)';
alphat = rx-ones(T,1)*(X*lambda)';
alphas = mean(rx)'-(X*lambda);
rmsea = mean(alphas.^2)^0.5;
meanabsa = mean(abs(alphas)); 
selambda = inv(X'*X)*X'*cov(alphat)/T*X*inv(X'*X);
selambda = diag(selambda).^0.5;

fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s','rmse_a','|a|');
fprintf('cross sectional regression on -- characteristics \n'); 
fprintf('        b %8.3f %8.3f %8.3f %8.3s %8.3f %8.3f \n',lambda(1),lambda(2:end),' ',rmsea,meanabsa); 
fprintf('        t %8.3f %8.3f %8.3f                   \n', lambda(1)./selambda(1),lambda(2:end)./selambda(2:end)); 

fitmean = (X*lambda);

% 2b size and bm and interaction term,

X = [ones(N,1)  mean(log(beme))' mean(logmc2)' mean(log(beme))'.*mean(logmc2)' ];
lambda = X \ mean(rx)';
alphat = rx-ones(T,1)*(X*lambda)';
alphas = mean(rx)'-(X*lambda);
rmsea = mean(alphas.^2)^0.5;
meanabsa = mean(abs(alphas)); 
selambda = inv(X'*X)*X'*cov(alphat)/T*X*inv(X'*X);
selambda = diag(selambda).^0.5;

lambdatfb = zeros(T,size(X,2)); 
alphatfb = zeros(T,N); 
for t = 1:T;     
    lambdatfb(t,:) = (X\rx(t,:)')';
    alphatfb(t,:) = rx(t,:)-(X*lambdatfb(t,:)')' ;
end; 
lambdafb = mean(lambdatfb)'; 
selambdafb = std(lambdatfb)'/T^0.5; 
alphasfb = mean(alphatfb);
rmseafb = mean(alphasfb.^2)^0.5;
meanabsafb = mean(abs(alphasfb)); 


fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s','rmse_a','|a|');
fprintf('cross sectional regression on -- characteristics \n'); 
fprintf('        b %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n',lambda(1),lambda(2:end),rmsea,meanabsa); 
fprintf('        t %8.3f %8.3f %8.3f %8.3f                  \n', lambda(1)./selambda(1),lambda(2:end)./selambda(2:end)); 
fprintf('same, by FMB with same right hand variable at each date \n'); 
fprintf('        b %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n',lambdafb(1),lambdafb(2:end),rmseafb,meanabsafb); 
fprintf('        t %8.3f %8.3f %8.3f %8.3f                  \n', lambdafb(1)./selambdafb(1),lambdafb(2:end)./selambdafb(2:end)); 

fitmean_cross = X*lambda; 


% 2b interact variables not mean?,

X = [ones(N,1)  mean(log(beme))' mean(logmc2)' mean(log(beme).*logmc2)' ];
lambda = X \ mean(rx)';
alphat = rx-ones(T,1)*(X*lambda)';
alphas = mean(rx)'-(X*lambda);
rmsea = mean(alphas.^2)^0.5;
meanabsa = mean(abs(alphas)); 
selambda = inv(X'*X)*X'*cov(alphat)/T*X*inv(X'*X);
selambda = diag(selambda).^0.5;

fprintf('                     cross variables not means \n'); 
fprintf('        b %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n',lambda(1),lambda(2:end),rmsea,meanabsa); 
fprintf('        t %8.3f %8.3f %8.3f % 8.3f                  \n', lambda(1)./selambda(1),lambda(2:end)./selambda(2:end)); 

% 2c less than three ? 

X = [ones(N,1) mean(log(beme))' mean(log(beme))'.*mean(logmc2)' ];
lambda = X \ mean(rx)';
alphat = rx-ones(T,1)*(X*lambda)';
alphas = mean(rx)'-(X*lambda);
rmsea = mean(alphas.^2)^0.5;
meanabsa = mean(abs(alphas)); 
selambda = inv(X'*X)*X'*cov(alphat)/T*X*inv(X'*X);
selambda = diag(selambda).^0.5;

fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s','rmse_a','|a|');
fprintf('cross sectional regression on -- characteristics \n'); 
fprintf('        b %8.3f %8.3f %8.3s %8.3f %8.3f %8.3f \n',lambda(1),lambda(2), ' ',lambda(3),rmsea,meanabsa); 
fprintf('        t %8.3f %8.3f %8.3s %8.3f                  \n', lambda(1)./selambda(1),lambda(2)./selambda(2), ' ',lambda(3)./selambda(3)); 

X = [ones(N,1) mean(logmc2)' mean(log(beme))'.*mean(logmc2)' ];
lambda = X \ mean(rx)';
alphat = rx-ones(T,1)*(X*lambda)';
alphas = mean(rx)'-(X*lambda);
rmsea = mean(alphas.^2)^0.5;
meanabsa = mean(abs(alphas));
selambda = inv(X'*X)*X'*cov(alphat)/T*X*inv(X'*X);
selambda = diag(selambda).^0.5;
fitmean_b_x = X*lambda;

fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s','rmse_a','|a|');
fprintf('cross sectional regression on -- characteristics \n');
fprintf('        b %8.3f %8.3s %8.3f %8.3f %8.3f %8.3f \n',lambda(1),' ', lambda(2),lambda(3),rmsea,meanabsa);
fprintf('        t %8.3f %8.3s %8.3f %8.3f                  \n', lambda(1)./selambda(1),' ', lambda(2)./selambda(2),lambda(3)./selambda(3));

X = [ones(N,1) mean(log(beme))'.*mean(logmc2)' ];
lambda = X \ mean(rx)';
alphat = rx-ones(T,1)*(X*lambda)';
alphas = mean(rx)'-(X*lambda);
rmsea = mean(alphas.^2)^0.5;
meanabsa = mean(abs(alphas));
selambda = inv(X'*X)*X'*cov(alphat)/T*X*inv(X'*X);
selambda = diag(selambda).^0.5;

fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s','rmse_a','|a|');
fprintf('cross sectional regression on -- characteristics \n');
fprintf('        b %8.3f %8.3s %8.3s %8.3f %8.3f %8.3f \n',lambda(1),' ', ' ',lambda(2),rmsea,meanabsa);
fprintf('        t %8.3f %8.3s %8.3s %8.3f                  \n', lambda(1)./selambda(1),' ', ' ',lambda(2)./selambda(2));


% -------------------------
% Pictures
% ---------------------------

figure;
subplot(2,2,1);
bar3(reshape(mean(rx),5,5)');
ylabel('size');
xlabel('value');
zlabel('E(Re)');
title('E(Re)')

subplot(2,2,2);
bar3(reshape(fitmean,5,5)');
ylabel('size');
xlabel('value');
zlabel('E(Re)');
title('Fit mean, no cross term')

subplot(2,2,3);
bar3(reshape(fitmean_cross,5,5)');
ylabel('size');
xlabel('value');
zlabel('E(Re)');
title('Fit mean, with cross term')

subplot(2,2,4);
bar3(reshape(mean(rx)'-fitmean_cross,5,5)');
ylabel('size');
xlabel('value');
zlabel('E(Re)-fit');
title('errors, with cross term')

print -depsc2 bar_fit_means.eps;


if 0
    figure ;
    subplot(1,2,1);
    plot(reshape(fitmean,5,5)',reshape(mean(rx),5,5)','.-','linewidth',2);
    hold on
    errorbar(reshape(fitmean,5,5)',reshape(fitmean,5,5)',reshape(std(rx)/size(rx,1)^0.5,5,5)');
    plot([0 1.6],[0 1.6]);
    axis([0 1.6 0 1.6]);
    xlabel('Fitted mean return');
    ylabel('mean return') ;
    title('ERe = s + bm cross sectional');
    % lines connect size groups within value
    
    subplot(1,2,2)
    plot(reshape(fitmean_cross,5,5)',reshape(mean(rx),5,5)','.-','linewidth',2);
    hold on
    plot(reshape(fitmean_b_x,5,5)',reshape(mean(rx),5,5)','v-');
    errorbar(reshape(fitmean_b_x,5,5)',reshape(fitmean_b_x,5,5)',reshape(std(rx)/size(rx,1)^0.5,5,5)');
    plot([0 1.6],[0 1.6]);
    axis([0 1.6 0 1.6]);
    xlabel('Fitted mean return');
    ylabel('mean return') ;
    title('ERe = s (+ bm) + s x bm cross sectional');
end

figure ; 
subplot(1,2,1);
plot(reshape(mean(rx),5,5)',reshape(fitmean,5,5)','.-','linewidth',2);
hold on
errorbar(reshape(mean(rx),5,5)',reshape(mean(rx),5,5)',reshape(std(rx)/size(rx,1)^0.5,5,5)');
plot([0 1.6],[0 1.6]);
axis([0 1.6 0 1.6]);
xlabel('Mean return'); 
ylabel('Fitted mean return') ; 
title('ERe = s + bm cross sectional'); 
% lines connect size groups within value

subplot(1,2,2)
plot(reshape(mean(rx),5,5)',reshape(fitmean_cross,5,5)','.-','linewidth',2);
hold on
plot(reshape(mean(rx),5,5)',reshape(fitmean_b_x,5,5)','v-');
errorbar(reshape(mean(rx),5,5)',reshape(mean(rx),5,5)',reshape(std(rx)/size(rx,1)^0.5,5,5)');
plot([0 1.6],[0 1.6]);
axis([0 1.6 0 1.6]);
xlabel('Mean return'); 
ylabel('Fitted mean return') ; 
title('ERe = s (+ bm) + s x bm cross sectional'); 
set(gcf,'paperposition',[0.25 2.5 8 3])

print -depsc2 xy_fit_means.eps;


% ---------------------------
% betas? 
% ---------------------------


betas = [ones(T,1) rmrf hml smb]\rx; 
X = [ones(N,1)  mean(log(beme))' mean(logmc2)'  ];
cd = X \ betas';
lambda = X \ mean(rx)';

disp(' '); 
disp(' betas')
disp('using size and beme'); 
fprintf('          %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s');
fprintf('cross sectional regression of betas on characteristics \n'); 
fprintf('        a %8.3f %8.3f %8.3f   \n',cd(:,1));
fprintf('        b %8.3f %8.3f %8.3f   \n',cd(:,2)); 
fprintf('        h %8.3f %8.3f %8.3f   \n',cd(:,3)); 
fprintf('        s %8.3f %8.3f %8.3f   \n',cd(:,4)); 
fprintf('cross sectional regression of means on characteristics \n'); 
fprintf('       Er %8.3f %8.3f %8.3f   \n',lambda); 
fprintf('sum of betas weighted by factor risk premia \n'); 
fit = cd(:,2)*mean(rmrf)+cd(:,3)*mean(hml)+cd(:,4)*mean(smb);
fprintf('      b*l %8.3f %8.3f %8.3f  \n',fit); 
fprintf(' Er-  b*l %8.3f %8.3f %8.3f  \n',lambda - fit); 
fprintf('        a %8.3f %8.3f %8.3f   \n',cd(:,1));


figure;
subplot(2,2,1);
bar3(reshape(X*lambda,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('E(Re)'); 
title('Fit E(Re)')

subplot(2,2,2);
bar3(reshape(X*fit,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('E(Re)'); 
title('Fit beta*lambda')

subplot(2,2,3);
bar3(reshape(X*lambda-X*fit,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('alpha'); 
title('Fit E(Re)-beta*lambda')

subplot(2,2,4);
bar3(reshape(betas(1,:)',5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('alpha'); 
title('individual alphas')
print -depsc2 ff_table_fit_size_beme.eps

% ----------

disp(' '); 
disp('using size beme and cross term'); 
betas = [ones(T,1) rmrf hml smb]\rx; 
X = [ones(N,1)  mean(log(beme))' mean(logmc2)' mean(log(beme))'.*mean(logmc2)' ];
cd = X \ betas';
lambda = X \ mean(rx)';

fprintf('          %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s');
fprintf('cross sectional regression of betas on characteristics \n'); 
fprintf('        a %8.3f %8.3f %8.3f %8.3f  \n',cd(:,1));
fprintf('        b %8.3f %8.3f %8.3f %8.3f  \n',cd(:,2)); 
fprintf('        h %8.3f %8.3f %8.3f %8.3f  \n',cd(:,3)); 
fprintf('        s %8.3f %8.3f %8.3f %8.3f  \n',cd(:,4)); 
fprintf('cross sectional regression of means on characteristics \n'); 
fprintf('       Er %8.3f %8.3f %8.3f %8.3f  \n',lambda); 
fprintf('sum of betas weighted by factor risk premia \n');
fit = cd(:,2)*mean(rmrf)+cd(:,3)*mean(hml)+cd(:,4)*mean(smb);
fprintf('      b*l %8.3f %8.3f %8.3f %8.3f  \n',fit); 
fprintf('  Er- b*l %8.3f %8.3f %8.3f %8.3f  \n',lambda - fit); 
fprintf('        a %8.3f %8.3f %8.3f %8.3f  \n',cd(:,1));

% ---------------------

disp(' '); 
disp('using size and cross term'); 

betas = [ones(T,1) rmrf hml smb]\rx; 
X = [ones(N,1) mean(logmc2)' (mean(log(beme)).*mean(logmc2))'  ];
cd = X \ betas';
lambda = X \ mean(rx)';

fprintf('          %8.8s %8.8s %8.8s %8.8s \n', 'const','beme','size','bm x s');
fprintf('cross sectional regression of betas on characteristics \n'); 
fprintf('        a %8.3f %8.8s %8.3f %8.3f   \n',cd(1,1), ' ', cd(2:end,1)); 
fprintf('        b %8.3f %8.8s %8.3f %8.3f   \n',cd(1,2), ' ', cd(2:end,2)); 
fprintf('        h %8.3f %8.8s %8.3f %8.3f   \n',cd(1,3), ' ', cd(2:end,3)); 
fprintf('        s %8.3f %8.8s %8.3f %8.3f   \n',cd(1,4), ' ', cd(2:end,4)); 
fprintf('cross sectional regression of means on characteristics \n'); 
fprintf('       Er %8.3f %8.8s %8.3f %8.3f   \n',lambda(1), ' ', lambda(2:end)); 
fprintf('sum of betas weighted by factor risk premia \n'); 
fit = cd(:,2)*mean(rmrf)+cd(:,3)*mean(hml)+cd(:,4)*mean(smb);
fprintf('      B*l %8.3f %8.8s %8.3f %8.3f  \n',fit(1), ' ', fit(2:end)); 
fprintf('   Er-B*l %8.3f %8.8s %8.3f %8.3f  \n',lambda(1)-fit(1), ' ', lambda(2:end)-fit(2:end)); 
fprintf('        a %8.3f %8.8s %8.3f %8.3f   \n',cd(1,1), ' ', cd(2:end,1));


figure;
subplot(2,2,1);
bar3(reshape(X*lambda,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('E(Re)'); 
title('Fit E(Re)')

subplot(2,2,2);
bar3(reshape(X*fit,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('E(Re)'); 
title('Fit beta*lambda')

subplot(2,2,3);
bar3(reshape(X*lambda-X*fit,5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('alpha'); 
title('Fit E(Re)-beta*lambda')

subplot(2,2,4);
bar3(reshape(betas(1,:)',5,5)');
ylabel('size'); 
xlabel('value'); 
zlabel('alpha'); 
title('individual alphas')
print -depsc2 ff_table_fit.eps


% -------------------------------
% connection to forecasting regressions
% ---------------------------------

% ----------
% just bm
% -----------

% reminder: cross section 
X = [ones(N,1)  mean(log(beme))'  ];
lambda = X \ mean(rx)';


% pooled without cross term 
% for standard errors orgnaize by N not by T, i.e. 1, 2, 3, ...25 at time
% 1, then 1,2,3...25 at time 2 etc. Then covmat = [sigma 00; 0 sigma 00..];

y = reshape(rx',T*25,1); 
x = [ones(T*25,1)  reshape(log(beme)',T*25,1)  ]; 
b = x\y; 
fit_pooled = x*b; 
x_pooled = x; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;
disp(' ')
disp('FORECASTS -- using BM only')
fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s  \n', 'const','beme','size','bm x s','R2');
fprintf('cross-sec %8.3f %8.3f  \n',lambda); 

fprintf('pooled forecast estimates  \n'); 
fprintf('  coeffs   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b,' ',' ',100*R2); 
fprintf('       t   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b./seb,' ',' ', 100*R2); 



%time dummies
y = reshape(rx',T*25,1); 
beme_t = log(beme)-mean(log(beme),2)*ones(1,25); % take out average at each time
x = [ones(T*25,1) reshape(beme_t',T*25,1) ];
b = x\y; 
fit_t = x*b; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;

fprintf('time dummies  \n'); 
fprintf('  coeffs   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b,' ',' ',100*R2); 
fprintf('       t   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b./seb,' ',' ', 100*R2); 

% portfolio dummies 
y = reshape(rx',T*25,1); 
beme_p = log(beme)-ones(T,1)*mean(log(beme),1); % take out average for each portfolio
x = [ones(T*25,1) reshape(beme_p',T*25,1)];
b = x\y; 
fit_p = x*b; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;
fprintf('portfolio dummies  \n'); 
fprintf('  coeffs   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b,' ',' ',100*R2); 
fprintf('       t   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b./seb,' ',' ', 100*R2); 

% check on individual portfolios to see if coefficients and std errors are
% reasonable

fprintf('check -- individual portfolios (means portfolio dummy)  \n'); 
% check -- individual regressions
for i = 1:25; 
    [b(:,i),seb(:,i),R2(i,1),R2vadj,v,F] = olsgmm(rx(:,i),[ones(T,1) log(beme(:,i)) ],0,0);
end;
fprintf(' avg coeff %8.3f %8.3f %8.3s %8.3s %8.3f  \n',mean(b'),' ',' ',100*mean(R2));
fprintf(' avg t     %8.3f %8.3f                    \n',mean(b'./seb'));


figure; 
plot(x_pooled(:,2),fit_pooled,'-');
xlabel('log b/m'); 
ylabel('E_t(R_{t+1})'); 
axis([-2.2 2 -0.5 2.2]); 
hold on
plot(reshape(x_pooled(:,2),N,T),reshape(fit_t,N,T),'-g'); 
plot(reshape(x_pooled(:,2),N,T)',reshape(fit_p,N,T)','-r');
legend('pooled','time','portfolio',4);
print -depsc2 er_time_and_portfolio_dummy.eps


% ----------
% just bm, dynamics (attempt to get rid of time vs. portfolio dummy effects
% -----------
if 0 % needs work. Get this to work! idea: dynamics should give a relation that is stable with time, portfolio dummies
    
    % reminder: cross section
    X = [ones(N,1)  mean(log(beme))'  mean(log(beme(13:T,:)-log(beme(1:T-12,:)))) ];
    lambda = X \ mean(rx)';
    
    
    % pooled without cross term
    % for standard errors orgnaize by N not by T, i.e. 1, 2, 3, ...25 at time
    % 1, then 1,2,3...25 at time 2 etc. Then covmat = [sigma 00; 0 sigma 00..];
    
    y = reshape(rx(13:T,:)',(T-12)*25,1);
    x = [ones((T-12)*25,1)  reshape(log(beme(13:T,:)',(T-12)*25,1)) reshape((log(beme(13:T,:))-log(beme(1:T-12,:)))',(T-12)*25,1)  ];
        b = x\y;
        fit_pooled = x*b;
        x_pooled = x;
        e = y-x*b;
        R2 = 1-var(e)/var(y);
        err = reshape(e,N,T)';
        sigma = cov(err);
        xx = inv(x'*x);
        xcx = x(1:N,:)'*sigma*x(1:N,:);
        for i = 2:T-12;
        xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
        end;
        covmat = xx*xcx*xx;
        seb = diag(covmat).^0.5;
        disp(' ')
        disp('FORECASTS -- using BM dyanamics ')
        fprintf('          %8.8s %8.8s %8.8s  \n', 'const','beme','DBeme');
        fprintf('cross-sec %8.3f %8.3f %8.3f  \n',lambda);
        
        fprintf('pooled forecast estimates  \n');
        fprintf('  coeffs   %8.3f %8.3f %8.3f %8.3s %8.3f  \n',b,' ',' ',100*R2);
        fprintf('       t   %8.3f %8.3f %8.3f %8.3s %8.3f  \n',b./seb,' ',' ', 100*R2);
        
        
        
        %time dummies
        y = reshape(rx',T*25,1);
        beme_t = log(beme)-mean(log(beme),2)*ones(1,25); % take out average at each time
        x = [ones(T*25,1) reshape(beme_t',T*25,1) ];
        b = x\y;
        fit_t = x*b;
        e = y-x*b;
        R2 = 1-var(e)/var(y);
        err = reshape(e,N,T)';
        sigma = cov(err);
        xx = inv(x'*x);
        xcx = x(1:N,:)'*sigma*x(1:N,:);
        for i = 2:T;
            xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
        end;
        covmat = xx*xcx*xx;
        seb = diag(covmat).^0.5;
        
        fprintf('time dummies  \n');
        fprintf('  coeffs   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b,' ',' ',100*R2);
        fprintf('       t   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b./seb,' ',' ', 100*R2);
        
        % portfolio dummies
        y = reshape(rx',T*25,1);
        beme_p = log(beme)-ones(T,1)*mean(log(beme),1); % take out average for each portfolio
        x = [ones(T*25,1) reshape(beme_p',T*25,1)];
        b = x\y;
        fit_p = x*b;
        e = y-x*b;
        R2 = 1-var(e)/var(y);
        err = reshape(e,N,T)';
        sigma = cov(err);
        xx = inv(x'*x);
        xcx = x(1:N,:)'*sigma*x(1:N,:);
        for i = 2:T;
            xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
        end;
        covmat = xx*xcx*xx;
        seb = diag(covmat).^0.5;
        fprintf('portfolio dummies  \n');
        fprintf('  coeffs   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b,' ',' ',100*R2);
        fprintf('       t   %8.3f %8.3f %8.3s %8.3s %8.3f  \n',b./seb,' ',' ', 100*R2);
        
        % check on individual portfolios to see if coefficients and std errors are
        % reasonable
        
        fprintf('check -- individual portfolios (means portfolio dummy)  \n');
        % check -- individual regressions
        for i = 1:25;
            [b(:,i),seb(:,i),R2(i,1),R2vadj,v,F] = olsgmm(rx(:,i),[ones(T,1) log(beme(:,i)) ],0,0);
        end;
        fprintf(' avg coeff %8.3f %8.3f %8.3s %8.3s %8.3f  \n',mean(b'),' ',' ',100*mean(R2));
        fprintf(' avg t     %8.3f %8.3f                    \n',mean(b'./seb'));
        
        
        figure;
        plot(x_pooled(:,2),fit_pooled,'-');
        xlabel('log b/m');
        ylabel('E_t(R_{t+1})');
        axis([-2.2 2 -0.5 2.2]);
        hold on
        plot(reshape(x_pooled(:,2),N,T),reshape(fit_t,N,T),'-g');
        plot(reshape(x_pooled(:,2),N,T)',reshape(fit_p,N,T)','-r');
        legend('pooled','time','portfolio',4);
        print -depsc2 er_diff_time_and_portfolio_dummy.eps
end



% --------------------------
% bm and size
% ------------------------

X = [ones(N,1)  mean(log(beme))' mean(logmc2)'  ];
lambda = X \ mean(rx)';

% pooled without cross term 
% for standard errors orgnaize by N not by T, i.e. 1, 2, 3, ...25 at time
% 1, then 1,2,3...25 at time 2 etc. Then covmat = [sigma 00; 0 sigma 00..];

y = reshape(rx',T*25,1); 
x = [ones(T*25,1)  reshape(log(beme)',T*25,1)  reshape(logmc2',25*T,1)  ]; 
b = x\y; 
fit_pooled = x*b; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;
disp(' ')
disp('FORECASTS -- using BM and SIZE')
fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s  \n', 'const','beme','size','bm x s','R2');
fprintf('cross-sec %8.3f %8.3f %8.3f  \n',lambda); 
fprintf('pooled forecast estimates  \n'); 
fprintf('  coeffs   %8.3f %8.3f %8.3f %8.3s %8.3f  \n',b,' ',100*R2); 
fprintf('       t   %8.3f %8.3f %8.3f %8.3s %8.3f  \n',b./seb,' ',100*R2); 



%time dummies
y = reshape(rx',T*25,1); 
beme_t = log(beme)-mean(log(beme),2)*ones(1,25); % take out average at each time
mc_t= logmc2-mean(logmc2,2)*ones(1,25);
x = [ones(T*25,1) reshape(beme_t',T*25,1) reshape(mc_t',T*25,1)];
b = x\y; 
fit_t = x*b; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;

fprintf('time dummies  \n'); 
fprintf('  coeffs   %8.3f %8.3f %8.3f %8.3s %8.3f  \n',b,' ',100*R2); 
fprintf('       t   %8.3f %8.3f %8.3f %8.3s %8.3f  \n',b./seb,' ',100*R2); 

% portfolio dummies 
y = reshape(rx',T*25,1); 
beme_p = log(beme)-ones(T,1)*mean(log(beme),1); % take out average for each portfolio
mc_p= logmc2-ones(T,1)*mean(logmc2,1);
x = [ones(T*25,1) reshape(beme_p',T*25,1) reshape(mc_p',T*25,1)];
b = x\y; 
fit_p = x*b; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;
fprintf('portfolio dummies  \n'); 
fprintf('  coeffs   %8.3f %8.3f %8.3f %8.3s %8.3f  \n',b,' ',100*R2); 
fprintf('       t   %8.3f %8.3f %8.3f %8.3s %8.3f  \n',b./seb,' ',100*R2); 

% check on individual portfolios to see if coefficients and std errors are
% reasonable

fprintf('check -- individual portfolios  \n'); 
% check -- individual regressions
for i = 1:25; 
    [b(:,i),seb(:,i),R2(i,1),R2vadj,v,F] = olsgmm(rx(:,i),[ones(T,1) log(beme(:,i)) logmc2(:,i)],0,0);
end;
fprintf(' avg coeff %8.3f %8.3f %8.3f %8.3s %8.3f  \n',mean(b'),' ',100*mean(R2));
fprintf(' avg t     %8.3f %8.3f %8.3f        \n',mean(b'./seb'));


% --------------------------
%  size and  cross term 
% ------------------------
X = [ones(N,1) mean(logmc2)' mean(log(beme))'.*mean(logmc2)' ];
lambda = X \ mean(rx)';

% pooled 

y = reshape(rx',T*25,1); 
x = [ones(T*25,1)  reshape(logmc2',25*T,1) reshape((log(beme).*logmc2)',T*25,1) ]; 
b = x\y; 
fit_pooled = x*b; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;
disp(' ')
disp('FORECASTS -- using size, BMxsize ')
fprintf('           %8.8s %8.8s %8.8s %8.8s %8.8s  \n', 'const','beme','size','bm x s','R2');
fprintf('cross-sec  %8.3f %8.3s %8.3f %8.3f        \n',lambda(1),' ', lambda(2:end)); 
fprintf('pooled forecast estimates  \n'); 
fprintf('  coeffs   %8.3f %8.3s %8.3f %8.3f %8.3f  \n',b(1),' ', b(2:end),100*R2); 
fprintf('       t   %8.3f %8.3s %8.3f %8.3f %8.3f  \n',b(1)./seb(1),' ',b./seb,100*R2); 


%time dummies
y = reshape(rx',T*25,1); 
beme_t = log(beme)-mean(log(beme),2)*ones(1,25); % take out average at each time
mc_t= logmc2-mean(logmc2,2)*ones(1,25);
cross_t = (log(beme).*logmc2)-mean(log(beme).*logmc2,2)*ones(1,N);
x = [ones(T*25,1) reshape(mc_t',T*25,1) reshape(cross_t',T*25,1)];
b = x\y; 
fit_t = x*b; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;

fprintf('time dummies  \n'); 
fprintf('  coeffs   %8.3f %8.3s %8.3f %8.3f %8.3f  \n',b(1),' ', b(2:end),100*R2); 
fprintf('       t   %8.3f %8.3s %8.3f %8.3f %8.3f  \n',b(1)./seb(1),' ',b./seb,100*R2); 


% portfolio dummies 
y = reshape(rx',T*25,1); 
beme_p = log(beme)-ones(T,1)*mean(log(beme),1); % take out average for each portfolio
mc_p= logmc2-ones(T,1)*mean(logmc2,1);
cross_p = (log(beme).*logmc2)-ones(T,1)*mean(log(beme).*logmc2,1);
x = [ones(T*25,1) reshape(mc_p',T*25,1) reshape(cross_p',T*25,1)];
b = x\y; 
fit_p = x*b; 
e = y-x*b;
R2 = 1-var(e)/var(y);
err = reshape(e,N,T)';
sigma = cov(err); 
xx = inv(x'*x);
xcx = x(1:N,:)'*sigma*x(1:N,:); 
for i = 2:T; 
    xcx = xcx+x(N*(i-1)+1:N*i,:)'*sigma*x(N*(i-1)+1:N*i,:);
end; 
covmat = xx*xcx*xx; 
seb = diag(covmat).^0.5;
fprintf('portfolio dummies  \n'); 
fprintf('  coeffs   %8.3f %8.3s %8.3f %8.3f %8.3f  \n',b(1),' ', b(2:end),100*R2); 
fprintf('       t   %8.3f %8.3s %8.3f %8.3f %8.3f  \n',b(1)./seb(1),' ',b./seb,100*R2); 


% check on individual portfolios to see if coefficients and std errors are
% reasonable

clear b
fprintf('check -- individual portfolios  \n'); 
% check -- individual regressions
for i = 1:25; 
    [b(:,i),seb(:,i),R2(i,1),R2vadj,v,F] = olsgmm(rx(:,i),[ones(T,1) logmc2(:,i) log(beme(:,i)).*logmc2(:,i) ],0,0);
end;
if 0; % allows printout to diagnose explosion of three variable model 
    disp('coefficients')
    for i = 1:25; 
        fprintf('  %8.8s %8.3f %8.3s %8.3f %8.3f %8.3f  \n',labels{i},b(1,i),' ',b(2:end,i),100*R2(i,1));
    end
    disp('t'); 
    for i = 1:25; 
        fprintf('  %8.8s %8.3f %8.3s %8.3f %8.3f        \n',labels{i},b(1,1)./seb(1,i),' ',b(2:end,1)./seb(2:end,i))
    end
end;
meanb = mean(b');
tb = mean(b'./seb');
fprintf(' avg coeff %8.3f %8.3s %8.3f %8.3f %8.3f  \n',meanb(1),' ',meanb(2:end),100*mean(R2));
fprintf(' avg t     %8.3f %8.3s %8.3f %8.3f        \n',tb(1),' ',tb(2:end));


% pooled time series and cross sectional regression with cross term

y = reshape(rx,T*25,1); 
x = [ones(T*25,1)  reshape(log(beme),T*25,1)  reshape(logmc2,25*T,1)  reshape(logmc2.*log(beme),25*T,1)]; 
lambda_p = x\y; 
fit_p = x*lambda_p; % pooled return forecast 
fcst_err = y-fit_p; 
Erx_p = reshape(fit_p,T,25); 
alphas_p = mean(reshape(fcst_err,T,25))';
rmsea_p = mean(alphas_p.^2)^0.5;
meanabsa_p = mean(abs(alphas_p)); 

%fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'const','b','h','s','beme','size','bm x s','rmse_a','|a|');
fprintf('pooled forecast estimates  \n'); 
fprintf(' lambda   %8.3f %8.3s %8.3s %8.3s %8.3f %8.3f %8.3f %8.3f %8.3f \n',lambda_p(1),' ',' ',' ',lambda_p(2:end),rmsea_p,meanabsa_p); 
%fprintf('        t %8.3f %8.3s %8.3s %8.3s %8.3f %8.3f %8.3f  \n', lambda_p(1)./selambda_p(1),' ',' ',' ',lambda_p(2:end)./selambda_p(2:end)); 

% pooled time series and cross sectional regression, time dummies

y = reshape(rx,T*25,1); 
x = [ones(T*25,1)  reshape(log(beme),T*25,1)  reshape(logmc2,25*T,1)  reshape(logmc2.*log(beme),25*T,1)]; 
lambda_p = x\y; 
fit_p = x*lambda_p; % pooled return forecast 
fcst_err = y-fit_p; 
Erx_p = reshape(fit_p,T,25); 
alphas_p = mean(reshape(fcst_err,T,25))';
rmsea_p = mean(alphas_p.^2)^0.5;
meanabsa_p = mean(abs(alphas_p)); 

%fprintf('          %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s \n', 'const','b','h','s','beme','size','bm x s','rmse_a','|a|');
fprintf('pooled forecast estimates  \n'); 
fprintf(' lambda   %8.3f %8.3s %8.3s %8.3s %8.3f %8.3f %8.3f %8.3f %8.3f \n',lambda_p(1),' ',' ',' ',lambda_p(2:end),rmsea_p,meanabsa_p); 
%fprintf('        t %8.3f %8.3s %8.3s %8.3s %8.3f %8.3f %8.3f  \n', lambda_p(1)./selambda_p(1),' ',' ',' ',lambda_p(2:end)./selambda_p(2:end)); 

                           

