capm_with_constant;% The variables defined in the capm_with_constant.m are in matlab's workspace, 
                   % so the excess returns data, BETA, ALPH, etc. are known here.

% Set the number of observations.
n = length(Rfood);

% Compute the estimated residuals for each capm regression model.
capm_error_food = Rfood-BETA(1)*Rmrf-ALPH(1);
capm_error_durables = Rdur-BETA(2)*Rmrf-ALPH(2);
capm_error_construction = Rcon-BETA(3)*Rmrf-ALPH(3);

% Estimate the standard deviation of the errors for each capm regression model.
capm_error_food_std = sqrt(capm_error_food'*capm_error_food/(n-2));
capm_error_durables_std = sqrt(capm_error_durables'*capm_error_durables/(n-2));
capm_error_construction_std = sqrt(capm_error_construction'*capm_error_construction/(n-2));

capm_errors_std = [ capm_error_food_std ; capm_error_durables_std ; capm_error_construction_std ];

% Compute confidence bands for the constants (monte carlo)
B = 1000000; % Number of simulations
s = .950; % Probability mass covered by the confidence bands.
C = NaN(B,3); % Create a big matrix where all the estimated constants will be saved.


for b=1:B
    for model = 1:3
        % I am god and I create data under the null hypothesis (capm is true, the constant is zero)
        e = randn(n,1)*capm_errors_std(model);
        excess_returns = BETA(i)*Rmrf + e;
        % I am an econometrician and I estimate the capm regression model (with constant).
        params = regression(excess_returns,[ones(n,1), Rmrf]);
        C(b,model) = params(1);
    end
end

for model=1:3
    figure(model)
    histogram(C(:,model),1000,0.95);
end