function [alpha, beta, R2, siga, sigb, chi2stat, chi2vals, chi2pv, N2, Fstat, Fpvals, Fpv, TNK, u] = tsregress2(rx,f);
% function tsregress. Does time series regressions
% inputs: 
% rx: T x N matrix of excess returns in percent units for test assets
% f: T x K matrix of factors 
% 
% outputs: 
% alpha N x 1 vector of alphas
% beta N x K vector of betas
% R2 N x 1 vector of time series R2 
% siga = N x 1 vector of standard errors of alpha
% sigb = N x K matrix of standard errors of betas
% chi2stat = chi2 statistic for all alphas = 0
% chi2 vals = 3 x 1 vector of 10 5 1 % probability values for chi2
%            statistic, for comparison
% chi2pv = probability value of chi2 statistic .. chance of seeing a larger
%     value if true alphas are zero 
% N2 = degrees of freedom for chi2 distribution
% Fstat = F statistic for all alphas = 0
% Fvals = 3 x 1 vector of 10 5 1 % probability values for F
% Fpv = probability value of F statistic
% TNK = T - N2 - K = other degrees of freedom for the F statistic
% u = residuals. 

T = size(rx,1); 
N = size(rx,2); 
K = size(f,2); 
if size(f,1) ~= size(rx,1); 
    disp('tsregress error: factors and returns must be of same length'); 
end; 

X = [ones(T,1) f];
b = X\rx; 
u = rx-X*b; 
sigma = cov(u,1);              % use 1/T in covariance matrix as per GRS test, not 1/(T-1) or 1/(T-K)
s2 = diag(sigma);     
mx = inv(X'*X); 
dmx = diag(mx); 
siga = (dmx(1)*s2).^0.5;        % standard errors of alpha, beta 
sigb = (s2*dmx(2:end)').^0.5;   % s2 is a column vector of 25. dmx' is a row vector corresponding to factors. 
                                % this produces a matrix the same size as the b3 coefficients. 
                                
                                % explanation. If s2 is error variance in one regression, inv(X'*X)*s2 produces the covariance matrix of 
                                % the regression coefficients. Our trouble is we have say, 25, regressions going on side by side. We can't
                                % (easily) produce 25 matrices side by side. But we're not really interested in the whole covariance matrix.
                                % I.e. we want [var(alpha)]
                                %              [var(beta) ]
                                % for each regression. We don't want all of what inv(X'*X)*s2 produces, namely
                                %              [ var(alpha) cov(alpha,beta)] 
                                %              [ cov(alpha,beta) var(beta) ] 
                                % so , how do we produce 
                                % [var(alpha1) var(alpha2) ....var(alpha25) ] 
                                % [ var(beta1) var(beta2)  ....var(beta25)  ] 
                                % well, for a single regression diag(inv(X'*X))*s2 would produce waht we want. So all this produces 
                                % diag(inv(X'*X))*[s2_1  s2_2 .....s2_25]. 

alpha = b(1,:)'; 
beta = b(2:end,:)';           % is now 25 x 3 with portfolios up and down 
R2 = 1-s2./(std(rx).^2)'; 

% ---------------------------
% chi2 and F tests on alphas
% ---------------------------

Ef = mean(f)'; 
sigmaf = cov(f,1); 

N2 = rank(sigma); 

chi2stat = T*(1+Ef'*pinv(sigmaf)*Ef)^-1*alpha'*pinv(sigma)*alpha;  
chi2vals = chi2inv([0.9 0.95 0.99],N2); % makes the cutoff points for 10,5,1% tests
chi2pv = 1-chi2cdf(chi2stat,N2);          % displays the probability value for this test

Fstat = (T-N2-K)/N2*(1+Ef'*pinv(sigmaf)*Ef )^(-1)*alpha'*pinv(sigma)*alpha; 
Fpvals = finv([0.9 0.95 0.99],N2,T-N2-K); 
Fpv = 1-fcdf(Fstat,N2,T-N2-K); 

TNK = T-N2-K; 

% note N2 and pinv allow for singular (less than full rank) sigma matrices.
% This happens if some of the assets are redundant, for example, if you
% test the model using R1, R2 and (R1+R2) -- N should really be 2 in this
% case. It also happens if the factor is a linear combination of the test
% assets. For example, if the market portfolio is the sum of the test asset
% returns, then the equally weighted portfolio of tests assets has no
% residual in the regression Re = a + bRem + error. This means there is a
% linear combination of the test assets with no error, i.e. the error
% covariance matrix is not full rank. pinv is a "psuedo inverse" that can
% invert less than full rank matrices. In these situations, there are
% really less degrees of freedom than you think, so we use N2 rather than
% N. I've avoided this situation in class and in the problem set, 
% but it will come up in the real world and a later problem set. 

return; 