? ----- PROGRAM c:\racd\racd3p1.lim  July 1998 ------
 
? This LIMDEP program was \booklim\ch3r4.lim
 
? It does all the analysis in Chapter 3 (except bootstrap in Table 3.3) of
? A.C. Cameron and Pravin K. Trivedi (1998)
? REGRESSION ANALYSIS OF COUNT DATA,
? Econometric Society Monograph No.30, Cambridge University Press.
 
? A title such as Table 3.3, col.3, p.69 means Table 3.2 column 3 page 69
 
? The program runs on DOS version of Limdep 7.
? It may not be completely backward compatible with Limdep 6,
? but most of the program was originally written for Limdep 6.
? The program is not written as efficiently written as possible,
? partly because it began with Limdep 6.
? It is written in, hopefully, a transparent manner.
? Also in places intermediate output is given to enable checks.
 
? The program takes a long time to run.
? This is mainly due to use of MINIMIZE (with numerical derivatives)
? to estimate models, and use of NLLSQ and ORDERED PROBIT.
? To speed up the program the Table 3.4 estimators NB1 MLE and NB2QGPMLE
? and the Table 3.6 binary Poisson estimator have been commented out.
? To run these eliminate ? at start of line(s).

? The data set racd3.asc is the same data as originally used in
? (1) A.C. Cameron and P.K. Trivedi (1986), "Econometric Models Based on
? Count Data: Comparisons and Applications of  Some Estimators and Tests",
? Journal of Applied Econometrics, Vol. 1, pp. 29-54.
? and in Table 4
? (2) A.C. Cameron, P.K. Trivedi, F. Milne and J. Piggott (1988),
? "A Microeconometric Model of the Demand for Health Care and Health
? Insurance in Australia", Review of Economic Studies, Vol.55, pp. 85-106.
? though in places this paper used more variables than those here
? and in
? (3) A.C. Cameron, P.K. Trivedi (1993),
? "Tests of Independence in Parametric Models with Applications and
?  Illustrations", Journal of Business & Economics Statistics,
?  Vol.11, pp.29-43.
? and in a number of other subsequent papers.
 
? This data is not a representative sample of Australians as it oversamples
? young and old. In particular, use of health services may be overstated.
? This is because while the original sample of 40,650 individuals
? from the 1977-78 Australian Health Survey is representative,
? the sample used here is restricted to single people over 18 years of age.
 
? See the R.E.Stud. (1988, pp.85-106) section 3 for more detailed
? discussion of the data than that given in the RACD book.
 
