clear;
%part a--------------------------------------------------
hold on;
disp('a');
S=1000;     beta1=zeros(1,S);    ste1=zeros(1,S);     t1=zeros(1,S);    
            beta2=zeros(1,S);    ste2=zeros(1,S);     t2=zeros(1,S);    
            beta3=zeros(1,S);    ste3=zeros(1,S);     t3=zeros(1,S);    
            beta4=zeros(1,S);    ste4=zeros(1,S);     t4=zeros(1,S);    

for j=1:9

if j<4
    sigma=0.5; 
elseif j<7
    sigma=1; 
else
    sigma=2; 
end    
if j==1 | j==4 | j==7
    n=100; 
elseif j==2 | j==5 | j==8
    n=400; 
else
    n=1600; 
end    

out15=0; out110=0;  out25=0; out210=0;  out35=0; out310=0;  out45=0; out410=0;  
for i=1:S
    X=randn(n,1);     e=randn(n,1);         Y=X+e*sigma;    
%OLS
    beta1(i)=inv(X'*X)*(X'*Y);
    V=(((Y-X*beta1(i))'*(Y-X*beta1(i)))/(size(X,1)-size(X,2)))*inv(X'*X);
    ste1(i)=sqrt(V);             t1(i)=(beta1(i)-1)/ste1(i);
    pv=2-2*tcdf(abs(t1(i)),size(X,1)-size(X,2));
    if pv>0.975 | pv<0.025
        out15=out15+1;
    end
    if pv>0.95 | pv<0.05
        out110=out110+1;
    end
%White errors
    beta2=beta1;
    er=Y-X*beta2(i);
    V=inv(X'*X)*(X'*(er.*er.*X))*inv(X'*X);
    ste2(i)=sqrt(V);             t2(i)=(beta2(i)-1)/ste2(i);
    pv=2-2*tcdf(abs(t2(i)),size(X,1)-size(X,2));
    if pv>0.975 | pv<0.025
        out25=out25+1;
    end
    if pv>0.95 | pv<0.05
        out210=out210+1;
    end
%GLS=OLS
    X3=X;       Y3=Y;   
    beta3(i)=inv(X3'*X3)*(X3'*Y3);
    V=(((Y3-X3*beta3(i))'*(Y3-X3*beta3(i)))/(size(X3,1)-size(X3,2)))*inv(X3'*X3);
    ste3(i)=sqrt(V);             t3(i)=(beta3(i)-1)/ste3(i);
    pv=2-2*tcdf(abs(t3(i)),size(X3,1)-size(X3,2));
    if pv>0.975 | pv<0.025
        out35=out35+1;
    end
    if pv>0.95 | pv<0.05
        out310=out310+1;
    end
%FGLS
    er=(Y-X*beta1(i)).*(Y-X*beta1(i));
    G=[ones(size(X)) X];
    gamma=inv(G'*G)*(G'*er);
    erf=sqrt(G*gamma);

    X4=X./erf;  Y4=Y./erf;
    beta4(i)=inv(X4'*X4)*(X4'*Y4);
    V=(((Y4-X4*beta4(i))'*(Y4-X4*beta4(i)))/(size(X4,1)-size(X4,2)))*inv(X4'*X4);
    ste4(i)=sqrt(V);             t4(i)=(beta4(i)-1)/ste4(i);
    pv=2-2*tcdf(abs(t4(i)),size(X4,1)-size(X4,2));
    if pv>0.975 | pv<0.025
        out45=out45+1;
    end
    if pv>0.95 | pv<0.05
        out410=out410+1;
    end
end
figure(1)
grid=0.1:0.025:1.9;  beta1=sort(beta1);      beta2=sort(beta2);
beta3=sort(beta3);      beta4=sort(beta4);
subplot(3,6,2*j-1)
hist(beta3,grid)
title(['GLS ' num2str([n,sigma])])
subplot(3,6,2*j)
hist(beta4,grid)
title(['FGLS ' num2str([n,sigma])])
%mean bias
mb1=ones(size(beta1))*beta1'/size(beta1,2)-1;
mb2=ones(size(beta2))*beta2'/size(beta2,2)-1;
mb3=ones(size(beta3))*beta3'/size(beta3,2)-1;
mb4=ones(size(beta4))*beta4'/size(beta4,2)-1;
%median bias
meb1=(beta1(S/2)+beta1(1+S/2))/2-1;
meb2=(beta2(S/2)+beta2(1+S/2))/2-1;
meb3=(beta3(S/2)+beta3(1+S/2))/2-1;
meb4=(beta4(S/2)+beta4(1+S/2))/2-1;
%root mean square error
rmse1=sqrt(ones(size(beta1))*(((beta1-1).*(beta1-1))')/S);
rmse2=sqrt(ones(size(beta2))*(((beta2-1).*(beta2-1))')/S);
rmse3=sqrt(ones(size(beta3))*(((beta3-1).*(beta3-1))')/S);
rmse4=sqrt(ones(size(beta4))*(((beta4-1).*(beta4-1))')/S);
%median absolute error
ab=sort(abs(beta1-1));
medae1=(ab(S/2)+ab(1+S/2))/2;
ab=sort(abs(beta2-1));
medae2=(ab(S/2)+ab(1+S/2))/2;
ab=sort(abs(beta3-1));
medae3=(ab(S/2)+ab(1+S/2))/2;
ab=sort(abs(beta4-1));
medae4=(ab(S/2)+ab(1+S/2))/2;
disp(['    n/10       sigma      Fr5%    Frac10% Mean bias Median bias  RMSE      MedAE']);
disp([n/10 sigma out15/S out110/S mb1 meb1 rmse1 medae1; 
      n/10 sigma out25/S out210/S mb2 meb2 rmse2 medae2; 
      n/10 sigma out35/S out310/S mb3 meb3 rmse3 medae3; 
      n/10 sigma out45/S out410/S mb4 meb4 rmse4 medae4]);
end
hold off; print -dbmp a.bmp
clear;
%part b--------------------------------------------------
hold on;
disp('b');
S=1000;     beta1=zeros(1,S);    ste1=zeros(1,S);     t1=zeros(1,S);    
            beta2=zeros(1,S);    ste2=zeros(1,S);     t2=zeros(1,S);    
            beta3=zeros(1,S);    ste3=zeros(1,S);     t3=zeros(1,S);    
            beta4=zeros(1,S);    ste4=zeros(1,S);     t4=zeros(1,S);    
            
for j=1:9

if j<4
    c=1; 
elseif j<7
    c=2; 
else
    c=4; 
end    
if j==1 | j==4 | j==7
    n=100; 
elseif j==2 | j==5 | j==8
    n=400; 
else
    n=1600; 
end    

out15=0; out110=0;  out25=0; out210=0;  out35=0; out310=0;  out45=0; out410=0;  
for i=1:S
    X=randn(n,1);     e=randn(n,1);         Y=X+e.*sqrt((X.*X)/c);    
%OLS
    beta1(i)=inv(X'*X)*(X'*Y);
    V=(((Y-X*beta1(i))'*(Y-X*beta1(i)))/(size(X,1)-size(X,2)))*inv(X'*X);
    ste1(i)=sqrt(V);             t1(i)=(beta1(i)-1)/ste1(i);
    pv=2-2*tcdf(abs(t1(i)),size(X,1)-size(X,2));
    if pv>0.975 | pv<0.025
        out15=out15+1;
    end
    if pv>0.95 | pv<0.05
        out110=out110+1;
    end
%White errors
    beta2=beta1;
    er=Y-X*beta2(i);
    V=inv(X'*X)*(X'*(er.*er.*X))*inv(X'*X);
    ste2(i)=sqrt(V);             t2(i)=(beta2(i)-1)/ste2(i);
    pv=2-2*tcdf(abs(t2(i)),size(X,1)-size(X,2));
    if pv>0.975 | pv<0.025
        out25=out25+1;
    end
    if pv>0.95 | pv<0.05
        out210=out210+1;
    end
%GLS
    X3=X./sqrt(X.*X);       Y3=Y./sqrt(X.*X);   
    beta3(i)=inv(X3'*X3)*(X3'*Y3);
    V=(((Y3-X3*beta3(i))'*(Y3-X3*beta3(i)))/(size(X3,1)-size(X3,2)))*inv(X3'*X3);
    ste3(i)=sqrt(V);             t3(i)=(beta3(i)-1)/ste3(i);
    pv=2-2*tcdf(abs(t3(i)),size(X3,1)-size(X3,2));
    if pv>0.975 | pv<0.025
        out35=out35+1;
    end
    if pv>0.95 | pv<0.05
        out310=out310+1;
    end
%FGLS
    er=(Y-X*beta1(i)).*(Y-X*beta1(i));
    G=[ones(size(X)) X];
    gamma=inv(G'*G)*(G'*er);
    erf=sqrt(G*gamma);

    X4=X./erf;  Y4=Y./erf;
    beta4(i)=inv(X4'*X4)*(X4'*Y4);
    V=(((Y4-X4*beta4(i))'*(Y4-X4*beta4(i)))/(size(X4,1)-size(X4,2)))*inv(X4'*X4);
    ste4(i)=sqrt(V);             t4(i)=(beta4(i)-1)/ste4(i);
    pv=2-2*tcdf(abs(t4(i)),size(X4,1)-size(X4,2));
    if pv>0.975 | pv<0.025
        out45=out45+1;
    end
    if pv>0.95 | pv<0.05
        out410=out410+1;
    end
end
figure(2)
grid=0.1:0.025:1.9;  beta1=sort(beta1);      beta2=sort(beta2);
beta3=sort(beta3);      beta4=sort(beta4);
subplot(3,6,2*j-1)
hist(beta3,grid)
title(['GLS ' num2str([n,c])])
subplot(3,6,2*j)
hist(beta4,grid)
title(['FGLS ' num2str([n,c])])
%mean bias
mb1=ones(size(beta1))*beta1'/size(beta1,2)-1;
mb2=ones(size(beta2))*beta2'/size(beta2,2)-1;
mb3=ones(size(beta3))*beta3'/size(beta3,2)-1;
mb4=ones(size(beta4))*beta4'/size(beta4,2)-1;
%median bias
meb1=(beta1(S/2)+beta1(1+S/2))/2-1;
meb2=(beta2(S/2)+beta2(1+S/2))/2-1;
meb3=(beta3(S/2)+beta3(1+S/2))/2-1;
meb4=(beta4(S/2)+beta4(1+S/2))/2-1;
%root mean square error
rmse1=sqrt(ones(size(beta1))*(((beta1-1).*(beta1-1))')/S);
rmse2=sqrt(ones(size(beta2))*(((beta2-1).*(beta2-1))')/S);
rmse3=sqrt(ones(size(beta3))*(((beta3-1).*(beta3-1))')/S);
rmse4=sqrt(ones(size(beta4))*(((beta4-1).*(beta4-1))')/S);
%median absolute error
ab=sort(abs(beta1-1));
medae1=(ab(S/2)+ab(1+S/2))/2;
ab=sort(abs(beta2-1));
medae2=(ab(S/2)+ab(1+S/2))/2;
ab=sort(abs(beta3-1));
medae3=(ab(S/2)+ab(1+S/2))/2;
ab=sort(abs(beta4-1));
medae4=(ab(S/2)+ab(1+S/2))/2;

disp(['    n/10         c        Fr5%    Frac10% Mean bias Median bias  RMSE      MedAE    ']);
disp([n/10 c out15/S out110/S mb1 meb1 rmse1 medae1; 
      n/10 c out25/S out210/S mb2 meb2 rmse2 medae2; 
      n/10 c out35/S out310/S mb3 meb3 rmse3 medae3; 
      n/10 c out45/S out410/S mb4 meb4 rmse4 medae4]);
end
hold off; print -dbmp b.bmp