Rets = dlmread("~/data/2015/french_data/10_Industry_Portfolios_monthly_ew.txt","",12,0);
Factors = dlmread ("~/data/2015/french_data/F-F_Research_Data_Factors_monthly.txt","",4,0);

RMRF = Factors(:,2);
Rf = Factors(:,5);
                             # subtract Rf from each row
eR = Rets(:,2)-Rf;
eR =[eR Rets(:,3)-Rf];
eR =[eR Rets(:,4)-Rf];
eR =[eR Rets(:,5)-Rf];
eR =[eR Rets(:,6)-Rf];
eR =[eR Rets(:,7)-Rf];
eR =[eR Rets(:,8)-Rf];
eR =[eR Rets(:,9)-Rf];
eR =[eR Rets(:,10)-Rf];
eR =[eR Rets(:,11)-Rf];
                           # challenge: find a more compact way
T = rows(eR)
N = columns(eR)

eR = eR/100.0;
eRm  = RMRF/100.0;
[b1, sigma1, e1] = ols(eR(:,1),[ones(T,1) eRm] );
[b2, sigma2, e2] = ols(eR(:,2),[ones(T,1) eRm] );
[b3, sigma3, e3] = ols(eR(:,3),[ones(T,1) eRm] );
[b4, sigma4, e4] = ols(eR(:,4),[ones(T,1) eRm] );
[b5, sigma5, e5] = ols(eR(:,5),[ones(T,1) eRm] );
[b6, sigma6, e6] = ols(eR(:,6),[ones(T,1) eRm] );
[b7, sigma7, e7] = ols(eR(:,7),[ones(T,1) eRm] );
[b8, sigma8, e8] = ols(eR(:,8),[ones(T,1) eRm] );
[b9, sigma9, e9] = ols(eR(:,9),[ones(T,1) eRm] );
[b10, sigma10, e10] = ols(eR(:,10),[ones(T,1) eRm] );

alpha=[b1(1),b2(1),b3(1),b4(1),b5(1),b6(1),b7(1),b8(1),b9(1),b10(1)]';
e = [e1, e2, e3, e4, e5, e6, e7 ,e8, e9, e10];
Sigma = cov(e)
SigmaInv = inv(Sigma)
Sm = mean(eRm)/std(eRm)
J0 = T * inv(1+Sm^2)* (alpha')*SigmaInv*alpha
J1 = ((T-N-1)/N) * inv(1+Sm^2)*alpha'*SigmaInv*alpha
