clear; close; clc;
%Define the population parameters
alpha=0.5;
beta=0.9;
% Sample size n
n=10000;
% Number of samples m
m=1000;
% Store the estimated beta in a vector
beta_hat=zeros(m,1);
%by hand
for i=1:m
    %Generate independent variable randomly
    x= randn(n,1);
    %Generate errors in the population
    e=randn(n,1);
    %Generate the dependent variable
    y=alpha+beta*x+e;
    %Generate the LS estimates of alpha and beta using matrix formulation
    X=[ones(n,1) x];
    beta_hatvec=(inv((X'*X)))*X'*y;
    %define the residuals
    resid = y-(X*beta_hatvec);
    %estimate of sigma_2
    sigma2_hat=(resid'*resid) / (size(X,1)-size(X,2));
    %estimate of Vhat
    vcov_beta_hat = (sigma2_hat*inv((X.'*X)));
    stderror_hatvec = sqrt(diag(vcov_beta_hat));
    % Pull out only the second component - the estimate of beta
    beta_hat(i)=beta_hatvec(2);
    std_error(i)=stderror_hatvec(2);
end
        t_stat=(beta_hat)./std_error.';
        sort_t_stat=sort(t_stat);

        crit10_actual=sort_t_stat(900);
        crit5_actual=sort_t_stat(950);
        crit1_actual=sort_t_stat(995);

        critical = [crit1_actual crit5_actual crit10_actual]