options crt;  ? 8e7
in cigad;
freq a;
smpl 30,78;
dot salespc incpc rprice astock;
  l. = log(.);
enddot;
d53 = (time >= 53);
d64 = (time >= 64);
smpl 30,70;
olsq lsalespc c lincpc lrprice lastock df d53 d64;  ? (a.1)
smpl 71,78;
forcst lsfit;
expfit = exp(lsfit);
salesf = exp(lsfit + @s2/2);  ? unbiased exponential forecast, based on
                              ? mean of a lognormal variable
smpl 78,78;
eoact = expfit/salespc;
uoact = salesf/salespc;
print salespc expfit salesf eoact uoact;  ? (a.2)

smpl 30,78;
d71 = (time >= 71);
olsq lsalespc c lincpc lrprice lastock df d53 d64 d71;  ? (b)

lcigtob = log(cigtob);
tovy = 2/incpc;
frml eq861 lcigtob = alpha + log(1 + yc*tovy) - yc*tovy;
param alpha,0 yc,1;  ? start yc > 0, to avoid saddlepoint
msd(noprint) lcigtob;  set alpha = @mean;  ? starting value for alpha
lsq eq861;  ? (c)

?corr lni lincpc;  ? same results repeated just below
lincpc2 = lincpc**2;
corr lni lincpc lincpc2;
olsq lsalespc c lincpc lincpc2 lrprice lastock df d53 d64 d71;
incelast = @coef(2) + 2*@coef(3)*lincpc;  ? derivative w.r.t. lincpc
print incelast;  ? (d) - print income elasticity for all values of lincpc
? compute a standard error for the income elasticity, at the mean
? value of lincpc, to make sure it's not too noisy
msd lincpc;
frml eqe incel = lincpc + 2*lincpc2*@mean;
analyz eqe;

ltpc = log(tpercig);
frml eq862 lsalespc = b0 + by*lincpc + bp*lrprice + bi*lni +
     ba*log(astock1 + r*astock2) + bdf*df + bf*f + bl*l + btpc*ltpc;
param b0,2.243 by,.462 bp,-1.218 bi,.971
      ba,.046 r,.264 bdf,-.075 bf,-.0021 bl,-.0235 btpc,-1.386;
lsq(step=golden) eq862;  ? (e.1) -- only 15 iterations
?lsq(maxit=200) eq862;  ? default step takes 117 iterations
const by,.462 ;  ? hold value of by fixed
lsq(step=golden) eq862;  ? (e.2) -- 8 iterations to LogL=54.0595
?lsq(maxit=50) eq862;  ? default step takes 39 iterations, to LogL=53.7732,
                       ? an inferior optimum

? (f) is rather open-ended; here are some of the suggestions
? from the book

?  Convert f from a percent (0-100) to a fraction (0-1)
f = f/100;
salesnf = salespc*(1-f);
lcnf = log(salesnf);
olsq lcnf c lincpc lincpc2 lrprice lastock df d53 d64 d71;

smpl 55,78;
salesf = salespc*f;
lcf = log(salesf);
olsq lcf c lincpc lincpc2 lrprice lastock df d64 d71;

logitf = log(f/(1-f));
wt = 1/(1-f);  ? this may not be the right specification for the
               ? weight in this type of regression...
olsq(weight=wt) logitf c lincpc lincpc2 lrprice lastock df d64 d71;