
[TC, Q, P1, P3, P2] = textread('nerlove.txt','%f %f %f %f %f');
ON=ones(size(TC));
%part b--------------------------------------------------
disp('b');
X=[ON log(Q) log(P1) log(P2) log(P3)];
Y=[log(TC)];

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

%part c----------------------------------------------
disp('c');
X=[ON log(Q) log(P1./P3) log(P2./P3)];
Y=[log(TC./P3)];

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))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c          Q       P1/P3     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 1 0 0];
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)))]);

%part d-----------------------------------------------
disp('d');
b=0; ster=0;
for i=1:5
X=[ON log(Q) log(P1./P3) log(P2./P3)];
Y=[log(TC./P3)];
k=i*29; m=i*29-28;
X=X(m:k,:);
Y=Y(m:k,:);

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))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c          Q       P1/P3     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'*X(:,1))/(X(:,1)'*X(:,1)))'*(Y-(Y'*X(:,1))/(X(:,1)'*X(:,1))));
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 1 0 0];
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)))]);
disp(['SSR        ' num2str((Y-X*beta)'*(Y-X*beta))]);
b(1)=b(1)+(Y-X*beta)'*(Y-X*beta);
b=[b beta(2)];  ster=[ster (((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))];
end
disp(['sum of SSR        ' num2str(b(1))]);
b=b(2:6);   ster=ster(2:6);     rts=1./b;
disp('Summary: beta2, Error variance, returns to scale');
disp([b; ster; rts]);

%part e-----------------------------------------------------
disp('e');
X=[ON log(Q) log(P1./P3) log(P2./P3)];
Y=[log(TC./P3)];
X=[            X(1:29,:) zeros(29,16)
   zeros(29,4) X(30:58,:) zeros(29,12)
   zeros(29,8) X(59:87,:) zeros(29,8)
   zeros(29,12) X(88:116,:) zeros(29,4)
   zeros(29,16) X(117:145,:) 
];

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))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c          Q       P1/P3     P2/P3           X 5');
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'*X(:,1))/(X(:,1)'*X(:,1)))'*(Y-(Y'*X(:,1))/(X(:,1)'*X(:,1))));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);
disp(['SSR         ' num2str((Y-X*beta)'*(Y-X*beta))]);

%part f---------------------------------------------------
disp('f')
GE0=[0 0 0 0]; GE1=[1 0 0 0];  GE2=[-1 0 0 0];  
GE3=[GE1 GE2 GE0 GE0 GE0
    GE1 GE0 GE2 GE0 GE0
    GE1 GE0 GE0 GE2 GE0
    GE1 GE0 GE0 GE0 GE2];
GE= [ GE3(:,1:20);    
    zeros(4,1) GE3(:,1:19);     
    zeros(4,2) GE3(:,1:18);     
    zeros(4,3) GE3(:,1:17)];
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)))]);

%part g-----------------------------------------------------
disp('g');
GE= [ zeros(4,2) GE3(:,1:18);     
      zeros(4,3) GE3(:,1:17)];
F=(GE*beta)'*inv(GE*V*GE')*(GE*beta)/size(GE,1);
disp(['F-stat of hypothesis: coefs 3 and 4 dont differ across groups ' num2str(F)]);
disp(['P-value of hypothesis: coefs 3 and 4 dont differ across groups ' num2str(1-Fcdf(F,size(GE,1),size(X,1)-size(X,2)))]);

X=[ON log(Q) log(P1./P3) log(P2./P3)];
Y=[log(TC./P3)];
X=[            X(1:29,1:2) zeros(29,8)   X(1:29,3:4)
   zeros(29,2) X(30:58,1:2) zeros(29,6)  X(30:58,3:4)
   zeros(29,4) X(59:87,1:2) zeros(29,4)  X(59:87,3:4)
   zeros(29,6) X(88:116,1:2) zeros(29,2) X(88:116,3:4)
   zeros(29,8) X(117:145,:)
];

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))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c          Q     X 5             P1/P3     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'*X(:,1))/(X(:,1)'*X(:,1)))'*(Y-(Y'*X(:,1))/(X(:,1)'*X(:,1))));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));
disp(R_2);

%part e-----------------------------------------------------
disp('h');
X=[ON log(Q) log(Q).*log(Q) log(P1./P3) log(P2./P3)];
Y=[log(TC./P3)];

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))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c       log(Q)     log(Q)^2   P1/P3     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);

%addon
disp('correct for heteroscedasticity');
figure(1)
plot(Y-X*beta);
Z=((Y-X*beta).*(Y-X*beta));
figure(2)
plot(1./Q, Z);
X=[ON 1./Q];
Y=[Z];

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))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c          1/Q  ');
disp([beta ste t pv]');
ERR=X*beta;
H=inv(diag(ERR));
plot(1./Q,Z,1./Q, ERR);

%correction for heteroscedasticity
ER2=ERR.^0.5;
X=[ON./ER2 log(Q)./ER2 log(Q).*log(Q)./ER2 log(P1./P3)./ER2 log(P2./P3)./ER2];
Y=[log(TC./P3)./ER2];

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))));
t=GE*beta./ste;
pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c       log(Q)     log(Q)^2   P1/P3     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);
figure(3);
plot(Y-X*beta)
