*
* Example 11.4 from p 231
* Feasible GLS for heteroscedasticity
*
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 / resids
# constant age ownrent income incomesq
linreg(spread=income) avgexp
# constant age ownrent income incomesq
linreg(spread=incomesq) avgexp
# constant age ownrent income incomesq
*
* This is mislabeled in the 1st printing of the text. The explanatory variables
* for the log variance model are 1, I and I**2.
*
set logusq = log(resids**2)
linreg logusq
# constant income incomesq
prj logvar
set spread = exp(logvar)
linreg(spread=spread) avgexp
# constant age ownrent income incomesq
*
* Two step, power of income
*
set logi   = log(income)
linreg logusq
# constant logi
compute alpha=%beta(2)
prj logvar
set spread = exp(logvar)
linreg(spread=spread) avgexp / resids
# constant age ownrent income incomesq
*
* Iterated. We'll allow up to forty iterations. Each time, we'll check to see if
* the new value of alpha is very close to the last one. Once the (absolute value)
* of the gap is small enough, we break the loop. Inside the loop, we suppress the
* regression output. We show the result with the last value of alpha once we exit
* the loop.
*
do iters=1,40
   set logusq = log(resids**2)
   linreg(noprint) logusq
   # constant logi
   if abs(%beta(2)-alpha)<.00001
      break
   compute alpha=%beta(2)
   prj logvar
   set spread = exp(logvar)
   linreg(spread=spread,noprint) avgexp / resids
   # constant age ownrent income incomesq
end do iters
disp "Convergence achieved after" iters "iterations"
disp "Estimated alpha is" alpha
*
linreg(spread=spread) avgexp / resids
# constant age ownrent income incomesq
*
* This defines a function which computes the log likelihood function for an
* regression with the SPREAD option generated for a given value of "a".
*
function HLikely a
linreg(spread=income**a,noprint) avgexp
# constant age ownrent income incomesq
compute Hlikely=%logl
end
*
* First, we're going to use FIND. The one parameter over which we're searching is
* alpha, so that's the only thing which goes on the NONLIN instruction.
*
compute alpha=1.0
nonlin alpha
find max HLikely(alpha)
end find
*
* Graphical examination. We evaluate the function on a grid of points. How fine a
* grid is needed will depend upon the smoothness of the function. Note, however,
* that there is no theoretical limit on how large alpha can be, and it could also
* be negative (unlikely in this case, given the variables, but residual variances
* can often be inversely related to a regressor). We'll start with a fairly
* coarse grid running from -5 to 5.
*
* Note that the actual function values in the graph in the next differ from these
* by a constant. There are terms in the log likelihood which don't depend upon
* alpha which can be omitted or included without changing the results.
*
set agrid 1 100 = -5+.1*(t-1)
set mgrid 1 100 = HLikely(agrid)
scatter(style=lines)
# agrid mgrid
*
* Since the negatives really do seem to be unrealistic and the peak looks like
* it's short of 5.0, we'll make a finer grid from 0 to 5
*
set agrid 1 100 = .05*(t-1)
set mgrid 1 100 = HLikely(agrid)
scatter(style=lines,footer="Figure 11.2 Plot of Log-Likelihood Function",hlabel="Alpha")
# agrid mgrid
*
* Using NLLS. We need a FRML which computes the regression function, and another
* which computes the variance. The first three lines prepare the regression part.
* This defines a formula based upon the regression, naming the free parameters
* b(1),...,b(5). alpha and b are the parameters to be estimated. (The sigma
* squared is concentrated out).
*
linreg(spread=spread) avgexp
# constant age ownrent income incomesq
frml varfunc = income**alpha
frml(lastreg,vector=b) regfrml
nonlin alpha b
nlls(sv=varfunc,frml=regfrml,trace) avgexp