options crt;  ? 3e5
in nerlov;
smpl 1,145;
lncp3 = log(costs/pf);
lnp13 = log(pl/pf);
lnp23 = log(pk/pf);
lnkwh = log(kwh);
group = int(order/100);  ? 1,2,3,4,5 for groups
dummy group;  ? creates dummies group1-group5
frml eqr r = 1/lnkwh;
set ssrtot = 0;
dot 1-5;
  select group.;
  olsq lncp3 c lnkwh lnp13 lnp23;  ? (a)
   set ssrtot = ssrtot + @ssr;
  analyz eqr;  ? (b)
  select 1;
  lny. = group. * lnkwh;  ? makes lny1-lny5 for below
  frml eqr. r. = 1/lny.;    ? makes eqr1-eqr5 for below
enddot;
olsq lncp3 lny1-lny5 lnp13 lnp23 group1-group5;  ? (c)
analyz eqr1-eqr5;  ? (d)
set df1 = 2*4;
set df2 = @nob-4*5;
set fstat = (@ssr-ssrtot)/df1/(ssrtot/df2);
cdf(f,df1=df1,df2=df2) fstat;  ? (e)
? Note:  the coefficient of lny2 *does* depend on the base of logarithms
? used, because the multiplicative factor log(10) enters it twice, so
? this factor no longer cancels from all but the intercept.
lny2 = lnkwh*lnkwh;
olsq lncp3 c lnkwh lny2 lnp13 lnp23;
  rename @ssr ssrunr;
  unmake @coef b0 by byy b1 b2;  param by byy;  ? save results for later
  copy @vcov vcvyy;
?
?  There are 2 different ways to make the F-test for whether
?  lnkwh and lny2 belong in the above equation.
? 1. The easiest way is to use the F version of ANALYZ in TSP 4.4 and later:
  frml r1 lnkwh;
  frml r2 lny2;
  analyz r1 r2;
  copy @fst feasy;
?
? 2. A harder way is to run the restricted regression and form the
?    F statistic by hand:
olsq(silent) lncp3 c lnp13 lnp23;
  rename @ssr ssrest;
set nob5 = @nob-5;
set fstat = (ssrest-ssrunr)/2/(ssrunr/(@nob-5));
cdf(f,df1=2,df2=nob5) fstat;
print feasy fstat;

? Compute returns to scale at median observations.
? Since both lnkwh and its regression coefficient enter the frml,
? we have to save the coefficient under a different name (by).

frml eqryy r = 1/(by+2*byy*lnkwh);
set jmed = 15;
do i=1,5;
  smpl jmed,jmed;  ? select median obs 15, 15+29, etc.
  analyz(names=(b0 by byy b1 b2),vcov=vcvyy) eqryy;    ? (f)
  set jmed=jmed+29;
enddo;