% problem set 7 
% reads in mutual fund data and processes it. 

close all; 
clear all; 
addpath('..\Data'); 
addpath('..\Programs'); 

% -------------------
% load data , make carhart portfolios
% ------------------

if 0; % If 1 is how I create the data, then save it for students. Commented out so it runs on student data
    x = load('fund_data_2011');
    x = x.output;
    date = x(1:end-9,1);
    x = x(1:end-9,2:end).*100; % program wants data that end in december 2011 update ends in september
    fprintf('using sample %8.0i to %8.0i \n',date(1), date(end))
    y = load('F-F_Research_Data_Factors_monthly.txt');
    I = find(and(y(:,1)>=date(1),y(:,1)<=date(end)));
    datesff = y(I,1);
    rmrf = y(I,2);
    smb = y(I,3);
    hml = y(I,4);
    rf = y(I,5);
    y = load('F-F_Momentum_Factor_monthly.txt');
    I = find(and(y(:,1)>=date(1),y(:,1)<=date(end)));
    umd = y(I,2);
    
    save fund_factor_data_2011 date x rmrf smb hml umd rf
    dlmwrite('fund_data_2011.csv',[date rmrf smb hml umd rf x],'precision',8)
    
    clear all % check to make sure it works with student file
end;

load  fund_factor_data_2011 date x rmrf smb hml umd rf
datesff = date; % remaining program used both, they're identical 

x = x - rf*ones(1,size(x,2));

% ------------------------------------
% form carhart portfolios. 
% programming note: I use two kinds of logical subscripting here to pick which returns to look at. 
% For example, we have to ignore the NaN (not available) returns. If you have a matrix x = [ 3 5 8 2 NaN]' 
% and you only want the 2,4 elements, there are two ways to do it. First, you can say y = x([2 4]). This generates y = [5 2]'; 
% Second you can make a matrix conformable with x of ones and zeros, sel = [0 1 0 1 0]. Then y = x(sel) produces the same thing. 
% The latter is particularly useful for the results of logical operations. For example isnan(x) produces a matrix with 1 in the position 
% of each NaN in x and a zero in the place of each number; ~isnan(x) (not isnan) produces a matrix with a 1 in the position of each number 
% and a zero in the position of each NaN. Thus, y = x(~isnan(x)) produces [ 3 5 8 2 ]';

'making carhart portfolios'

