[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)