options crt;  ? 4e7
in waugh;
smpl 1,200;
lp = log(price);
lg = log(green);
ln = log(nostalks);
select disperse > 0;
ld = log(disperse);
msd(noprint) lp;
rename @sum sumlny;
olsq price c green nostalks disperse;
 set logl1 = @logl;
 unmake @coef b1-b4;  set sigi = 1/@s;
olsq lp c lg ln ld;
 set logl2 = @logl - sumlny;
olsq lp c green nostalks disperse;
 set logl3 = @logl - sumlny;
olsq price c lg ln ld;
 set logl4 = @logl;
print logl1-logl4;
?
?  Box-Tidwell regressions are done most simply with the FIML
?  command, which will use the proper Jacobian automatically.
?
param ly 1 lx 1 sigi b1-b4;
?
frml btres (price**ly-1)/ly - (b1 + b2*(green**lx-1)/lx +
                            b3*(nostalks**lx-1)/lx +
                            b4*(disperse**lx-1)/lx ) ;
?
?  FIML should yield the following results:
?  LogL = -784.015,  LY = .616,  LX = 1.25 .
?  The default HITER=G will actually work fine for this model, but HITER=B
?  also works OK, and HITER=B is preferred in part (b) of this exercise.
?  This FIML command takes about 12 seconds on a Pentium 166Mhz
fiml(endog=price,hiter=b,tol=1.e-7,maxit=100) btres;
 rename @logl loglml;
dot 1-4;
  set test. = 2*(loglml-logl.);
enddot;
print test1-test4;  ? LR tests, chisq(2)
?
frml eqlam diff = ly-lx;
analyz eqlam;  ? wald test
delete price lp;

smpl 1,91;
in cole;
lp = log(price);
ls = log(speed);
lc = log(cap);
msd(noprint) lp;
rename @sum sumlny;
olsq price speed cap c d73-d84;
 set logl1 = @logl;
 unmake @coef b1-b15;  set sigi = 1/@s;
olsq lp ls lc c d73-d84;
 set logl2 = @logl - sumlny;
olsq lp speed cap c d73-d84;
 set logl3 = @logl - sumlny;
olsq price ls lc c d73-d84;
 set logl4 = @logl;
print logl1-logl4;
?
?  FIML version
?
frml btres (price**ly-1)/ly - (b1*(speed**lx-1)/lx +
                       b2*(cap**lx-1)/lx +
                       b3+b4*d73+b5*d74+b6*d75+b7*d76+
                       b8*d77+b9*d78+b10*d79+b11*d80+
                       b12*d81+b13*d82+b14*d83+b15*d84 );
param ly 1 lx 1 sigi b1-b15;
set b3 = b3 - 1/ly + (b1+b2)/lx;  ? adjust constant term for -1/lx stuff
?
?  FIML should yield the following results:
?  LogL = -919.131 (says failure to improve, but is actually converged)
?  ly = .414, lx = .021
?  This FIML command takes about 20 seconds on a Pentium 166Mhz
? For this particular model/dataset,
? it is important to *not* use the default HITER=G for FIML -- this leads
? to false convergence, with LogL = -924, ly = .9 or so.
fiml(endog=price,tol=1.e-7,maxit=300,hiter=b) btres;
 rename @logl loglml;
dot 1-4;
  set test. = 2*(loglml-logl.);
enddot;
print test1-test4;  ? LR tests, chisq(2)
?
frml eqlam diff = ly-lx;
analyz eqlam;  ? wald test