%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Xavier Martin G. Bautista
% Econometrics 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The following estimates the coefficients of a regression using OLS and
% the standard errors by bootstrap.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

% Set the true parameters and placeholders for results.

T = 100;                                                                    % Sample size.

beta0 = 0.5;                                                                % Intercept.
beta1 = 3;                                                                  % Coefficient on X1.
beta2 = 1;                                                                  % Coefficient on X2.

sigma = 1;                                                                  % Standard deviation.

pdf = 0;                                                                    % Distribution of error terms. 0 = normal distribution, 1 = cauchy distribution.

bootsim = 1500;                                                             % Number of bootstrap repetitions.
sim = 200;                                                                  % Number of simulations.

b_ols = zeros(sim,3);                                                       % OLS estimates.
se_ols = zeros(sim,3);                                                      % OLS standard errors.                
se_boot = zeros(sim,3);                                                     % Bootstrap standard errors.

% Generate X data.

x1 = 2 + ((1:1:T)'/T).*normrnd(0,1,T,1);
x2 = 3 + 0.5*x1 + sigma.*normrnd(0,1,T,1);

x = [ones(T,1) x1 x2];

for s = 1:sim
    
    % Generate y data and u.
    
    if pdf == 0    
        u = sigma.*normrnd(0,1,T,1);                                        % Normal errors.
    elseif pdf == 1
        u = (sigma.*randn(n,1))./(sigma.*randn(n,1));                       % Cauchy errors.
    end
    
    y = beta0 + beta1.*x1 + beta2.*x2 + u;
    
    % OLS estimates.

    b_ols(s,:) = (inv(x'*x)*x'*y)';                                         % OLS estimates.

    % OLS standard errors.

    uhat = y - x*b_ols(s,:)';                                               % Residuals.
    s2 = (uhat'*uhat)/(T - 3);                                              % s squared.

    se_ols(s,:) = sqrt(s2*diag(inv(x'*x)))';                                % OLS standard errors.
    
    % Bootstrap standard errors.
    
    b_boot = zeros(bootsim,3);
    
    for i = 1:bootsim
        
        yb = x*b_ols(s,:)';                                                 % yhat from original estimation.
    
        yboot = zeros(T,1);
        e = zeros(T,1);
            
        % Draw the error terms from the empirical distribution and
        % construct the y's for the bootstrap.
        
        draw = ones(T,1) + T*rand(T,1);    
        draw = fix(draw);
    
        for j = 1:size(x,1)
            e(j,1) = uhat(draw(j,1),1);                                     % Resampled errors.
            yboot(j,1) = yb(j,1) + e(j,1);                                  % Bootstrap y's.
        end
    
        % Use OLS to compute the estimated coefficients with yprime.
    
        b_boot(i,:) = (inv(x'*x)*x'*yboot)';                                % OLS estimates from each bootstrap repetition.
        
    end

    bootse = (std(b_boot,0,1))';                             
    
    se_boot(s,1:size(se_boot,2)) = bootse';                                 % Bootstrap standard errors.

end

% Display the results of the last estimation.

fprintf('--------------------------------------------------------------\n')
fprintf('Results of Last Simulation \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('Sample Size: %d \n', T)
fprintf('Bootstrap Repetitions: %d \n', bootsim)
fprintf('Number of Simulations: %d \n', sim)
fprintf('\n')

fprintf('OLS Estimates \n')
fprintf('\n')
fprintf('            b           SE(OLS)      SE(Boot)          \n')
fprintf('b0        %0.4f        %0.4f        %0.4f       \n', b_ols(sim,1), se_ols(sim,1), se_boot(sim,1))
fprintf('b1        %0.4f        %0.4f        %0.4f       \n', b_ols(sim,2), se_ols(sim,2), se_boot(sim,2))
fprintf('b2        %0.4f        %0.4f        %0.4f       \n', b_ols(sim,3), se_ols(sim,3), se_boot(sim,3))
fprintf('\n')

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('Sample Size: %d \n', T)
fprintf('Bootstrap Repetitions: %d \n', bootsim)
fprintf('Number of Simulations: %d \n', sim)
fprintf('\n')

fprintf('Empirical Means and Standard Deviations \n')
fprintf('\n')
fprintf('Coefficients \n')
fprintf('          Means      Std. Devs      \n')
fprintf('b0        %0.4f        %0.4f        \n', mean(b_ols(:,1)), std(b_ols(:,1),0,1))
fprintf('b1        %0.4f        %0.4f        \n', mean(b_ols(:,2)), std(b_ols(:,2),0,1))
fprintf('b2        %0.4f        %0.4f        \n', mean(b_ols(:,3)), std(b_ols(:,3),0,1))
fprintf('\n')
fprintf('OLS Standard Errors \n')
fprintf('          Means      Std. Devs      \n')
fprintf('SE(b0)    %0.4f        %0.4f        \n', mean(se_ols(:,1)), std(se_ols(:,1),0,1))
fprintf('SE(b1)    %0.4f        %0.4f        \n', mean(se_ols(:,2)), std(se_ols(:,2),0,1))
fprintf('SE(b2)    %0.4f        %0.4f        \n', mean(se_ols(:,3)), std(se_ols(:,3),0,1))
fprintf('\n')
fprintf('Bootstrap Standard Errors \n')
fprintf('          Means      Std. Devs      \n')
fprintf('SE(b0)    %0.4f        %0.4f        \n', mean(se_boot(:,1)), std(se_boot(:,1),0,1))
fprintf('SE(b1)    %0.4f        %0.4f        \n', mean(se_boot(:,2)), std(se_boot(:,2),0,1))
fprintf('SE(b2)    %0.4f        %0.4f        \n', mean(se_boot(:,3)), std(se_boot(:,3),0,1))
fprintf('\n')
