function  [bb2 bsum2] = asness(rx,f,rf,lbls,flbls,lead,lag)
% rx = excess returns
% f = factors (do not include constant)
% lbls = return lables, e.g. {'first return','second return'}
% flbls = factor lables, e.g. {'rmrf','hml'}
% lead, lag = nubmer of leads and lags in time series regression 0 = contemporaneous only
% plsmin = 1 if you want separate up and down betas

T = size(rx,1);
if T ~= size(f,1);
    disp('error -- factors and returns are different lengths');
end;
if any(any(isnan(rx))');
    disp('missing values in rx');
end;
if any(any(isnan(f))');
    disp('missing values in f');
end;
K = size(f,2);
N = size(rx,2);

if lead>0;
    disp('NEED TO REPROGRAM TO ALLOW LEADS');
end;

% do regression with no leads/lags or up/down betas as a benchmark

X = [ones(T-lag,1) f(1+lag:end,:)]; % keep sample same. NO leads allowed by the way.
b = X\rx(1+lag:end,:);
e = rx(1+lag:end,:)-X*b;
se = var(e).^0.5;
R2 = 1-var(e)./var(rx(1+lag:end,:));
seb = (diag(inv(X'*X))*(var(e))).^0.5;

disp(' ' );
fprintf(' %20s %10s \n',' ','Monthly returns');
fprintf(' %20s %6s %6s %6s %6s (%6s) \n','  ', 'E(rx)','s(rx)','Sharpe','Annual','tstat');
for i = 1:N;
    fprintf(' %20s ', char(lbls(i)));
    fprintf('%6.2f ',mean(rx(:,i)), std(rx(:,i)), [1 12^0.5].*mean(rx(:,i))/std(rx(:,i)) );
    fprintf('(%6.2f) \n',mean(rx(:,i))/(std(rx(:,i))/(size(rx(:,i),1)^0.5)));
end;
for i = 1:K;
    fprintf(' %20s ', char(flbls(i)));
    fprintf('%6.2f ',mean(f(:,i)), std(f(:,i)), [1 12^0.5].*mean(f(:,i))/std(f(:,i)) );
    fprintf('(%6.2f) \n',mean(f(:,i))/(std(f(:,i))/(size(f(:,i),1)^0.5)));
end;


disp(' ' );
fprintf(' %20s %10s \n',' ','betas');
fprintf(' %20s %6s %6s %6s ','  ', 'E(rx)','alpha','t(a)');
for i = 1:K;
    fprintf('%6s ', char(flbls(i)))
end;
fprintf('%6s %6s %6s \n','R2','a-Sharpe','s(e)')
for i = 1:N;
    fprintf(' %20s ', char(lbls(i)));
    fprintf('%6.2f ',mean(rx(1+lag:end,i)), b(1,i));
    fprintf('(%5.2f)',b(1,i)/seb(1,i));
    fprintf('%6.2f ',b(2:end,i));
    fprintf('%6.2f %6.2f %6.2f\n', R2(i),b(1,i)/se(i),se(i));
end;

if 1 % standard errors of betas can be commented out to keep output from being overwhelming
    fprintf(' %20s %15s \n',' ','std errors');
    fprintf(' %20s %6s %6s','  ', 'E(rx)','alpha');
    for i = 1:K;
        fprintf('%6s ', char(flbls(i)))
    end;
    fprintf('%6s \n','R2')
    for i = 1:N;
        fprintf(' %20s ', char(lbls(i)));
        fprintf('%6.2f ',std(rx(1+lag:end,i))/(T^0.5), seb(1,i));
        fprintf('%6.2f ',seb(2:end,i));
        fprintf('%6.2f \n', R2(i));
    end;
end


% now lead and lagged betas too.  -- lead=0, lag = 0 corresponds to only contemporaneous betas
if ((lead > 0) | (lag > 0)) ;
    
    rxb = rx(1+lag:T-lead,:); % leave room for leads and lags in regressions
    premium = mean(f); % saved to multiply betas
    fb = f(1+lag:T-lead,:);
    for i = 1:lag;  % add lead and lagged factors
        fb = [f(1+lag-i:T-lead-i,:) fb]; % organization of f  : 1 rmrf-lag hml-lag....rmrf-1 hml-1...rmrf hml...rmrf+1 hml+1 ...r30y+lead
    end;
    for i = 1:lead;
        fb = [fb f(1+lag+i:T-lead+i,:) ];
    end;
    
    X = [ones(T-lead-lag,1) fb];
    bb = X\rxb;
    e = rxb-X*bb;
    se = var(e).^0.5;
    R2b = 1-var(e)./var(rxb);
    seb = diag(inv(X'*X))*var(e); % many regressions at once here.
    alphas = bb(1,:);
    talphas = alphas./seb(1,:).^0.5;
    bb = bb(2:end,:);
    bsum = bb(1:K,:);
    for i = 1:lag+lead;
        bsum = bsum+bb(i*K+1:i*K+K,:);
    end;
    
    disp(' ' );
    fprintf(' %20s %10s %6.2f %10s %6.2f %10s\n',' ',' betas with', lead, ' leads and', lag, 'lags');
    fprintf(' %20s %6s %6s %6s','  ', 'E(Rx)','alpha','t(a)');
    for i = 1:K;
        fprintf('%6s ', char(flbls(i)))
    end;
    fprintf('%6s %6s %6s \n','R2','a-Sharpe','s(e)')
    Erx = mean(rxb);
    for i = 1:N;
        fprintf(' %20s ', char(lbls(i)));
        fprintf('%6.2f ',Erx(i), alphas(i));
        fprintf('(%5.2f)',talphas(i));
        fprintf('%6.2f ',bsum(:,i));
        fprintf('%6.2f %6.2f %6.2f \n', R2b(i),alphas(i)/se(i),se(i))
    end;
    fprintf(' %20s %6s %6s %6s','  ', 'E(Rx)','alpha','t(a)');
    for i = 1:K;
        fprintf('%6s ', char(flbls(i)) ) ; % eg rmrf
        for j = 1:lag;
            fprintf('  -%1i   ', j ); % eg -1
        end;
    end;
    fprintf('%6s \n','R2')
    for i = 1:N;
        fprintf(' %20s ', char(lbls(i)));
        fprintf('%6.2f ',Erx(i), alphas(i));
        fprintf('(%5.2f)',talphas(i));
        % bb is K*(lag+1) long. thus rmrf_0 is in the K*(lag)+1 position
        indis = K*(lag:-1:0)+1; % indices of rmrf_0 rmrf_-1,....
        for j = 2:K;
            indis = [ indis K*(lag:-1:0)+j ];
        end;
        %indis
        fprintf('%6.2f ',bb(indis,i));
        fprintf('%6.2f \n', R2b(i))
    end;
    
    
    disp(' ' );
    fprintf(' %20s %10s %6.2f %10s %6.2f %10s\n',' ',' betas x premiums/E(rx) with', lead, ' leads and', lag, 'lags');
    fprintf(' %20s %6s %6s %6s','  ', 'E(Rx)','alpha','t(a)');
    for i = 1:K;
        fprintf('%6s ', char(flbls(i)))
    end;
    fprintf('%6s %6s %6s \n','R2','a-Sharpe','s(e)')
    Erx = mean(rxb);
    for i = 1:N;
        fprintf(' %20s ', char(lbls(i)));
        fprintf('%6.2f ',Erx(i), alphas(i));
        fprintf('(%5.2f)',talphas(i));
        fprintf('%6.2f ',bsum(:,i).*premium'/Erx(i)); %% rest of this section is all copied from above except this line
        fprintf('%6.2f %6.2f %6.2f \n', R2b(i),alphas(i)/se(i),se(i))
    end;
    fprintf(' %20s %6s %6s %6s','  ', 'E(Rx)','alpha','t(a)');
    for i = 1:K;
        fprintf('%6s ', char(flbls(i)) ) ; % eg rmrf
        for j = 1:lag;
            fprintf('  -%1i   ', j ); % eg -1
        end;
    end;
    fprintf('%6s \n','R2')
    for i = 1:N;
        fprintf(' %20s ', char(lbls(i)));
        fprintf('%6.2f ',Erx(i), alphas(i));
        fprintf('(%5.2f)',talphas(i));
        % bb is K*(lag+1) long. thus rmrf_0 is in the K*(lag)+1 position
        indis = K*(lag:-1:0)+1; % indices of rmrf_0 rmrf_-1,....
        for j = 2:K;
            indis = [ indis K*(lag:-1:0)+j ];
        end;
        %indis
        fprintf('%6.2f ',bb(indis,i));
        fprintf('%6.2f \n', R2b(i))
    end;
    
    
    
end;


% up/down betas

rxb = rx(1+lag:T-lead,:); % leave room for leads and lags in regressions
% organization of f  : 1 rmrf-lag+ rmrf-lag-  hml-lag+ hml-lag-.......rmrf-1+ rmrf-1- ...

fb = f(1+lag:T-lead,:);
X = [ones(T-lead-lag,1) fb.*(fb>=0) fb.*(fb<0)];
bb2 = X\rxb;
e = rxb-X*bb2;
R2b2 = 1-var(e)./var(rxb);
%seb = (diag(inv(f'*f))*(var(e))).^0.5;
alphas = bb2(1,:);
bb2 = bb2(2:end,:);

if lead > 0 ;
    disp('error: not programmed for leads') ;
end;

disp(' ' );
fprintf(' %20s up/down betas with 0 leads and lags \n',' ');
fprintf(' %20s %6s %6s ','  ', 'E(Rx)','alpha');
for i = 1:K;
    fprintf('%6s ', [char(flbls(i)) '+'])
    fprintf('%6s ', [char(flbls(i)) '-'])
end;
fprintf('%6s \n','R2')
Erx = mean(rxb);
for i = 1:N;
    fprintf(' %20s ', char(lbls(i)));
    fprintf('%6.2f ',Erx(i), alphas(i));
    % for k = 3 want b in the order 1 4 2 5 3 6 so up and down alternate
    indx = [ 1 K+1];
    for j = 2:K;
        indx = [ indx j j+K ];
    end;
    fprintf('%6.2f ',bb2(indx,i));
    fprintf('%6.2f \n', R2b2(i))
end;




% finally, lead lag and up/down
if 1
    
    rxb = rx(1+lag:T-lead,:); % leave room for leads and lags in regressions
    % organization of f  : 1 rmrf-lag+ rmrf-lag-  hml-lag+ hml-lag-.......rmrf-1+ rmrf-1- ...
    
    fb = f(1+lag:T-lead,:);
    for i = 1:lag;  % add lead and lagged factors
        fb = [f(1+lag-i:T-lead-i,:) fb]; % organization of f  : 1 rmrf-lag hml-lag....rmrf-1 hml-1...rmrf hml...rmrf+1 hml+1 ...r30y+lead
    end;
    for i = 1:lead;
        fb = [fb f(1+lag+i:T-lead+i,:) ];
    end;
    
    X = [ones(T-lead-lag,1) fb.*(fb>=0) fb.*(fb<0)];
    % for capm orgainzation is rm-3+ rm-2+ rm-1+ rm+ rm-3- rm-2- rm-1- rm-
    % for ff3 organization is rm-3+ hml-3+ smb-3+ rm-2+ hml-2+ smb-2+ ...
    bb2 = X\rxb;
    e = rxb-X*bb2;
    R2b2 = 1-var(e)./var(rxb);
    %seb = (diag(inv(f'*f))*(var(e))).^0.5;
    alphas = bb2(1,:);
    bb2 = bb2(2:end,:);
    
    if lead > 0 ;
        disp('error: not programmed for leads') ;
    end;
    % diagnostic to see if this works.
    %disp('bb2 lags added to sum');
    %disp([(1:K) ((lag+1)*K+1:(lag+1)*K+K)])
    bsum2 = bb2([(1:K) ((lag+1)*K+1:(lag+1)*K+K)],:);
    for i = 1:lag+lead;
        %disp([i*K+(1:K) i*K+((lag+1)*K+1:(lag+1)*K+K)])
        bsum2 = bsum2 + bb2([i*K+(1:K) i*K+((lag+1)*K+1:(lag+1)*K+K)],:);
    end;
    
    disp(' ' );
    fprintf(' %10s %20s %6.2f %10s %6.2f %10s \n',' ','up/down betas with', lead, 'leads and', lag, 'lags');
    fprintf(' %20s %6s %6s ','  ', 'E(Rx)','alpha');
    for i = 1:K;
        fprintf('%6s ', [char(flbls(i)) '+'])
        fprintf('%6s ', [char(flbls(i)) '-'])
    end;
    fprintf('%6s \n','R2')
    Erx = mean(rxb);
    for i = 1:N;
        fprintf(' %20s ', char(lbls(i)));
        fprintf('%6.2f ',Erx(i), alphas(i));
        % for k = 3 want b in the order 1 4 2 5 3 6 so up and down alternate
        indx = [ 1 K+1];
        for j = 2:K;
            indx = [ indx j j+K ];
        end;
        fprintf('%6.2f ',bsum2(indx,i));
        fprintf('%6.2f \n', R2b2(i))
    end;
    
    if (size(f,2) == 1) & (lag==3);
        disp('                    showing all  lags')
        fprintf(' %20s %6.6s %6.6s %6.6s %6.6s %6.6s %6.6s %6.6s %6.6s\n',' ','rm-3+','rm-2+','rm-1+','rm+','rm-3-','rm-2-','rm-1-1','rm-');
        for i = 1:N;
            fprintf(' %20s ', char(lbls(i)));
            fprintf('%6.2f ',bb2(:,i));
            fprintf('\n');
        end;
        
        disp('                     showing all lags together and current ')
        fprintf(' %20s %6.6s %6.6s %6.6s %6.6s \n',' ','rlag+','rm+','rlag-','rm-');
        for i = 1:N;
            fprintf(' %20s ', char(lbls(i)));
            fprintf('%6.2f ',sum(bb2(1:3,i)), bb2(4,i), sum(bb2(5:7,i)), bb2(8,i));
            fprintf('\n');
        end;
        
    end;
    
    % scatter plot with up/down betas too
    if (K == 1) & (N<5) & (lead==0) & (lag==0) ;
        for i = 1:N;
            figure;
            plot(f,rx(:,i),'.');
            hold on;
            xy = [f alphas + f.*(f>=0)*bsum2(1) + f.*(f<0)*bsum2(2)];
            xy = sort(xy,1);
            plot(xy(:,1),xy(:,2),'-k')
            title(['scatterplot for ' char(lbls(i)) ' ' char(flbls(1))]);
            print('-depsc2',['scatter_' char(lbls(i)) '.eps']);
        end;
    end;
    
    % annual horizon betas
    disp(' ');
    rfa = 100*( exp(filter(ones(12,1),1,log(1+rf/100))) -1 );
    X = [ones(T,1) 100*(exp(filter(ones(12,1),1,log(1+(f+rf*ones(1,size(f,2)))/100)))-1)-rfa*ones(1,size(f,2))]; % keep sample same. NO leads allowed by the way.
    y=100*(exp(filter(ones(12,1),1,log(1+(rx+rf*ones(1,size(rx,2)))/100)))-1)-rfa*ones(1,size(rx,2));
    X = X(13:end,:);
    y = y(13:end,:);
    [b,seb,R2,R2vadj,v,F] = olsgmm(y,X,12,0); % necessary with overlapping data
    %b = X\y;
    e = y-X*b;
    %se = var(e).^0.5;
    %R2 = 1-var(e)./var(y);
    %seb = (diag(inv(X'*X))*(var(e))).^0.5;
    
    disp(' ' );
    fprintf(' %20s %10s \n',' ','Annual betas');
    fprintf(' %20s %6s %6s %6s %6s ','  ', 'E(rx)','s(rx)','alpha','t(a)');
    for i = 1:K;
        fprintf('%6s ', char(flbls(i)))
    end;
    fprintf('%6s %6s %6s \n','R2','a-Sharpe','s(e)')
    for i = 1:N;
        fprintf(' %20s ', char(lbls(i)));
        fprintf('%6.2f ',mean(y(:,i)),std(y(:,i)), b(1,i));
        fprintf('(%5.2f)',b(1,i)/seb(1,i));
        fprintf('%6.2f ',b(2:end,i));
        fprintf('%6.2f %6.2f %6.2f \n', R2(i),b(1,i)/se(i),std(e(:,i)));
    end;
    
    
    
end;


