ex3_finecx.m

计算最小二乘法ordinary least squares

clear all;
clc;

load econm2029_US_industries.mat;

whos;

T = size(r,1);
K = size(r,2);


% Create excess monthly returns

xsr = r - repmat(rf,1,K);
xsmkt = mkt;


% Plot out mean monthly returns and return standard deviations

xsrMean = nanmean(xsr);
xsrStd = std(xsr);

xsmktMean = nanmean(xsmkt);
xsmktStd = std(xsmkt);

figure;
bar([xsrMean xsmktMean])
title('Mean monthly returns: 30 industry portfolios and market');

figure;
bar([xsrStd xsmktStd])
title('Monthly return standard deviations: 30 industry portfolios and market');


% Scatter plot means and standard deviations

figure;
plot([xsrStd xsmktStd],[xsrMean xsmktMean],'k+');


% Run regressions for each of the 30 portfolios and store alphas and betas

x = [ones(T,1) xsmkt];
alphas = zeros(1,K);
betas = zeros(1,K);
alphaTStats = zeros(1,K);
betaTStats = zeros(1,K);
residuals = zeros(T,K);

fprintf('\n\nIndustry return regression analysis:\n');

for i=1:K;
    
    tmpcoeff = inv(x'*x)*x'*xsr(:,i);
    tmpres = xsr(:,i) - x*tmpcoeff;
    tmps2 = sum(tmpres.^2)/(T-1);
    tmpvcv = tmps2*inv(x'*x);
    tmpse = sqrt(diag(tmpvcv));
    
    residuals(:,i) = tmpres;
    
    alphas(i) = tmpcoeff(1);
    betas(i) = tmpcoeff(2);
    alphaTStats(i) = tmpcoeff(1)./tmpse(1);
    betaTStats(i) = tmpcoeff(2)./tmpse(2);

    fprintf('Portfolio %g: alpha (t-stat) : %4.2f (%4.2f): beta (t-stat) : %4.2f (%4.2f)\n', ...
        i,alphas(i),alphaTStats(i),betas(i),betaTStats(i));
    
    clear tmp*;

end;

clear x;
    
subplot(2,2,1);
bar(alphas)
title('Estimated industry alphas');
subplot(2,2,2);
bar(alphaTStats)
title('Estimated industry alpha t-stats');
subplot(2,2,3);
bar(betas)
title('Estimated industry betas');
subplot(2,2,4);
bar(betaTStats)
title('Estimated industry beta t-stats');


figure;
plot(betas,xsrMean,'k+')
title('Estimated Security Market Line');


% Construct J0 and J1

sigma = cov(residuals);

fprintf('\nWald tests of all alphas=0:\n');

J0 = T*alphas*inv(sigma)*alphas'/(1+xsmktMean^2/xsmktStd^2);
fprintf('J0 = %5.3f, 5 pct critical value = %5.3f\n',J0,chi2inv(0.95,K));

J1 = (T-K-1)*alphas*inv(sigma)*alphas'/(K*(1+xsmktMean^2/xsmktStd^2));
fprintf('J1 = %5.3f, 5 pct critical value = %5.3f\n',J1,finv(0.95,K,T-K-1));