*
* Example 11.5 from p 235
*
open data tablef9-1[1].txt
data(format=prn,org=columns) 1 100 mdr acc age income avgexp ownrent selfempl
*
set incomesq = income**2
set posexp   = avgexp>0
*
smpl(series=posexp)
linreg avgexp
# constant age ownrent income incomesq
*
set logusq = log(%resids**2)
linreg logusq
# constant income incomesq
compute gamma=%beta,gamma(1)=gamma(1)+1.2704
prj fitted
set spread = exp(fitted)
*
* This does up to 100 iterations, but breaks the loop if the squared length of
* the adjustment to the scedastic function is small.
*
do iters=1,100
   linreg(spread=spread,noprint) avgexp
   # constant age ownrent income incomesq
   set usqgap = (%resids**2/spread)-1
   linreg(noprint) usqgap
   # constant income incomesq
   compute gamma=gamma+%beta
   if %normsqr(%beta)<=1.e-12
      break
   prj adjust
   set spread = spread*exp(adjust)
end do iters
disp iters "Iterations taken"
*
* The estimates here will depend upon the criterion used to exit the iteration
* loop.
*
linreg(spread=spread) avgexp
# constant age ownrent income incomesq
*
* Compute the covariance matrix of the parameters in the scedastic function. The
* "%xx" matrix from the auxiliary regression is multiplied by 2.
*
linreg(noprint) usqgap
# constant income incomesq
linreg(create,lastreg,title="Estimates of Scedastic Function",coeffs=gamma,covmat=2.0*%xx,noscale) usqgap