clear all; close all; clc; [Y1 Y2 Y3 X1 X2 X3 X4 X5 X6 X7 X8 X9 X10] = textread('hw3p2.txt','%f %f %f %f %f %f %f %f %f %f %f %f %f'); %ON=ones(size(Y1)); %part b-------------------------------------------------- disp('1,3. OLS and SE'); X=[X1 X2 X3 X4]; Y=Y1; beta=inv(X'*X)*(X'*Y); V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X); GE=eye(size(X,2)); ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2)))); beta1=beta; err1=Y-X*beta; ste1=ste; X=[X5 X6 X7 X8]; Y=Y2; beta=inv(X'*X)*(X'*Y); V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X); GE=eye(size(X,2)); ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2)))); beta2=beta; err2=Y-X*beta; ste2=ste; X=[X1 X2 X5 X6 X9 X10]; Y=Y2; beta=inv(X'*X)*(X'*Y); V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X); GE=eye(size(X,2)); ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2)))); beta3=beta; err3=Y-X*beta; ste3=ste; disp(' X1 X2 X3 X4'); disp(['Y1 ' num2str(beta1')]); disp(['se ' num2str(ste1')]); disp(' X5 X6 X7 X8'); disp(['Y2 ' num2str(beta2')]); disp(['se ' num2str(ste2')]); disp(' X1 X2 X5 X6 X9 X10'); disp(['Y3 ' num2str(beta3') ]); disp(['se ' num2str(ste3')]); %sigma estimate E=[err1 err2 err3]; disp('2. Matrix sigma'); SIG=E'*E/size(E,1); disp(SIG); N=size(X1,1); %SUR estimates X=[X1, X2, X3, X4, zeros(N,4+6); zeros(N,4), X5, X6, X7, X8, zeros(N,6); zeros(N,4+4), X1, X2, X5, X6, X9, X10]; Y=[Y1; Y2; Y3]; iOm=kron(inv(SIG),eye(N)); iXOX=inv(X'*iOm*X); beta=iXOX*(X'*iOm*Y); V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*iOm*X); GE=eye(size(X,2)); ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2)))); beta0=[beta1; beta2; beta3]; ste0=[ste1; ste2; ste3]; %test hypothesis GE=[1 zeros(1,7) -1 zeros(1,5); 0 1 zeros(1,7) -1 zeros(1,4); zeros(1,4) 1 zeros(1,5) -1 zeros(1,3) zeros(1,5) 1 zeros(1,5) -1 zeros(1,2) ]; F=(GE*beta)'*inv(GE*V*GE')*(GE*beta)/size(GE,1); disp(['F-stat of hypothesis: coefs dont differ across groups ' num2str(F)]); disp(['P-value of hypothesis: coefs dont differ across groups ' num2str(1-Fcdf(F,size(GE,1),size(X,1)-size(X,2)))]); beta2=beta-iXOX*GE'*inv(GE*iXOX*GE')*GE*beta; disp(' steOLS betaOLS betaGLS steGLS betaGLSR'); disp([ste0 beta0 beta ste beta2]); % t=GE*beta./ste; % pv=2-2*tcdf(abs(t),size(X,1)-size(X,2)); % disp(' c Q P1 P2 P3 '); % disp([beta ste t pv]'); % disp('R-squared'); % R_2=zeros(1,2); % R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON))); % R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2)); % disp(R_2); % % GE=[0 0 1 1 1]; % t=(GE*beta-1)./sqrt(spdiags((GE*V*GE'))); % disp(['P-value of hypothesis: CRS ' num2str(2-2*tcdf(abs(t),size(X,1)-size(X,2)))]);