ta = 0; 
carhart_port = NaN*zeros(size(x,1),10); % set up space for carhart portfolio returns 
biglag = [];                            % set this up so you can use biglag = [biglag; new ]
for pickdate = 12:12:(size(x,1)-12);    % form portfolios on dec 31 of each date starting with end of year 2. 
    ta = ta+1; 
    alive_1yr = sum(~isnan(x(pickdate-11:pickdate,:)))==12; % 1 if fund is alive at pickdate, and has 1 year of return history
    ret1yr = mean(x(pickdate-11:pickdate,alive_1yr));       % average monthly return for last year, for all the alive1yr funds
        
    % save variables for fund-level regression -- this idea isn't working,
    % produces way too big a coefficient. Why? 
    if 0
        lag_and_next_month = [ret1yr;x(pickdate+1,alive_1yr)];
        lag_and_next_month = lag_and_next_month(:,~any(isnan(lag_and_next_month))); %eliminate funds with no next return
        biglag = [biglag lag_and_next_month]; % build huge matrix for fund-level regression
    end
    
    lag_and_next = [ret1yr;x(pickdate+1:pickdate+12,alive_1yr)];  % for all alive1yr funds, keep next year with last year
    lag_and_next = sortrows(lag_and_next',1)';                    % sort horizontally by lagged average return
    breakpoints = floor((0.1:0.1:1)'*size(lag_and_next,2)); 
    breakpoints = [0;breakpoints];                                 % indices of funds at decile breakpoints
    for i = 1:10;                                 % now form decile portfolio returns for the next year. loop through portfolios 
        retkeep(ta,i) = mean(lag_and_next(1,breakpoints(i)+1:breakpoints(i+1))'); % last year return, saved for portfolio means, regression
        carhart_port(pickdate+1:pickdate+12,i)=nanmean(lag_and_next(2:end,breakpoints(i)+1:breakpoints(i+1)),2); 
        frac_beat_mkt(ta,i) = sum((mean(lag_and_next(2:13,breakpoints(i)+1:breakpoints(i+1)))-mean(rmrf(pickdate+1:pickdate+12)))' >0 );
        frac_less_mkt(ta,i) = sum((mean(lag_and_next(2:13,breakpoints(i)+1:breakpoints(i+1)))-mean(rmrf(pickdate+1:pickdate+12)))' <0 );
        total = frac_beat_mkt(ta,i)+frac_less_mkt(ta,i);
        frac_beat_mkt(ta,i)=frac_beat_mkt(ta,i)./total;
        frac_less_mkt(ta,i)=frac_less_mkt(ta,i)./total;
    end;
end; 

% chop to available dates 

carhart_port = carhart_port(13:end,:); 
datesff_1 = datesff(13:end,:);
rmrf_1 = rmrf(13:end,:); 
hml_1 = hml(13:end,:); 
smb_1 = smb(13:end,:); 
umd_1 = umd(13:end,:);
rf_1 = rf(13:end,:); 


% --------------------------------------------------------
% form carhart-style portfolios using 5 years of returns. 
% --------------------------------------------------------
'making carhart portfolios with 5 years of returns'

ta = 0; 
for pickdate = 60:12:(size(x,1)-12); % form portfolios on dec 31 of each date starting with end of year 6.
    ta = ta+1; 
    alive_5yr = sum(~isnan(x(pickdate-59:pickdate,:)))==60; % 1 if fund is alive at pickdate, and has 5 years of return history
    ret5yr = mean(x(pickdate-59:pickdate,alive_5yr)); % return from t-5 to t-1 
    lag_and_next = [ret5yr;x(pickdate+1:pickdate+12,alive_5yr)];  % keep next year with last 5 year average
    lag_and_next = sortrows(lag_and_next',1)'; % sort horizontally by lagged average return
    breakpoints = floor((0.1:0.1:1)'*size(lag_and_next,2)); 
    breakpoints = [0;breakpoints]; % indices of funds at decile breakpoints
    for i = 1:10; % now form decile portfolio returns for the next year 
        ret5yrkeep(ta,i) = mean(lag_and_next(1,breakpoints(i)+1:breakpoints(i+1))');
        for t = 1:12; 
            %sels = find(~isnan(lag_and_next(t+1,breakpoints(i)+1:breakpoints(i+1))')); % find non NAN returns in this month 
            %num_in_port_5a(pickdate+t,i) = length(sels); 
            %carhart_port_5a(pickdate+t,i) = nanmean(lag_and_next(t+1,breakpoints(i)+sels-1),2);
            carhart_port_5a(pickdate+t,i) = nanmean(lag_and_next(t+1,breakpoints(i)+1:breakpoints(i+1)),2);
        end; 
    end; 
end; 

% chop to available dates 

datesff_5a = datesff(61:end); 
carhart_port_5a = carhart_port_5a(61:end,:); 
rmrf_5a = rmrf(61:end); 
hml_5a = hml(61:end); 
smb_5a = smb(61:end); 
umd_5a = umd(61:end);
rf_5a = rf(61:end); 

% --------------------------------------------------------
% form carhart-style portfolios using 5 years of alphas
% --------------------------------------------------------

'making carhart portfolios with 5 years of alphas'

f = [rmrf hml smb umd];
ta = 0; 
for pickdate = 60:12:(size(x,1)-12); % form portfolios on dec 31 of each date starting with end of year 6.
    ta = ta+1; 
    alive_5yr = sum(~isnan(x(pickdate-59:pickdate,:)))==60; % 1 if fund is alive at pickdate, and has 5 years of return history
    retsx_past = x(pickdate-59:pickdate,alive_5yr);
    f_past = f(pickdate-59:pickdate,:);
    %[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK] = tsregress2(retsx_past,f_past);
   
    X = [ones(size(retsx_past,1),1) f_past];
    b = X\retsx_past; 
    alpha = b(1,:)'; % rest of tsregress was taking forever
        
    lag_and_next = [alpha';x(pickdate+1:pickdate+12,alive_5yr)];  % keep next year with last 5 year alpha
    lag_and_next = sortrows(lag_and_next',1)'; % sort horizontally by lagged alpha
    breakpoints = floor((0.1:0.1:1)'*size(lag_and_next,2)); 
    breakpoints = [0;breakpoints]; % indices of funds at decile breakpoints
    for i = 1:10; % now form decile portfolio returns for the next year 
        alpha5yrkeep(ta,i) = mean(lag_and_next(1,breakpoints(i)+1:breakpoints(i+1))');
        for t = 1:12; 
            %sels = find(~isnan(lag_and_next(t+1,breakpoints(i)+1:breakpoints(i+1))')); % find non NAN returns in this month 
            %num_in_port_5alpha(pickdate+t,i) = length(sels); 
            carhart_port_5alpha(pickdate+t,i) = nanmean(lag_and_next(t+1,breakpoints(i)+1:breakpoints(i+1)),2);
        end; 
    end; 
end; 

% chop to available dates 

datesff_5alpha = datesff(61:end); 
carhart_port_5alpha = carhart_port_5alpha(61:end,:); 
rmrf_5alpha = rmrf(61:end); 
hml_5alpha = hml(61:end); 
smb_5alpha = smb(61:end); 
umd_5alpha = umd(61:end);
rf_5alpha = rf(61:end); 

'done making data'

% --------------------------
% place code from graveyard here to revive those ideas
% may need debugging on breakpoints and nanmean
% --------------------


% -------------------------------------------------------
% Start analysis of carhart data
% -----------------------------------------------------


disp('Carhart style data'); 
rx = [carhart_port carhart_port(:,end)-carhart_port(:,1)]; 
retkeep = [retkeep retkeep(:,10)-retkeep(:,1)];
T = size(carhart_port,1); 
f = [rmrf_1];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK,u ] = tsregress2(rx,f);
disp(' '); 
disp('Carhart portfolios'); 
disp(' '); 
fprintf('                  %7s %7s %7s %7s %7s %7s %7s %7s %7s %7s %7s \n','loser','2','3','4','5','6','7','8','9','winner','w-l');  
fprintf('    Mean lag retx %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', mean(retkeep)); 
fprintf('   Mean xs return %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', mean(rx)); 
fprintf('     std err Mean %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', std(rx)./T^0.5); 
fprintf('       CAPM alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', alpha'); 
fprintf('         t(alpha) %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', alpha'./siga'); 
fprintf('        CAPM beta %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', beta(:,1)'); 
fprintf('          CAPM R2 %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', R2'); 
fprintf(' Mean prev return %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', mean(retkeep)); 
fprintf('   Avg frac > mkt %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', mean(frac_beat_mkt)); 
fprintf('   Avg frac < mkt %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', mean(frac_less_mkt)); 

% ------------------
% returns and lagged returns 
% -------------------------

figure; 
plot(mean(retkeep),mean(rx),'v','markersize',5)
hold on;
axis([-1.5 2.5 -1.5 2.5]);
xlabel('Avg monthly return, last year'); 
ylabel('Avg Monthly return, next year'); 
text(mean(retkeep),mean(rx),{' loser',' 2',' 3',' 4',' 5',' 6',' 7',' 8',' 9',' winner', 'w-l'})

b = [ones(10,1) mean(retkeep(:,1:10))']\mean(rx(:,1:10))';
plot([-1.5 2.5]',[1 -1.5; 1 2.5]*b); 
text(0, -0.5, [num2str(b(1),'%4.2f') ' + ' num2str(b(2),'%4.2f') ' R_{t-1} ']); 


if 0;
    % this idea isn't working -- produces way too big slope
    b = [ones(size(biglag,2),1) biglag(:,1)']\biglag(:,2)'; % fund-level regression
    plot([-1.5 2.5]',[1 -1.5; 1 2.5]*b);
    text(0, -0.5, [num2str(b(1),'%4.2f') ' + ' num2str(b(2),'%4.2f') ' R_{t-1} '],'color','r');
   
    % and this one produces a negative slope. Time series and cross section
    % variation again!
    
    T = size(retkeep,1);
    K = size(retkeep,2);
    for t = 1:T
        retnext(t,:) = mean(rx(12*(t-1)+1:12*t,:)); % mean of carhart portfolio in the next year
    end;
    
    b = [ones(T*K,1) reshape(retkeep,T*K,1)]\reshape(retnext,T*K,1);
    plot([-1.5 2.5]',[1 -1.5; 1 2.5]*b,'-r');
    text(0, -0.5, [num2str(b(1),'%4.2f') ' + ' num2str(b(2),'%4.2f') ' R_{t-1} '],'color','r');   
end
print -depsc2 ps7a.eps;


fprintf('chi2 statistic, prob value       %7.2f %7.2f \n' , [chi2stat 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)    %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, prob value          %7.2f %7.2f \n', [Fstat 100*Fpv]);
fprintf('10, 5, 1%% p values of F          %7.2f %7.2f %7.2f \n', Fpvals);


f = [rmrf_1 hml_1 smb_1 ];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK, u] = tsregress2(rx,f);

disp('carhart portfolios, 3 factor model'); 
fprintf('         3F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', alpha'); 
fprintf('    t on 3F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', alpha'./siga'); 
fprintf('             3F b %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', beta(:,1)'); 
fprintf('             3F h %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', beta(:,2)'); 
fprintf('             3F s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', beta(:,3)'); 
fprintf('            3F R2 %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', R2'); 

fprintf('chi2 statistic, prob value       %7.2f %7.2f \n' , [chi2stat 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)    %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, prob value          %7.2f %7.2f \n', [Fstat 100*Fpv]);
fprintf('10, 5, 1%% p values of F          %7.2f %7.2f %7.2f \n', Fpvals);

disp('carhart portfolios, 4 factor model'); 
f = [rmrf_1 hml_1 smb_1 umd_1];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK, u] = tsregress2(rx,f);

fprintf('         4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', alpha'); 
fprintf('    t on 4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', alpha'./siga'); 
fprintf('             4F b %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', beta(:,1)'); 
fprintf('             4F h %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', beta(:,2)'); 
fprintf('             4F s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', beta(:,3)'); 
fprintf('             4F u %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', beta(:,4)'); 
fprintf('            4F R2 %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', R2'); 

fprintf('chi2 statistic, prob value       %7.2f %7.2f \n' , [chi2stat 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)    %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, prob value          %7.2f %7.2f \n', [Fstat 100*Fpv]);
fprintf('10, 5, 1%% p values of F          %7.2f %7.2f %7.2f \n', Fpvals);


% -------------------------------------------
% Do analysis of 5 - 0 year data 
% -------------------------------------------

disp('Carhart style data -- sort on year -5 to year 0 win/lose'); 

rx = [carhart_port_5a carhart_port_5a(:,end)-carhart_port_5a(:,1)]; 
T = size(carhart_port_5a,1); 
f = [rmrf_5a hml_5a smb_5a umd_5a];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK, u] = tsregress2(rx,f);
disp('5-0 portfolios, 4 factor model.'); 
fprintf('                  %7s %7s %7s %7s %7s %7s %7s %7s %7s %7s %7s \n','loser','2','3','4','5','6','7','8','9','winner','w-l');
fprintf('             mean %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', mean(rx)); 
fprintf('         4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', alpha'); 
fprintf('    t on 4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', alpha'./siga'); 
fprintf('             4F b %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,1)'); 
fprintf('             4F h %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,2)'); 
fprintf('             4F s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,3)'); 
fprintf('             4F u %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,4)'); 
fprintf('            4F R2 %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', R2'); 

fprintf('chi2 statistic, prob value       %7.2f %7.2f \n' , [chi2stat 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)    %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, prob value          %7.2f %7.2f \n', [Fstat 100*Fpv]);
fprintf('10, 5, 1%% p values of F          %7.2f %7.2f %7.2f \n', Fpvals);



% -------------------------------------------
% Do analysis of 5 - 0 year alpha sort data 
% -------------------------------------------

disp('Carhart style data -- sort on year -5 to year 0 4F alphas'); 

rx = carhart_port_5alpha; 
T = size(carhart_port_5alpha,1); 
f = [rmrf_5alpha hml_5alpha smb_5alpha umd_5alpha];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK, u] = tsregress2(rx,f);
disp('5-0 alpha portfolios, 4 factor model.'); 
  
fprintf('             mean %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', mean(rx)); 
fprintf('         4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', alpha'); 
fprintf('    t on 4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', alpha'./siga'); 
fprintf('             4F b %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,1)'); 
fprintf('             4F h %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,2)'); 
fprintf('             4F s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,3)'); 
fprintf('             4F u %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,4)'); 
fprintf('            4F R2 %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', R2'); 

fprintf('chi2 statistic, prob value       %7.2f %7.2f \n' , [chi2stat 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)    %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, prob value          %7.2f %7.2f \n', [Fstat 100*Fpv]);
fprintf('10, 5, 1%% p values of F          %7.2f %7.2f %7.2f \n', Fpvals);





% ----------------
% alpha t statistics
% ---------------

if 1; % for development, allows you to skip this while debugging below
    f = [rmrf hml smb umd];
    N = size(x,2);
    alpha = zeros(N,1);
    betas = zeros(4,N);
    alphat = zeros(N,1);
    for  i = 1:size(x,2); % loop through funds
        indx = find(~isnan(x(:,i))); % indices of valid dates
        if size(indx,1)>6;
            [alpha(i), betas(:,i), R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK] = ...
                tsregress2(x(indx,i),f(indx,:));
            alphat(i) = alpha(i)/siga;
        else;
            alpha(i) = NaN;
            alphat(i) = NaN;
            betas(:,i) = NaN*ones(4,1);
        end
    end;
    alphat = sort(alphat);
    alphat2 = alphat(find(~isnan(alphat)));
    alpha = sort(alpha);
    alpha2 = alpha(find(~isnan(alpha)));
    
    N2 = size(alphat2,1);
    
    medianalphat = alphat2(floor(N2/2));
    alpha5 = alphat2(floor(0.05*N2));
    alpha95 = alphat2(floor(0.95*N2));
    t5 = tinv(0.05,24);
    t50 = tinv(0.5,24);
    t95 = tinv(0.95,24);
    n5 = norminv(0.05);
    n50 = norminv(0.5);
    n95 = norminv(0.95);
   
    disp(' alpha 5, 50 95');
    disp([alpha5 medianalphat alpha95]);
    disp(' t 5, 95, 24 degrees of freedom');
    disp([t5 t50 t95]);
    disp(' N 5, 50, 95, ');
    disp([n5 n50 n95]);
    
    
    figure;
    plot(alphat2,(1:N2)'/N2*100,'Linewidth',2);
    hold on
    plot([-10 10],[0 0],'-k')
    plot(medianalphat*[1 1],[-10 110],'-k');
    text(medianalphat+.1,20,['median = ' num2str(medianalphat,'%3.2f')]);
    plot(alpha5*[1 1],[-10 110],'-k');
    text(alpha5+.1,90,['5% = ' num2str(alpha5,'%3.2f')]);
    plot(t5*[1 1],[-10 110],'-r');
    text(t5+.1,80,['t 5% = ' num2str(t5,'%3.2f')],'color','r');
    plot(alpha95*[1 1],[-10 110],'-k');
    text(alpha95+.1,70,['95% = ' num2str(alpha95,'%3.2f')]);
    plot(t95*[1 1],[-10 110],'-r');
    text(alpha95+.1,60,['t 95% = ' num2str(t95,'%3.2f')],'color','r');
    plot(alphat2,100*tcdf(alphat2,24),'-r','linewidth',2)
    axis([-5 5 -5 105])
    title('alpha t statistic');
    print -depsc2 ps7_alphatstat.eps;
   
    
   % version that compares to normal rather than t(24) 
    figure;
    plot(alphat2,(1:N2)'/N2*100,'Linewidth',2);
    hold on
    plot([-10 10],[0 0],'-k')
    plot(medianalphat*[1 1],[-10 110],'-k');
    text(medianalphat+.1,20,['median = ' num2str(medianalphat,'%3.2f')]);
    plot(alpha5*[1 1],[-10 110],'-k');
    text(alpha5+.1,90,['5% = ' num2str(alpha5,'%3.2f')]);
    plot(n5*[1 1],[-10 110],'-r');
    text(n5+.1,80,['N 5% = ' num2str(n5,'%3.2f')],'color','r');
    plot(alpha95*[1 1],[-10 110],'-k');
    text(alpha95+.1,70,['95% = ' num2str(alpha95,'%3.2f')]);
    plot(t95*[1 1],[-10 110],'-r');
    text(alpha95+.1,60,['N 95% = ' num2str(n95,'%3.2f')],'color','r');
    plot(alphat2,100*normcdf(alphat2),'-r','linewidth',2)
    axis([-5 5 -5 105])
    title('alpha t statistic');
    print -depsc2 ps7_alphatstat_normal.eps;
    
    
    
    medianalpha = alpha2(floor(N2/2));
    alpha5 = alpha2(floor(0.05*N2));
    alpha95 = alpha2(floor(0.95*N2));
    disp(' alpha 5, 50 95');
    disp([alpha5 medianalpha alpha95]);
    
    figure;
    plot(alpha2,(1:N2)'/N2*100,'Linewidth',2);
    hold on
    plot(alpha2,0*(1:N2)'/N2*100,'-k');
    plot(medianalpha*[1 1],[-10 110],'-k');
    text(medianalpha+.1,20,['median = ' num2str(medianalpha,'%3.2f')]);
    plot(alpha5*[1 1],[-10 110],'-k');
    text(alpha5+.1,50,['5% = ' num2str(alpha5,'%3.2f')]);
    plot(alpha95*[1 1],[-10 110],'-k');
    text(alpha95+.1,10,['95% = ' num2str(alpha95,'%3.2f')]);
    axis([-1.5 1.5 -5 105])
    title('alpha, percent per month')
    print -depsc2 ps7_alpha.eps;
 
    
end; % ends if 0



stop;


% ---------------------
% graveyard of old carhart portfolio ideas
% ---------------------


% --------------------------------------------------------
% form carhart-style portfolios using 5-1 returns. 
% --------------------------------------------------------

ta = 0; 
for pickdate = 60:12:(size(x,1)-12); % form portfolios on dec 31 of each date starting with end of year 6.
    ta = ta+1; 
    alive_5yr = sum(~isnan(x(pickdate-59:pickdate,:)))==60; % 1 if fund is alive at pickdate, and has 5 years of return history
    ret5yr = mean(x(pickdate-59:pickdate-12,alive_5yr)); % return from t-5 to t-1 
    lag_and_next = [ret5yr;x(pickdate+1:pickdate+12,alive_5yr)];  % keep next year with last 5 year average
    lag_and_next = sortrows(lag_and_next',1)'; % sort horizontally by lagged average return
    breakpoints = floor((0.1:0.1:1)'*size(lag_and_next,2)); 
    breakpoints = [1;breakpoints]; % indices of funds at decile breakpoints
    for i = 1:10; % now form decile portfolio returns for the next year 
        ret5yrkeep(ta,i) = mean(lag_and_next(1,breakpoints(i):breakpoints(i+1))');
        for t = 1:12; 
            sels = find(~isnan(lag_and_next(t+1,breakpoints(i):breakpoints(i+1))')); % find non NAN returns in this month 
            num_in_port_5(pickdate+t,i) = length(sels); 
            carhart_port_5(pickdate+t,i) = mean(lag_and_next(t+1,breakpoints(i)+sels-1)')';
        end; 
    end; 
end; 

% chop to available dates 

datesff_5 = datesff(61:end); 
carhart_port_5 = carhart_port_5(61:end,:); 
rmrf_5 = rmrf(61:end); 
hml_5 = hml(61:end); 
smb_5 = smb(61:end); 
umd_5 = umd(61:end);
rf_5 = rf(61:end); 

% write out data set

if print_data; 
    fid = fopen('carhart_portfolios_5.txt','w'); 
    fprintf(fid,'%% Carhart-style portfolios \r');
    fprintf(fid,'%% all stock funds are sorted in to portfolios based on year - 5 to year - 1 return ending in dec. Then we track  \r');
    fprintf(fid,'%% monthly returns of the portfolios for the next year. \r');
    fprintf(fid,'%% all returns are percent monthly. These are ALREADY EXCESS -- no need to subtract rf  \r'); 
    fprintf(fid,'%13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s \r','% Date','Past loser','2','3','4','5','6','7','8','9','Winner','rmrf','hml','smb','umd','rf');
    fprintf(fid,'%13i %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f \r', ...
        [datesff_5 carhart_port_5 rmrf_5 hml_5 smb_5 umd_5 rf_5]'); 
    fclose(fid);
    dlmwrite('carhart_portfolios_5.csv',[datesff_5 carhart_port_5 rmrf_5 hml_5 smb_5 umd_5 rf_5]); 
end;
   


% --------------------------------------------------------
% form carhart-style portfolios using 5 years of returns and extra filters. 
% --------------------------------------------------------

ta = 0; 
for pickdate = 60:12:(size(x,1)-12); % form portfolios on dec 31 of each date starting with end of year 6.
    ta = ta+1; 
    alive_5yr = sum(~isnan(x(pickdate-59:pickdate,:)))==60; % 1 if fund is alive at pickdate, and has 5 years of return history
    
    ret5yr = x(pickdate-59:pickdate,alive_5yr);
    ret5yr = ret5yr - rmrf(pickdate-59:pickdate)*ones(1,size(ret5yr,2)); % subtract market returns; 
    tests1 = mean(ret5yr)'; 
    tests2 = mean(ret5yr)'./(std(ret5yr)*5^0.5)' ;
    tests3 = sum(ret5yr>0)';
    tests_all =  tests2./std(tests2)+tests3./std(tests3);  % divide by std dev to give equal weight
        
    lag_and_next = [tests_all';x(pickdate+1:pickdate+12,alive_5yr)];  % keep next year with last 5 year average
    lag_and_next = sortrows(lag_and_next',1)'; % sort horizontally by lagged average return
    breakpoints = floor((0.1:0.1:1)'*size(lag_and_next,2)); 
    breakpoints = [1;breakpoints]; % indices of funds at decile breakpoints
    for i = 1:10; % now form decile portfolio returns for the next year 
        ret5yrkeep(ta,i) = mean(lag_and_next(1,breakpoints(i):breakpoints(i+1))');
        for t = 1:12; 
            sels = find(~isnan(lag_and_next(t+1,breakpoints(i):breakpoints(i+1))')); % find non NAN returns in this month 
            num_in_port_5b(pickdate+t,i) = length(sels); 
            carhart_port_5b(pickdate+t,i) = mean(lag_and_next(t+1,breakpoints(i)+sels-1)')';
        end; 
    end; 
end; 

% chop to available dates 

datesff_5b = datesff(61:end); 
carhart_port_5b = carhart_port_5b(61:end,:); 
rmrf_5b = rmrf(61:end); 
hml_5b = hml(61:end); 
smb_5b = smb(61:end); 
umd_5b = umd(61:end);
rf_5b = rf(61:end); 

% write out data set

if print_data; 
    fid = fopen('carhart_portfolios_5b.txt','w'); 
    fprintf(fid,'%% Carhart-style portfolios -- 5 year returns, 3 sorts\r');
    fprintf(fid,'%13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s \r','% Date','Past loser','2','3','4','5','6','7','8','9','Winner','rmrf','hml','smb','umd','rf');
    fprintf(fid,'%13i %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f %13.8f \r', ...
        [datesff_5b carhart_port_5b rmrf_5b hml_5b smb_5b umd_5b rf_5b]'); 
    fclose(fid);
    dlmwrite('carhart_portfolios_5b.csv',[datesff_5b carhart_port_5b rmrf_5b hml_5b smb_5b umd_5b rf_5b]); 
end;



% -------------------------------------------
% Do analysis of 5 - 1 year data 
% -------------------------------------------

disp('Carhart style data -- sort on year -5 to year -1 year win/lose'); 

rx = carhart_port_5; 
T = size(carhart_port_5,1); 
f = [rmrf_5 hml_5 smb_5 umd_5];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK] = tsregress2(rx,f);
disp('5-1 portfolios, 4 factor model.'); 
fprintf('             mean %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', mean(rx)); 
fprintf('         4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', alpha'); 
fprintf('    t on 4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', alpha'./siga'); 
fprintf('             4F b %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,1)'); 
fprintf('             4F h %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,2)'); 
fprintf('             4F s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,3)'); 
fprintf('             4F u %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,4)'); 
fprintf('            4F R2 %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', R2'); 

fprintf('chi2 statistic, prob value       %7.2f %7.2f \n' , [chi2stat 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)    %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, prob value          %7.2f %7.2f \n', [Fstat 100*Fpv]);
fprintf('10, 5, 1%% p values of F          %7.2f %7.2f %7.2f \n', Fpvals);



% -------------------------------------------
% Do analysis of 5 - 0 year data 
% -------------------------------------------

disp('Carhart style data -- sort on year -5 to year 0 win/lose with extra filters'); 

rx = carhart_port_5b; 
T = size(carhart_port_5b,1); 
f = [rmrf_5b hml_5b smb_5b umd_5b];
[alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N, Fstat, Fpvals, Fpv, TNK] = tsregress2(rx,f);
disp('5-0 portfolios with filters, 4 factor model.'); 
fprintf('             mean %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', mean(rx)); 
fprintf('         4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', alpha'); 
fprintf('    t on 4F alpha %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', alpha'./siga'); 
fprintf('             4F b %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,1)'); 
fprintf('             4F h %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,2)'); 
fprintf('             4F s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,3)'); 
fprintf('             4F u %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', beta(:,4)'); 
fprintf('            4F R2 %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', R2'); 

fprintf('chi2 statistic, prob value       %7.2f %7.2f \n' , [chi2stat 100*chi2pv]);
fprintf('10, 5, 1%% p values of chi2(N)    %7.2f %7.2f %7.2f \n', chi2vals);
fprintf('F statistic, prob value          %7.2f %7.2f \n', [Fstat 100*Fpv]);
fprintf('10, 5, 1%% p values of F          %7.2f %7.2f %7.2f \n', Fpvals);
