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)))]);

