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