%% Day4 - Econometrics
%
% 1) Check the empirical coverage of bootstrap test for equality of means
% in two samples of unequal sizes
% 2) We will look at the residual bootstrap to get CI for beta in a regression
%
%
%
%%
% modified example from 
% http://courses.washington.edu/matlab1/Bootstrap_examples.html
% 

% Bootstrapping instead of a t-test (with unequal sample sizes)
nReps = 1000;       %depends on your computer and your patience :)
n1 = 30;            %sample size 1
n2 = 15;            %sample size 2
alpha = .05;        %alpha value


x1 = randn(n1,1);
x2 = randn(n2,1);

%define the statistic as the difference between means
myStatistic = @(x1,x2) mean(x1)-mean(x2);

sampStat = myStatistic(x1,x2);
bootstrapStat = zeros(nReps,1);
for i=1:nReps
    sampX1 = x1(ceil(rand(n1,1)*n1));
    sampX2 = x2(ceil(rand(n2,1)*n2));
    bootstrapStat(i) = myStatistic(sampX1,sampX2);
end


CI = prctile(bootstrapStat,[100*alpha/2,100*(1-alpha/2)]);

%Hypothesis test: Does the confidence interval cover zero?
H = CI(1)>0 | CI(2)<0;

%Draw a histogram of the sampled statistic
clf
xx = min(bootstrapStat):.01:max(bootstrapStat);
hist(bootstrapStat,xx);
hold on
ylim = get(gca,'YLim');
h1=plot(sampStat*[1,1],ylim,'y-','LineWidth',2);
h2=plot(CI(1)*[1,1],ylim,'r-','LineWidth',2);
plot(CI(2)*[1,1],ylim,'r-','LineWidth',2);
h3=plot([0,0],ylim,'b-','LineWidth',2);
xlabel('Difference between means');

decision = {'Fail to reject H0','Reject H0'};
title(decision(H+1));
legend([h1,h2,h3],{'Sample mean',sprintf('95%% CI',100*alpha),'H0 mean'},'Location','NorthWest');


%%%%%%%%%%%% your task is described here %%%%%%%%%%%%
% try to modify this code to find out what is the
% coverage probability of the bootstrap CI
% So, re-run the code many times (1000? maybe 10000?) times. If the CI's coverage
% is correct, the CI should cover the H0 mean approximately
% in 95% cases (that is, 950 or 9500 times)
% Don't forget to comment out the histogram part
%%%%%%%%%%%% the end of your task %%%%%%%%%%%%%%%%%%%


%%
% Residual Bootstrap to get CI in regression
clear all
clc

rng('default')

% data generating process
beta_1  = 1;
beta_2  = 2;
nSample = 25;
x       = randn(nSample,1);
X       = [ones(nSample,1) x];
epsilon = randn(nSample,1);
y       = beta_1 + beta_2*x + epsilon;

resLeSage = ols(y,X);
bOLS    = resLeSage.beta;
bstd    = resLeSage.bstd;

% save residuals
resid   = y - X*bOLS;

% we want to test H0: beta_2 = 2
beta2H0 = 2;
% the t-statistic from our data is
tStat   = (bOLS(2)-2)/bstd(2);

% we want a test without normality assumption on errors
% for this, we will make use of the Bootstrap

%number of bootstrap replications
nBoot   = 2000; 

% pre-allocate tStatBoot - t-values for bootstrapped
tStatBoot = zeros(nBoot,1);

%%%%%%%%%%%% your task is described here %%%%%%%%%%%%
% draw nBoot samples of residuals
for iBoot = 1:nBoot
    % (i) Draw nSample random values from the vector of residuals
    % First, generate uniform random variables from 1,....,nSample
    % These uniform RVs will tell you which residual to pick
    % Save those to residBoot
    

    % (ii) Compute nSample y values according to the model of H0 hypothesis
    % multiply residuals by sqrt(N/(N-K)), N is sample size, K is # of
    % regressors (including constant)
    

    % (iii) Regress yBoot on X and get OLS estimator bOLSBoot and standard
    % errors of this estimator bstdBoot
    

    % (iv) Calculate t-statistic under H0 that beta_2 = 2
    
end
%%%%%%%%%%%% the end of your task %%%%%%%%%%%%%%%%%%%

% sort the bootstrapped t-statistic
sortedTStat = sort(tStatBoot);
alpha       = 0.05;
tStat
CIBoot      = [sortedTStat(round(nBoot*(alpha/2))), sortedTStat(round(nBoot*(1-alpha/2)))]
% The t-statistic from our data is within the bootstraped CI, so we do not
% reject H0: beta = 2

%%
%

