** Chapter 23.4 ** for W.H. Greene, Econometric Analysis 6th ed. *****************
* (c) Noel Roy 2003, 2008
*
*                         MODELS FOR DISCRETE CHOICE
*
*===============================================================================
* 23.4 ESTIMATION AND INFERENCE IN BINARY CHOICE MODELS
*
* Example 23.3 (p. 781) Probability Models
*
* Read Spector-Mazzeo Program Effectiveness Data
READ (TableF16-1.prn) /NAMES
*
* Estimate the probability models in Table 23.1.
* Linear Model:
OLS GRADE GPA TUCE PSI /COEF=b
* Logistic Model:
LOGIT GRADE GPA TUCE PSI /COEF=bL COV=VL
GEN1 K=$K
* Probit Model:
PROBIT GRADE GPA TUCE PSI /COEF=BP COV=VP
GEN1 LLFPROB=$LLF
* Marginal effects and scale factors are reported by both commands. See Chapter
* 25 of the SHAZAM Manual for more information on the LOGIT and PROBIT commands.
*
* The Gumbel model is not specifically supported but can be estimated using the
* NL command with the LOGDEN option. The OLS estimates b are used as starting values.
NL 1 /NCOEF=4 LOGDEN START=b COEF=BETA
* The likelihood function exp(-exp(-x'Beta)) for the Gumbel model is given on p. 774. 
EQ (GRADE-1)*exp(b1*GPA+b2*TUCE+b3*PSI+b4) + GRADE*log(1-exp(-exp(b1*GPA+b2*TUCE+b3*PSI+b4)))
END
* Use (23-9) to calculate the marginal effects (labelled "Slope" in Table 23.1)
* of the Gumbel model. 
GENR ONE=1
* Calculate value of x at the sample means.
?STAT GPA TUCE PSI ONE  /MEAN=XBAR
* Cross multiply vector of means by Beta.
MATRIX BXBAR=BETA'XBAR
* Calculate the Gumbel density function f(xbar'Beta) (obtained by differentiating Prob(Y=1|x) )
MATRIX SCALE=exp(-exp(bxbar))*exp(bxbar)
PRINT SCALE
* Drop the constant term from Beta, since it has no significance as a "marginal" effect.
COPY BETA BETA13 /FROW=1;3 /TROW=1;3
MATRIX SLOPE=SCALE*BETA13
PRINT SLOPE
DELETE ONE
*
* Replicate Figure 23.2 using the GRAPH command. First generate the probability
* functions. The NCDF(x) function with the GENR command returns the cumulative
* normal probability of a standard normal variable. The probability will be
* in the zero-one range. Then sort the data in ascending order using the SORT
* command with GPA as the sort variable. Include the OBS variable (i.e. a time
* trend variable) to allow for unsorting the data back to its original order.
*
?STAT TUCE / MEAN=MEANTUCE
GENR WITHPSI=NCDF(BP:4+BP:1*GPA+BP:2*MEANTUCE+BP:3*1)
GENR WOPSI=NCDF(BP:4+BP:1*GPA+BP:2*MEANTUCE)
SORT GPA WITHPSI WOPSI OBS
GRAPH WITHPSI WOPSI GPA /LINEONLY
*
* Now unsort the data.
*
SORT OBS GPA
*
* The coefficients, t ratios, and slopes (Marginal effects) are automatically
* generated by the PROBIT and LOGIT commands, but not the standard errors and t ratios
* of the marginal effects, as in Table 23.2.
* Use the formulas in Section 23.4.2 to calculate these.
* Compute marginal effects at the means using the formal results.
* For the moment, ignore the fact that the 4th variable in X is a
* dummy variable.
*
* Logit Model
* Calculate Beta'x.
MATRIX BLX=bL'XBAR
* Calculate the Logstic Probabilities Prob(GRADE=1|x)
MATRIX PL=1/(1+EXP(-BLX))
* Logistic density function=P(Beta x) (1-P(Beta x))
MATRIX FL=PL*(1-PL)
* Asymptotic covariance matrix of the marginal effects Gamma_hat.
MATRIX GL=FL(IDEN(K)+(1-2*PL)*BL*XBAR')
MATRIX VGAMMAL=GL*VL*GL'
* Cacculate the standard errors from the diagonal elements.
MATRIX SLOPESEL=SQRT(DIAG(VGAMMAL))
MATRIX SLOPETL=FL*BL/SLOPESEL
PRINT SLOPESEL SLOPETL
*
* Now do the same for the Probit Model.
MATRIX BPX=BP'XBAR
MATRIX PP=NCDF(BPX)
?DISTRIB BPX /PDF=FP BEG=1 END=1
MATRIX GP=FP*(IDEN(K)-BPX*BP*XBAR')
MATRIX VGAMMAP=GP*VP*GP'
MATRIX SLOPESEP=SQRT(DIAG(VGAMMAP))
MATRIX SLOPETP=FP*BP/SLOPESEP
PRINT SLOPESEP SLOPETP

* Allow for fact PSI is binary, using equations (23-25) and (23-26).
* Xbar when PSI=0
MATRIX  XBAR0 = XBAR 
MATRIX XBAR0(3)=0 
* Xbar when PSI=1
MATRIX  XBAR1 = XBAR 
MATRIX XBAR1(3)=1 
* Logit model
* Redo calculations with PSI=0
MATRIX BLX0=bL'XBAR0
MATRIX PL0=1/(1+EXP(-BLX0))
MATRIX FL0=PL0*(1-PL0)
* Ditto for PSI=1.
MATRIX BLX1=bL'XBAR1
MATRIX PL1=1/(1+EXP(-BLX1))
MATRIX FL1=PL1*(1-PL1)
* Calculate the "marginal effict" for a change in PSI from 0 to 1.
MATRIX DELFL = PL1-PL0
* Calculate asy var deltaF as per (23-26).
MATRIX SEDELFL=SQRT(FL1**2*XBAR1'VL*XBAR1+FL0**2*XBAR0'VL*XBAR0)
MATRIX TRDELFL=DELFL/SEDELFL
PRINT DELFL SEDELFL TRDELFL
* Do same with Probit Model.
MATRIX BPX0=bP'XBAR0
MATRIX PP0=NCDF(BPX0)
?DISTRIB BPX0 /PDF=FP0 BEG=1 END=1
MATRIX BPX1=bP'XBAR1
MATRIX PP1=NCDF(BPX1)
?DISTRIB BPX1 /PDF=FP1 BEG=1 END=1
MATRIX DELFP = PP1-PP0
MATRIX SEDELFP=SQRT(FP1**2*XBAR1'VP*XBAR1+FP0**2*XBAR0'VP*XBAR0)
MATRIX TRDELFP=DELFP/SEDELFP
PRINT DELFP SEDELFP TRDELFP
*
STOP
*
*===============================================================================
Updated November 26, 2008