clear;
tri=1000;
n=500;
range=-1:0.1:1;
power=zeros(length(range),1);
power2=zeros(length(range),1);

for j=1:length(range);
    reject=0;
    reject2=0;
for i=1:tri
    x=randn(n,1);
    u= (x.^2).*randn(n,1);
    y=1+range(j)*x+u;
    X=[ones(n,1) x];
    beta=X\y;
    e=y-X*beta;
    var=inv(X'*X)*(e'*e)/(n-2); %wrong OLS var matrix
    var2=(inv(X'*X)*X'*diag(e.^2)*X*inv(X'*X)); %White var matrix
    t=(beta(2))/(var(2,2)^0.5);
    if abs(t)>1.96;
        reject=reject+1;
    end
        t2=(beta(2))/(var2(2,2)^0.5);
    if abs(t2)>1.96;
        reject2=reject2+1;
    end
    
    
end
    power(j,:)=reject/tri;
    power2(j,:)=reject2/tri;
end

plot(range, power, range, power2)
legend('t Test without Correction','t Test with White Correction');
title(['Why Heteroskedasticity Sucks']);
xlabel('Degree of Heteroskedasticity')
ylabel('Rejection Rate')