. di "stmle.do by Colin Cameron: Basic Stata MLE example using Poisson" stmle.do by Colin Cameron: Basic Stata MLE example using Poisson . . . ********** MEMORY MANAGEMENT . set maxvar 100 width 1000 (maxvar and maxobs no longer need be set with this version of Stata) . . . ********** READ DATA . * You need file jaggia.asc in your directory . * Infile: FREE FORMAT WITHOUT DICTIONARY . * As there is space between each observation data is also space-delimited . * free format and then there is no need for a dictionary file . * The following command spans more that one line so use /* and */ . . infile docno weeks numbids takeover bidprem insthold size leglrest /* > */ realrest finrest regulatn whtknght using jaggia.asc (126 observations read) . . . ********** DATA TRANSFORMATIONS . gen sizesq = size*size . label variable sizesq "size squared" . . . ******** CHECK DATA: DESCRIPTIVE STATISTICS . describe Contains data obs: 126 vars: 13 size: 7,056 (99.2% of memory free) ------------------------------------------------------------------------------- 1. docno float %9.0g 2. weeks float %9.0g 3. numbids float %9.0g 4. takeover float %9.0g 5. bidprem float %9.0g 6. insthold float %9.0g 7. size float %9.0g 8. leglrest float %9.0g 9. realrest float %9.0g 10. finrest float %9.0g 11. regulatn float %9.0g 12. whtknght float %9.0g 13. sizesq float %9.0g size squared ------------------------------------------------------------------------------- Sorted by: Note: dataset has changed since last saved . summarize Variable | Obs Mean Std. Dev. Min Max ---------+----------------------------------------------------- docno | 126 82174.41 2251.783 78001 85059 weeks | 126 11.44898 7.711424 2.857 41.429 numbids | 126 1.738095 1.432081 0 10 takeover | 126 1 0 1 1 bidprem | 126 1.346806 .189325 .9426754 2.066366 insthold | 126 .2518175 .1856136 0 .904 size | 126 1.219031 3.096624 .017722 22.169 leglrest | 126 .4285714 .4968472 0 1 realrest | 126 .1825397 .3878308 0 1 finrest | 126 .1031746 .3054011 0 1 regulatn | 126 .2698413 .4456492 0 1 whtknght | 126 .5952381 .4928054 0 1 sizesq | 126 10.99902 59.91479 .0003141 491.4646 . . . ********** POISSON REGRESSION . * . * Use /* and */ as command spans two lines . . poisson numbids leglrest realrest finrest whtknght /* > */ bidprem insthold size sizesq regulatn Iteration 0: log likelihood = -184.9518 Iteration 1: log likelihood = -184.94833 Iteration 2: log likelihood = -184.94833 Poisson regression Number of obs = 126 LR chi2(9) = 33.25 Prob > chi2 = 0.0001 Log likelihood = -184.94833 Pseudo R2 = 0.0825 ------------------------------------------------------------------------------ numbids | Coef. Std. Err. z P>|z| [95% Conf. Interval] ---------+-------------------------------------------------------------------- leglrest | .2601464 .1509594 1.723 0.085 -.0357286 .5560213 realrest | -.1956597 .1926309 -1.016 0.310 -.5732093 .1818899 finrest | .07403 .2165219 0.342 0.732 -.3503452 .4984053 whtknght | .4813822 .1588698 3.030 0.002 .170003 .7927613 bidprem | -.6776959 .3767373 -1.799 0.072 -1.416087 .0606956 insthold | -.3619913 .4243292 -0.853 0.394 -1.193661 .4696788 size | .1785026 .0600221 2.974 0.003 .0608614 .2961438 sizesq | -.0075693 .0031217 -2.425 0.015 -.0136878 -.0014509 regulatn | -.0294392 .1605682 -0.183 0.855 -.344147 .2852686 _cons | .9860599 .5339202 1.847 0.065 -.0604044 2.032524 ------------------------------------------------------------------------------ . . . ********** MAXIMUM LIKELIHOOD FOR POISSON . . * Based on Weibull example in [R] ml - maximum likelihood . * See especially Stata 6 Manual Reference H-O p.379 and p.388-399 . * The possible methods are lf, d0, d1 and d2. . * Here I use d0, d1 and d2 which give log-density etc for each observation. . * I do not use lf which gives the log-likeihood over all observations . . * The general form for the command is . * (1) Define the program, log-density, perhaps gradient and perhaps hessian > . * program define myprog . * : : . * end . * (2) Run the program providing dependent variable(s) and regressor(s) . * ml model d*** myprog (dep1 dep2 ... = x's for eqn1) (x's for eqn2) .. > . * ml maximize . * where d*** called by d0, d1, d1debug, d2 and d2debug . * and the debug variants check anayltical against numerical derivatives . * and for robust standard errors based on A-in*B*A-inv add , robust . . . * For the Poisson log-density . * ln L = Sum_i {-mu_i + y_i*ln(mu_i) - ln(y_i!)} . * where . * mu_i = exp(x_i'b) . * . * The first-order conditions are . * Sum_i (y_i - exp(x_i'b))x_i = 0 . * . * and the second-derivative matrix is . * Sum_i -exp(x_i'b))*x_i*x_i' . . . * (1) Define the program, log-density, perhaps gradient and perhaps hessian > . * The Poisson model is simple as only one index and one dependent variable . * There is one dependent variable y which is stored in $ML_y1 . * There is one index which is x'b and for the Poisson model is mu i.e. E[y|x] > . * For log-density -mu_i + y_i*ln(mu_i) - ln(y_i!) where mu_i = exp(x_i'b) . * For gradient d/db = (y_i - mu_i)*x_i . * but Stata wands d/dtheta_i = (y_i - mu_i) using theta_i = x_i'b . * For hessian mu_i*x_i*x_i' . * but Stata wants d2/dtheta_i^2 = mu_i using theta_i = x_i'b . . * (2) Example is . * ml model d1debug mypois1 (numbids = bidprem) . * ml maximize . . * For program debuggin which produces much output . * set trace on . . * The following programs are given . * mypois0 use numerical derivatives throughout . * mypois1 use analytical first derivative and numerical second . * mypois1r same as mypois1 but form robust variance matrix . * mypois2 use analytical first and second derivative . * mypois1 same as mypois1 but form robust variance matrix . . . ******* (1) MYPOIS0 PROGRAM: POISSON WITH NUMERICAL DERIVATIVES THROUGHOUT . . program define mypois0 1. version 6.0 2. args todo b lnf /* Need to use the names todo b and lnf > todo always contains 1 and may be ignored > b is parameters and lnf is log-density */ 3. tempvar theta1 /* create as needed to calculate lf, g, ... */ 4. mleval `theta1' = `b', eq(1) /* theta1 is theta1_i = x_i'b */ 5. local y "$ML_y1" /* create to make program more readable */ 6. tempvar lnyfact mu 7. quietly gen double `lnyfact' = lnfact(`y') 8. quietly gen double `mu' = exp(`theta1') 9. mlsum `lnf' = -`mu' + `y'*ln(`mu') - `lnyfact' 10. end . . ******* (1) OUTPUT: POISSON WITH NUMERICAL DERIVATIVES THROUGHOUT . . ml model d0 mypois0 (numbids = leglrest realrest finrest whtknght /* > */ bidprem insthold size sizesq regulatn) . ml maximize initial: log likelihood = -229.63257 alternative: log likelihood = -201.87145 rescale: log likelihood = -201.87145 Iteration 0: log likelihood = -201.87145 Iteration 1: log likelihood = -187.1169 Iteration 2: log likelihood = -184.95692 Iteration 3: log likelihood = -184.94833 Iteration 4: log likelihood = -184.94833 Number of obs = 126 Wald chi2(9) = 33.79 Log likelihood = -184.94833 Prob > chi2 = 0.0001 ------------------------------------------------------------------------------ numbids | Coef. Std. Err. z P>|z| [95% Conf. Interval] ---------+-------------------------------------------------------------------- leglrest | .2601464 .1509594 1.723 0.085 -.0357286 .5560214 realrest | -.1956597 .192631 -1.016 0.310 -.5732096 .1818901 finrest | .07403 .216522 0.342 0.732 -.3503452 .4984052 whtknght | .4813821 .1588698 3.030 0.002 .1700029 .7927613 bidprem | -.6776959 .3767378 -1.799 0.072 -1.416088 .0606967 insthold | -.3619913 .4243294 -0.853 0.394 -1.193662 .4696792 size | .1785026 .0600224 2.974 0.003 .0608608 .2961444 sizesq | -.0075694 .0031217 -2.425 0.015 -.0136878 -.0014509 regulatn | -.0294392 .1605682 -0.183 0.855 -.3441471 .2852686 _cons | .98606 .5339209 1.847 0.065 -.0604058 2.032526 ------------------------------------------------------------------------------ . . . ******* (2) MYPOIS1 PROGRAM: POISSON WITH ANAYLTICAL FIRST DERIVATIVES . . program define mypois1 1. version 6.0 2. args todo b lnf g /* Need to use the names todo b and lnf > todo always contains 1 and may be ignored > b is parameters, lnf is log-density > g is gradient */ 3. tempvar theta1 /* create as needed to calculate lf, g, ... */ 4. mleval `theta1' = `b', eq(1) /* theta1 is theta1_i = x_i'b */ 5. local y "$ML_y1" /* create to make program more readable */ 6. tempvar lnyfact mu 7. quietly gen double `lnyfact' = lnfact(`y') 8. quietly gen double `mu' = exp(`theta1') 9. mlsum `lnf' = -`mu' + `y'*ln(`mu') - `lnyfact' 10. /* Following code is extra for analytical first derivatives */ . if `todo'==0 | `lnf'==. {exit} 11. tempname d1 12. mlvecsum `lnf' `d1' = `y' - `mu' 13. matrix `g' = `d1' 14. end . . ******* (2) MYPOIS1 OUTPUT: POISSON WITH ANALYTICAL FIRST DERIVATIVES . . ml model d1debug mypois1 (numbids = leglrest realrest finrest whtknght /* > */ bidprem insthold size sizesq regulatn) . ml maximize initial: log likelihood = -229.63257 alternative: log likelihood = -201.87145 rescale: log likelihood = -201.87145 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = .0011843 d1debug: End derivative-comparison report ----------------------------------- Iteration 0: log likelihood = -201.87145 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = 6.29e-09 d1debug: End derivative-comparison report ----------------------------------- Iteration 1: log likelihood = -187.1169 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = 6.83e-07 d1debug: End derivative-comparison report ----------------------------------- Iteration 2: log likelihood = -184.95692 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = .0003457 d1debug: End derivative-comparison report ----------------------------------- Iteration 3: log likelihood = -184.94833 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = .0004541 d1debug: End derivative-comparison report ----------------------------------- Iteration 4: log likelihood = -184.94833 Number of obs = 126 Wald chi2(9) = 33.79 Log likelihood = -184.94833 Prob > chi2 = 0.0001 ------------------------------------------------------------------------------ numbids | Coef. Std. Err. z P>|z| [95% Conf. Interval] ---------+-------------------------------------------------------------------- leglrest | .2601464 .1509594 1.723 0.085 -.0357286 .5560214 realrest | -.1956597 .192631 -1.016 0.310 -.5732096 .1818901 finrest | .07403 .216522 0.342 0.732 -.3503452 .4984052 whtknght | .4813821 .1588698 3.030 0.002 .1700029 .7927613 bidprem | -.6776959 .3767378 -1.799 0.072 -1.416088 .0606967 insthold | -.3619913 .4243294 -0.853 0.394 -1.193662 .4696792 size | .1785026 .0600224 2.974 0.003 .0608608 .2961444 sizesq | -.0075694 .0031217 -2.425 0.015 -.0136878 -.0014509 regulatn | -.0294392 .1605682 -0.183 0.855 -.3441471 .2852686 _cons | .98606 .5339209 1.847 0.065 -.0604058 2.032526 ------------------------------------------------------------------------------ . . . ******* (3) MYPOIS1r PROGRAM: ROBUST SE'S FOR POISSON ANAYLTICAL FIRST . . program define mypois1r 1. version 6.0 2. args todo b lnf g negH g1 /* Need to use the names todo b and lnf > todo always contains 1 and may be ignored > b is parameters, lnf is log-density > g is gradient negH need not be provided here > g1 is to be used for robust se's */ 3. tempvar theta1 /* create as needed to calculate lf, g, ... */ 4. mleval `theta1' = `b', eq(1) /* theta1 is theta1_i = x_i'b */ 5. local y "$ML_y1" /* create to make program more readable */ 6. tempvar lnyfact mu 7. quietly gen double `lnyfact' = lnfact(`y') 8. quietly gen double `mu' = exp(`theta1') 9. mlsum `lnf' = -`mu' + `y'*ln(`mu') - `lnyfact' 10. /* Following code is extra for analytical first derivatives > and also changed due to robust se's */ . if `todo'==0 | `lnf'==. {exit} 11. tempname d1 12. quietly replace `g1' = `y' - `mu' 13. mlvecsum `lnf' `d1' = `g1' 14. matrix `g' = `d1' 15. end . . ******* (3) MYPOIS1r OUTPUT: ROBUST SE'S FOR POISSON ANAYLTICAL FIRST . . ml model d1debug mypois1r (numbids = leglrest realrest finrest whtknght /* > */ bidprem insthold size sizesq regulatn), robust . ml maximize initial: log likelihood = -229.63257 alternative: log likelihood = -201.87145 rescale: log likelihood = -201.87145 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = .0011843 d1debug: End derivative-comparison report ----------------------------------- Iteration 0: log likelihood = -201.87145 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = 6.29e-09 d1debug: End derivative-comparison report ----------------------------------- Iteration 1: log likelihood = -187.1169 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = 6.83e-07 d1debug: End derivative-comparison report ----------------------------------- Iteration 2: log likelihood = -184.95692 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = .0003457 d1debug: End derivative-comparison report ----------------------------------- Iteration 3: log likelihood = -184.94833 d1debug: Begin derivative-comparison report --------------------------------- d1debug: mreldif(gradient vector) = .0004541 d1debug: End derivative-comparison report ----------------------------------- Iteration 4: log likelihood = -184.94833 Number of obs = 126 Wald chi2(9) = 34.98 Log likelihood = -184.94833 Prob > chi2 = 0.0001 ------------------------------------------------------------------------------ | Robust numbids | Coef. Std. Err. z P>|z| [95% Conf. Interval] ---------+-------------------------------------------------------------------- leglrest | .2601464 .1250534 2.080 0.037 .0150462 .5052465 realrest | -.1956597 .1816168 -1.077 0.281 -.5516222 .1603027 finrest | .07403 .2635711 0.281 0.779 -.4425598 .5906198 whtknght | .4813821 .1064948 4.520 0.000 .2726562 .690108 bidprem | -.6776959 .297425 -2.279 0.023 -1.260638 -.0947536 insthold | -.3619913 .3231802 -1.120 0.263 -.9954127 .2714302 size | .1785026 .062355 2.863 0.004 .056289 .3007161 sizesq | -.0075694 .0027788 -2.724 0.006 -.0130158 -.0021229 regulatn | -.0294392 .1420508 -0.207 0.836 -.3078536 .2489752 _cons | .98606 .4137396 2.383 0.017 .1751453 1.796975 ------------------------------------------------------------------------------ . . . ******* (4) MYPOIS2 PROGRAM: POISSON WITH ANAYLTICAL SECOND DERIVATIVES . . program define mypois2 1. version 6.0 2. args todo b lnf g negH /* Need to use the names todo b and lnf > todo always contains 1 and may be ignored > b is parameters, lnf is log-density > g is gradient, negH is negative hessian */ 3. tempvar theta1 /* create as needed to calculate lf, g, ... */ 4. mleval `theta1' = `b', eq(1) /* theta1 is theta1_i = x_i'b */ 5. local y "$ML_y1" /* create to make program more readable */ 6. tempvar lnyfact mu 7. quietly gen double `lnyfact' = lnfact(`y') 8. quietly gen double `mu' = exp(`theta1') 9. mlsum `lnf' = -`mu' + `y'*ln(`mu') - `lnyfact' 10. /* Following code is extra for analytical first derivatives */ . if `todo'==0 | `lnf'==. {exit} 11. tempname d1 12. mlvecsum `lnf' `d1' = `y' - `mu' 13. matrix `g' = `d1' 14. /* Following code is extra for analytical second derivatives */ . if `todo'==0 | `lnf'==. {exit} 15. tempname d11 16. mlmatsum `lnf' `d11' = `mu' 17. matrix `negH' = `d11' 18. end . . ******* (4) MYPOIS2 OUTPUT: POISSON WITH ANALYTICAL SECOND DERIVATIVES . . ml model d2debug mypois2 (numbids = leglrest realrest finrest whtknght /* > */ bidprem insthold size sizesq regulatn) . ml maximize initial: log likelihood = -229.63257 alternative: log likelihood = -201.87145 rescale: log likelihood = -201.87145 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = .0011843 d2debug: mreldif(negative Hessian) = .0006026 d2debug: End derivative-comparison report ----------------------------------- Iteration 0: log likelihood = -201.87145 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = 6.29e-09 d2debug: mreldif(negative Hessian) = .0419201 d2debug: End derivative-comparison report ----------------------------------- Iteration 1: log likelihood = -187.1169 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = 6.83e-07 d2debug: mreldif(negative Hessian) = 3.48e-06 d2debug: End derivative-comparison report ----------------------------------- Iteration 2: log likelihood = -184.95692 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = .0003457 d2debug: mreldif(negative Hessian) = 6.50e-07 d2debug: End derivative-comparison report ----------------------------------- Iteration 3: log likelihood = -184.94833 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = .0004541 d2debug: mreldif(negative Hessian) = 6.54e-07 d2debug: End derivative-comparison report ----------------------------------- Iteration 4: log likelihood = -184.94833 Number of obs = 126 Wald chi2(9) = 33.79 Log likelihood = -184.94833 Prob > chi2 = 0.0001 ------------------------------------------------------------------------------ numbids | Coef. Std. Err. z P>|z| [95% Conf. Interval] ---------+-------------------------------------------------------------------- leglrest | .2601464 .1509594 1.723 0.085 -.0357286 .5560214 realrest | -.1956597 .192631 -1.016 0.310 -.5732096 .1818901 finrest | .07403 .216522 0.342 0.732 -.3503452 .4984052 whtknght | .4813821 .1588698 3.030 0.002 .1700029 .7927613 bidprem | -.6776959 .3767378 -1.799 0.072 -1.416088 .0606967 insthold | -.3619913 .4243294 -0.853 0.394 -1.193662 .4696792 size | .1785026 .0600224 2.974 0.003 .0608608 .2961444 sizesq | -.0075694 .0031217 -2.425 0.015 -.0136878 -.0014509 regulatn | -.0294392 .1605682 -0.183 0.855 -.3441471 .2852686 _cons | .98606 .5339209 1.847 0.065 -.0604058 2.032526 ------------------------------------------------------------------------------ . . . ******* (5) MYPOIS2r PROGRAM: ROBUST SE'S FOR POISSON ANAYLTICAL SECOND . . program define mypois2r 1. version 6.0 2. args todo b lnf g negH g1 /* Need to use the names todo b and lnf > todo always contains 1 and may be ignored > b is parameters, lnf is log-density > g is gradient negH need not be provided here > g1 is to be used for robust se's */ 3. tempvar theta1 /* create as needed to calculate lf, g, ... */ 4. mleval `theta1' = `b', eq(1) /* theta1 is theta1_i = x_i'b */ 5. local y "$ML_y1" /* create to make program more readable */ 6. tempvar lnyfact mu 7. quietly gen double `lnyfact' = lnfact(`y') 8. quietly gen double `mu' = exp(`theta1') 9. mlsum `lnf' = -`mu' + `y'*ln(`mu') - `lnyfact' 10. /* Following code is extra for analytical first derivatives > and also changed due to robust se's */ . if `todo'==0 | `lnf'==. {exit} 11. tempname d1 12. quietly replace `g1' = `y' - `mu' 13. mlvecsum `lnf' `d1' = `g1' 14. matrix `g' = `d1' 15. /* Following code is extra for analytical second derivatives */ . if `todo'==0 | `lnf'==. {exit} 16. tempname d11 17. mlmatsum `lnf' `d11' = `mu' 18. matrix `negH' = `d11' 19. end . . ******* (3) MYPOIS2r OUTPUT: ROBUST SE'S FOR POISSON ANAYLTICAL SECOND . . ml model d2debug mypois2r (numbids = leglrest realrest finrest whtknght /* > */ bidprem insthold size sizesq regulatn), robust . ml maximize initial: log likelihood = -229.63257 alternative: log likelihood = -201.87145 rescale: log likelihood = -201.87145 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = .0011843 d2debug: mreldif(negative Hessian) = .0006026 d2debug: End derivative-comparison report ----------------------------------- Iteration 0: log likelihood = -201.87145 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = 6.29e-09 d2debug: mreldif(negative Hessian) = .0419201 d2debug: End derivative-comparison report ----------------------------------- Iteration 1: log likelihood = -187.1169 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = 6.83e-07 d2debug: mreldif(negative Hessian) = 3.48e-06 d2debug: End derivative-comparison report ----------------------------------- Iteration 2: log likelihood = -184.95692 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = .0003457 d2debug: mreldif(negative Hessian) = 6.50e-07 d2debug: End derivative-comparison report ----------------------------------- Iteration 3: log likelihood = -184.94833 d2debug: Begin derivative-comparison report --------------------------------- d2debug: mreldif(gradient vector) = .0004541 d2debug: mreldif(negative Hessian) = 6.54e-07 d2debug: End derivative-comparison report ----------------------------------- Iteration 4: log likelihood = -184.94833 Number of obs = 126 Wald chi2(9) = 34.98 Log likelihood = -184.94833 Prob > chi2 = 0.0001 ------------------------------------------------------------------------------ | Robust numbids | Coef. Std. Err. z P>|z| [95% Conf. Interval] ---------+-------------------------------------------------------------------- leglrest | .2601464 .1250534 2.080 0.037 .0150462 .5052465 realrest | -.1956597 .1816168 -1.077 0.281 -.5516222 .1603027 finrest | .07403 .2635711 0.281 0.779 -.4425598 .5906198 whtknght | .4813821 .1064948 4.520 0.000 .2726562 .690108 bidprem | -.6776959 .297425 -2.279 0.023 -1.260638 -.0947536 insthold | -.3619913 .3231802 -1.120 0.263 -.9954127 .2714302 size | .1785026 .062355 2.863 0.004 .056289 .3007161 sizesq | -.0075694 .0027788 -2.724 0.006 -.0130158 -.0021229 regulatn | -.0294392 .1420508 -0.207 0.836 -.3078536 .2489752 _cons | .98606 .4137396 2.383 0.017 .1751453 1.796975 ------------------------------------------------------------------------------ . . . ******* EXTRAS . . * For models with two indexes see Weibull code in Stata 6 manual . . . ********** CLOSE OUTPUT . log close