n = 100;
B = 1000;
X = rand(n,1);
epsilon = rand(n,1);
epsilon = epsilon - mean(epsilon);
beta = [-10]';
Y = X*beta + epsilon;
[b,bint,r,rint,stats] = regress(Y,X);
disp('[beta,betaint]')
[b,bint]
disp('Std(beta)')
Std = sqrt(var(epsilon)*diag(inv(X'*X)))
disp('Bootstrap')
Z = bootstrp(B,'regress',Y,X);
disp('mean(Z)')
mean(Z)
disp('std(Z)')
std(Z)
hist(Z)