? The regressors are defined in RACD on page 68
? Socioeconomic:
?  SEX       1 if female, 0 if male
?  AGE       Age in years divided by 100
?            (measured as mid-point of 10 age groups from 15-19 years to
?             65-69 with 70 or more coded treated as 72)
?  AGESQ     AGE squared
?  INCOME    Annual income in Australian dollars divided by 1000
?            (measured as mid-point of coded ranges Nil, <200, 200-1000,
?             1001-, 2001-, 3001-, 4001-, 5001-, 6001-, 7001-, 8001-10000,
?             10001-12000, 12001-14000, with 14001- treated as 15000
? Health insurance:
?  LEVYPLUS  1 if covered by private health insurance fund for private
?            patient in public hospital (with doctor of choice), 0 otherwis
?  FREEPOOR  1 if covered by government because low income, recent immigran
?            unemployed, 0 otherwise
?  FREEREPA  1 if covered free by government because of old-age or disabili
?            pension, or because invalid veteran or family of deceased
?            veteran, 0 otherwise
?  (Omitted category LEVY is 1 if covered by Medibank health insurance)
? Health status:
?  ILLNESS   Number of illnesses in past 2 weeks with 5 or more coded as 5
?  ACTDAYS   Number of days of reduced activity in past two weeks due to
?            illness or injury
?  HSCORE    General health questionnaire score using Goldberg's method.
?            High score indicates bad health.
?  CHCOND1   1 if chronic condition(s) but not limited in activity, 0 other
?  CHCOND2   1 if chronic condition(s) and limited in activity, 0 otherwise
 
? Note that the R.E.Stud. (1988) article uses different names:
?       FREEREPA was called FREEOTHER
?       HSCORE   was called GHQ
?       CHCOND1  was called LIMCHRON
?       CHCOND2  was called NONLIMCHRON
?
? The count variables are defined in R.E.Stud. (1988, pp.91)
?  DVISITS  Number of consultations with a doctor or specialist
?           in the past 2 weeks
?           (same as DOCTORCON in R.E.Stud. (1988)
?                and NOCNSLT in J.A.E.(1986))
?  NONDOCCO Number of consultations with non-doctor health professionals
?           (chemist, optician, physiotherapist, social worker,
?            district community nurse, chiropodist or chiropractor)
?            in the past 2 weeks
?  HOSPADMI Number of admissions to a hospital, psychiatric hospital,
?           nursing or convalescent home in the past 12 months
?           (up to 5 or more admissions which is coded as 5)
?  HOSPDAYS Number of nights in a hospital, etc. during most recent
?           admission: taken, where appropriate, as the mid-point of the
?           intervals 1, 2, 3, 4, 5, 6, 7, 8-14, 15-30, 31-60, 61-79
?           with 80 or more admissions coded as 80. If no admission
?           in past 12 months then equals zero.
?  MEDICINE Total number of prescribed and nonprescribed medications used
?           in past 2 days
?  PRESCRIB Total number of prescribed medications used in past 2 days
?  NONPRESC Total number of nonprescribed medications used in past 2 days
 
? +++++++++ Open output file and read in the data
OPEN; OUTPUT=c:\racd\racd3p1.out $
READ; NOBS=5190;NVAR=19;
  FILE=c:\racd\racd3.asc;
   NAMES = SEX,AGE,AGESQ,INCOME,LEVYPLUS,FREEPOOR,FREEREPA,
           ILLNESS,ACTDAYS,HSCORE,CHCOND1,CHCOND2,
           DVISITS,NONDOCCO,HOSPADMI,HOSPDAYS,
           MEDICINE,PRESCRIB,NONPRESC;
       FORMAT= (F3.0,F5.2,F7.4,F5.2,15F4.0) $
 
? +++++++++ Setting up
? Next line defines the regressors. See above for definitions.
NAMELIST; X = ONE,SEX,AGE,AGESQ,INCOME,LEVYPLUS,FREEPOOR,FREEREPA,
              ILLNESS,ACTDAYS,HSCORE,CHCOND1,CHCOND2 $
? Next line defines dependent variable - number of doctor visits
NAMELIST; Y = DVISITS $
? Next line defines dependent variable for CREATE which does
? not recognize NAMELISTs so that Y will not do
CREATE; YCREATE = DVISITS $
 
 
? +++++++++ Descriptive Statistics and frequency distribution (page 68)
 
TITLE; * RACD Table 3.1, p.68: Actual frequency distribution for DVISITS. $
HISTOGRAM; RHS=Y; INT=10 $
 
TITLE; * RACD Table 3.2, p.68: Descriptive Statistics for Y and X. $
DSTATS; RHS=Y,X $
 
? +++++++++ Poisson PMLE (page 69)
 
? --------- Poisson PMLE and possibly incorrect ML standard erros
 
TITLE; ** RACD Table 3.3, cols.1-2, p.69 : Poisson PMLE with ML st.errors $
? Poisson gives both OLS and Poisson estimates.
? The estimator maximizes (3.3) with standard errors calculated by (3.6).
POISSON; LHS=Y; RHS=X $
? Save key Poisson output for diagnostics later on
? PPRED = Prediction = EXP(X'B)
? PRES  = Residual   = Y - EXP(X'B)
? PEARSRES = Pearson residual = (Y - EXP(X'B))/SQR(EXP(X'B))
? PB = Poisson parameter estimates
CREATE; PPRED = EXP(DOT(X,B));
        PPREDSQ = PPRED*PPRED;
        PRES = YCREATE - PPRED;
        PEARSRES = (YCREATE-PPRED)/SQR(PPRED) $
MATRIX; PB = B $
CALC; LIST; PLOGL = LOGL; PKREG = KREG; PN = N $
 
? +++++++++ Various estimates of overdispersion parameter (page 69) follow
 
? --------- McCullagh-Nelder and GMT estimates of alpha
? Standard estimate for NB1 is the GLM estimate eq. (3.17)
? This uses earlier output from Poisson regression
TITLE; * RACD CH.3, p.69: GLM Estimate of NB1 Overdisp parameter. $
CALC; LIST; DEGOFFR = PN - PKREG $
MATRIX; LIST; PEARSON=XDOT(PEARSRES) $
MATRIX; LIST; ANB1MN = {1/DEGOFFR}*PEARSON $
 
? And a similar estimator for NB2 is the estimator eq. (3.19)
TITLE; * RACD CH.3, p.69: GMT Estimate of NB2 Overdisp parameter. $
CREATE; YGMT = (PRES*PRES-PPRED)/(PPRED*PPRED) $
CALC; LIST; ANB2GMT = (PN/DEGOFFR)*XBR(YGMT)$
? Check
DSTAT; RHS = PRES,PPRED,Y,YGMT $
 
 
? -------- Cameron-Trivedi estimates of Overdispersion Parameter
TITLE; * RACD CH.3, p.69: CT 1986 p46 Estimate of NB2 Overdisp param $
CREATE; YSTARCT = PRES*PRES-PPRED $
REGRESS; LHS = YSTARCT; RHS = PPREDSQ $
CALC; LIST; ANB2CT = B(1) $
 
TITLE; * RACD CH.3, p.69: CT 1986 p46 Estimate of NB1 Overdisp param $
? Same dependent as NB2 but regress on PPRED not PPREDSQ
REGRESS; LHS = YSTARCT; RHS = PPRED $
CALC; LIST; ANB1CT = B(1) $
 
 
? +++++++++ Different Standard Errors for Poisson PML Table 3.3 (page 69).
 
? (1) Usual ML (2) MLOP (3) NB1 (4) NB2 (5) Robust Sandwich (6) Bootstrap
 
? (1) ML Standard Errors (as a check, should be same as given by Poisson)
? This is eq. (3.6)
TITLE; * RACD Table 3.3, col.2 p.69: Poisson PMLE ML standard errors. $
MATRIX; UNITY = INIT(PKREG,1,1) $
MATRIX; NOLIST; XWX = XDOT(X,PPRED); VML = SINV(XWX) $
MATRIX; NOLIST; DIAGV = DIAG(VML);
        LIST; PSEML = SQRT(DIAGV) | UNITY $
MATRIX; LIST; PTML = ISQR(DIAGV) | PB $
? Isqr above is inverse square root
 
? (2) ML Outer product standard errors
CREATE; PRESSQ = PRES*PRES $
MATRIX; NOLIST; XSHATX = XDOT(X,PRESSQ);
        VMLOP = SINV(XSHATX) $
MATRIX; NOLIST; DIAGV = DIAG(VMLOP); LIST; PSEMLOP = SQRT(DIAGV) | UNITY $
MATRIX; LIST; PTMLOP = ISQR(DIAGV) | PB $
 
? (3) NB1 standard errors (method of McCullagh-Nelder)
? This is eq. (3.16)
? Tau is usual GLM estimate of alpha in (3.17) obtained earlier.
TITLE; * RACD Table 3.3, cols.4&8, p.69: Poisson PMLE NB1 standard errors. $
CALC; TAU = ANB1MN $
MATRIX; LIST; TAUSQRT = SQRT(TAU) $
MATRIX; LIST; PSENB1 = TAUSQRT*PSEML $
CALC; LIST; INVTAUSQ = 1/TAUSQRT $
MATRIX; LIST; PTNB1 = INVTAUSQ*PTML $
 
? (4) NB2 standard errors (method of GMT and CT)
? This is eq. (3.18)
? Ppred and alpha estimate using (3.19) are obtained earlier.
TITLE; * RACD Table 3.3, col.5, p.69: Poisson PMLE NB2 standard errors. $
CALC; LIST; ALPHANB2=ANB2GMT $
CREATE; NB2VAR = PPRED+ALPHANB2*PPRED*PPRED $
MATRIX; NOLIST; XWX = XDOT(X,PPRED) $
MATRIX; NOLIST; XSNB2X = XDOT(X,NB2VAR);
        VNB2 = SINV(XWX) | XSNB2X | SINV(XWX) $
MATRIX; NOLIST; DIAGV = DIAG(VNB2);
        LIST; PSENB2 = SQRT(DIAGV) | UNITY $
MATRIX; LIST; PTNB2= ISQR(DIAGV) | PB  $
? Check that Negbin 2 errors are always positive, i.e. minimum NB2VAR > 0
? (In theory with underdispersion mean + alpha*mean^2 may be < 0)
DSTAT; RHS = NB2VAR $
 
? (5) Robust Sandwich (or Eicker-White or Huber) standard errors
? This is eq. (3.20)
? Ppred and Pres are obtained earlier.
TITLE; * RACD Table 3.3, col.6, p.69: Poisson PMLE RS standard errors. $
MATRIX; NOLIST; XWX = XDOT(X,PPRED) $
CREATE; PRESSQ = PRES*PRES $
MATRIX; NOLIST; XSHATX = XDOT(X,PRESSQ);
        VRS = SINV(XWX) | XSHATX | SINV(XWX) $
MATRIX; NOLIST; DIAGV = DIAG(VRS); LIST; PSERS = SQRT(DIAGV) | UNITY $
MATRIX; LIST; PTRS = ISQR(DIAGV) | PB $
 
? (6) Bootstrap errors in column 7 of Table 3.3 on page 69
? Use separate program racd3p2.lim
 
? +++++++++ Negative binomial estimators: NB1 MLE and NB2 QGPMLE (page 76)
 
? These take quite some time as some use MINIMIZE command which is slower.
 
? -------- Negative Binomial: NB2 MLE.
? This maximizes (3.28)
TITLE; * RACD Table 3.4, cols.1,5, p.76: NB2 MLE. $
POISSON;LHS=Y;RHS=X;MODEL=N $
 
CREATE; LOGYFACT = LGM(YCREATE+1) $
 
? -------- Following checks MINIMIZE against POISSON command results
 
? Poisson using MINIMIZE as check
? This maximizes (3.3)
? TITLE; * RACD Table 3.4: Poisson MLE using MINIMIZE. $
? This gives BHHH standard errors as uses DFP algorithm
? To get Newton standard errors add  ; ALG=NEWTON  which takes longer
? MINIMIZE; LABELS = B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,B13;
?           START = -2.2,.2,1.0,-.8,-.2,.1,-.4,.1,.2,.1,.03,.1,.1;
?           FCN = EXP(X'B1) - YCREATE*X'B1 + LOGYFACT $
 
? NB2 using MINIMIZE as check
? This maximizes (3.28)
? TITLE; * RACD Table 3.4: NB2 using MINIMIZE. $
? This gives BHHH standard errors as uses DFP algorithm
? To get Newton standard errors add  ; ALG=NEWTON  which takes longer
? MINIMIZE; LABELS = B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,B13,A;
?           START = -2.2,.2,1.0,-.8,-.2,.1,-.4,.1,.2,.1,.03,.1,.1,1.0;
?           FCN = -LGM(YCREATE+(1/A)) + LOGYFACT + LGM((1/A))
?                  +(YCREATE+(1/A))*LOG(1+A*EXP(X'B1))
?                  -YCREATE*LOG(A)-YCREATE*X'B1 $
 
? -------- Negative Binomial: NB1 MLE
? This maximizes unnumbered equation at top of page 74
? Note that need to use MINIMIZE as not part of Poisson
? This gives BHHH standard errors as uses DFP algorithm
? To get Newton standard errors add  ; ALG=NEWTON  which takes longer
? To run this eliminate ? at start of next six lines
? TITLE; * RACD Table 3.4, col.3,7, p.79: NB1 MLE using MINIMIZE. $
? MINIMIZE; LABELS = B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,B13,A;
?           START = -2.2,.2,1.0,-.8,-.2,.1,-.4,.1,.2,.1,.03,.1,.1,1.0;
?           FCN = -LGM(YCREATE+(1/A)*EXP(X'B1)) + LOGYFACT
?                 + LGM((1/A)*EXP(X'B1))
?               +(YCREATE+(1/A)*EXP(X'B1))*LOG(1+A) - YCREATE*LOG(A)$
 
? -------- Negative Binomial: NB2 QGPMLE.
? This maximizes (3.36)
? This gives BHHH standard errors as uses DFP algorithm
? To get Newton standard errors add  ; ALG=NEWTON  which takes longer
? To run this eliminate ? at start of next six lines
? TITLE; * RACD Table 3.4: NB2 QGPMLE using MINIMIZE. $
? MINIMIZE; LABELS = B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,B13;
?           START = -2.2,.2,1.0,-.8,-.2,.1,-.4,.1,.2,.1,.03,.1,.1;
?           FCN = (1/.286)*LOG(1+.286*EXP(X'B1))
?                - YCREATE*LOG(.286*EXP(X'B1))
?                + YCREATE*LOG(1+.286*EXP(X'B1)) $
 
? ---------- Negative Binomial NB1 QGPMLE.
 
? This is eq. (3.3) minimized with standard errors (3.16)
? Same as Poisson PMLE in Table 3.3 with corrected NB1 standard errors $
 
 
? +++++++++ LM Overdispersion Tests (page 79)
 
? Regression based tests for overdispersion
? Cameron and Trivedi J. Econometrics 1990
? For overdispersion do "t-test" of H0: alpha=0 vs. H1: alpha > 0
? For underdispersion do "t-test" of H0: alpha=0 vs. H1: alpha < 0
 
TITLE; * RACD CH.3, p.79: LM Overdispersion Test of NB2 form. $
? Negbin 2 form is var = mean + alpha*mean*mean
? Equation estimated is (3.39). LM test is t-statistic for the constant.
CREATE; YSTARLM = (PRES*PRES-YCREATE)/PPRED $
REGRESS; LHS = YSTARLM; RHS = PPRED $
? Following is saved as an alternative estimate of alpha
CALC; LIST; ANB2LM = B(1) $
 
TITLE; * RACD CH.3, p.79: LM Overdispersion Test of NB1 form. $
? Negbin 1 form is var = mean + alpha*mean
? Equation estimated is (3.40). LM test is t-statistic for the constant.
? Same dependent as NB2 but regress on ONE not PPRED
REGRESS; LHS = YSTARLM; RHS = ONE $
? Following is saved as an alternative estimate of alpha
CALC; LIST; ANB1LM = B(1) $
 
 
 
? +++++++++ Effect of regressors (page 83).
 
? Table 3.5, col.1, p.83 is just the Poisson PMLE estimates
 
TITLE; * RACD Table 3.5, col.2, p.83: Average of individual effects $
? This is (3.44) used to calculate (3.43)
? Since constant included, can do the shortcut method equation 3.44
? and use the earlier Poisson estimator scaled by Ybar.
CALC; LIST; YMEAN=XBR(YCREATE) $
MATRIX; LIST; PMEANEFF = YMEAN*PB $
 
TITLE; * RACD Table 3.5, col.3, p.83: Individual effect at Average value $
? This is (3.45)
MATRIX; NOLIST; XMEAN=MEAN(X) $
CALC; LIST; PRDXMEAN=EXP(XMEAN'PB) $
MATRIX; LIST; PEFFMEAN = PRDXMEAN*PB $
 
? RACD Table 3.5, col.4, p.83 is just the OLS estimates
 
TITLE; * RACD Table 3.5, col.5, p.83: Elasticity  $
? Semi-standardised coefficient. Scale beta_hat by mean of x.
MATRIX; LIST; PELAST = DIAG(XMEAN) | PB $
 
TITLE; * RACD Table 3.5, col.6, p.83: Semi-standardised coeff.  $
? Semi-standardised coefficient. Scale beta_hat by standard deviation of x.
MATRIX; NOLIST; XSTDEV = SDEV(X) $
MATRIX; LIST; PSEMIST = DIAG(XSTDEV) | PB $
 
TITLE; * RACD Table 3.5, col.7&8, p.83: Mean and variance of X $
MATRIX; NOLIST; XMEAN=MEAN(X) $
MATRIX; NOLIST; XSTDEV = SDEV(X) $
 
 
 
? +++++++++ Other Models for counts (page 92)
 
? -------- Binary probit and logit given to compare with binary Poisson
 
TITLE; ***** RACD CH.3: Binary Probit and logit to compare with bin Poiss. $
CREATE; YONE = 0;
        IF (YCREATE > 0) YONE = 1 $
DSTAT; RHS=Y,YONE $
PROBIT; LHS=YONE; RHS=X $
LOGIT; LHS=YONE; RHS=X $
 
? -------- Binary Poisson
 
TITLE; * RACD Table 3.6, col.1 (in theory), p.92: Binary Poisson. $
? This gave results different from Table 3.6 in the book!
? Most likely the book version is incorrect!
? Maximizes (3.54) with F( ) defined in unnumbered equation after (3.54)
? Pr(0) = exp(-lamda) = exp(-exp(x'b))
? So Pr(1) = 1-exp(-lamda) = 1-exp(-exp(x'b))
? This gives BHHH standard errors as uses DFP algorithm
? To get Newton standard errors add  ; ALG=NEWTON  which takes longer
MINIMIZE; LABELS = B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,B13;
          START = -.9,.2,-1.3,1.8,.0,.1,-.3,.2,.2,.1,.03,.1,.1;
          FCN = -YONE*LOG(1-EXP(-EXP(X'B1)))
                  - (1-YONE)*(-EXP(X'B1)) $
 
? -------- Ordered Probit and rescaling of ordered probit
 
? Next get estimate of OLS standard deviation to rescale ordered probit
TITLE; * RACD CH.3.6, col.3, p.92: OLS no transformation. $
REGRESS; LHS=Y; RHS=X; RES=RY; HETERO $
CALC; LIST; OLSSSQRD = SSQRD $
 
TITLE; * RACD Table 3.6: Ordered Probit before rescaling. $
? Since there are no 8's combine 8 and 9 into one category
CREATE; YOP = YCREATE;
        IF (YCREATE=8 | YCREATE=9) YOP = 8 $
DSTAT; RHS=Y,YOP $
ORDERED PROBIT; LHS=YOP; RHS=X $
 
TITLE; * RACD Table 3.6, col.2, p.92: Ordered Probit rescaled by OLS stdev.$
MATRIX; OPB = B $
CALC; LIST; OLSS = SQR(OLSSSQRD) $
MATRIX; LIST; RESCOPB = OLSS*OPB  $
 
? -------- OLS on Transformations
TITLE; * RACD Table 3.6, col.4, p.92: OLS of log(y+0.1) transformation. $
CREATE; LOGYPT1 = LOG(YCREATE + 0.1) $
REGRESS; LHS=LOGYPT1; RHS=X; RES=RLOGYPT1; KEEP=PLOGYPT1; HETERO $
 
TITLE; * RACD Table 3.6: OLS of log(y+0.2) transformation. $
CREATE; LOGYPT2 = LOG(YCREATE + 0.2) $
REGRESS; LHS=LOGYPT2; RHS=X; RES=RLOGYPT2; KEEP=PLOGYPT2; HETERO $
 
TITLE; * RACD Table 3.6: OLS of log(y+0.4) transformation. $
CREATE; LOGYPT4 = LOG(YCREATE + 0.4) $
REGRESS; LHS=LOGYPT4; RHS=X; RES=RLOGYPT4; KEEP=PLOGYPT4; HETERO $
 
TITLE; * RACD Table 3.6, col.5, p.92: OLS of Square root(y) transformation $
CREATE; SQRTY = SQR(YCREATE) $
REGRESS; LHS=SQRTY; RHS=X; RES=RSQRTY; KEEP=PSQRTY; HETERO $
 
? Table 3.6, col.6, p92 is just the Poisson estimates
 
TITLE; * RACD Table 3.6, cols.3-6: Skewness and Kurtosis of Residuals
DSTAT; RHS=Y,LOGYPT1,LOGYPT2,LOGYPT4,SQRTY,
           RY,RLOGYPT1,RLOGYPT2,RLOGYPT4,RSQRTY,PRES,NLSRES $
 
TITLE; * RACD Table 3.6: Properties of Predictions $
CREATE; TLOGYPT1 = EXP(PLOGYPT1) - 0.1;
        TLOGYPT2 = EXP(PLOGYPT2) - 0.2;
        TLOGYPT4 = EXP(PLOGYPT4) - 0.4;
        TSQRTY = PSQRTY*PSQRTY $
DSTAT; RHS=YCREATE,PPRED,TLOGYPT1,TLOGYPT2,TLOGYPT4,TSQRTY$
 
? -------- NLS with exponential mean
 
TITLE; * RACD Table 3.6, col.7, p.92: NLS with mean exp(xb). $
NLSQ; LHS=Y; START=PB;
      LABELS=BONE,BSEX,BAGE,BAGESQ,BINCOME,BLEVYPLU,BFREEPOO,
             BFREEREP,BILLNESS,BACTDAYS,BHSCORE,BCHCOND1,BCHCOND2;
      FCN=EXP(DOT[X]) $
? Get Robust Sandwich standard errors
? This is not currently working
MATRIX; NLSB = B $
CREATE; NLSPRED = EXP(DOT(X,NLSB));
        NLSRES = YCREATE - NLSPRED;
        NLSRESSQ = NLSRES*NLSRES;
        NLSPRESQ = NLSPRED*NLSPRED;
        NLSNEWSQ = NLSRESSQ*NLSPRESQ  $
MATRIX; NOLIST; XWX = XDOT(X,NLSPRESQ) $
MATRIX; NOLIST; XSHATX = XDOT(X,NLSNEWSQ);
        VEW = SINV(XWX) | XSHATX | SINV(XWX) $
MATRIX; NOLIST; DIAGV = DIAG(VEW);
        LIST; NLSSEEW = SQRT(DIAGV) | UNITY $
MATRIX; lIST; NLSTEW = ISQR(DIAGV) | NLSB $
 
STOP $