------------------------------------------------------------------------------------------------------ log: c:\Imbook\bwebpage\Section3\mma12p1integration.txt log type: text opened on: 18 May 2005, 21:17:14 . . ********** OVERVIEW OF MMA12P1INTEGRATION.DO ********** . . * STATA Program . * copyright C 2005 by A. Colin Cameron and Pravin K. Trivedi . * used for "Microeconometrics: Methods and Applications" . * by A. Colin Cameron and Pravin K. Trivedi (2005) . * Cambridge University Press . . * Chapter 12.3.3 pages 391-2 . * Computes integral numerically and by simulation . * (1) Illustrate Midpoint Rule (page 392) . * (2) Illustrate Monte Carlo integral (Table 12.1 page 392) . * . * for computing E[x] and E[exp(-exp(x))] for x ~ N[0,1] . . * No data need be read in. . . ********** SETUP ********** . . set more off . version 8.0 . . ********** (1) NUMERICAL INTEGRATION USING MIDPOINT RULE ********** . . * Midpoint rule for n evaluation points between a and b is . * Integral = Sum (j=1 to n) [(b-a)/n]*f(xbar_j) . * where xbar_j is midpoint between x_j-1 and x_j . . program midpointrule, rclass 1. version 8 2. /* define arguments. Here trueb2 = b2 in Phi(b1 + b2*x2) */ . args neval a b 3. drop _all 4. scalar increment = (`b'-`a') / `neval' 5. set obs `neval' 6. /* Compute the function of interest */ . gen xbar = `a' - 0.5*increment + increment*_n 7. gen density = exp(-xbar*xbar/2)/sqrt(2*_pi) 8. * Following is contribution to E[x] when x ~ N[0,1] . gen f1xbar = xbar*density 9. * Following is contribution to E[exp(-exp(x))] when x ~ N[0,1] . gen f2xbar = exp(-exp(x))*density 10. /* Compute the averages */ . quietly sum f1xbar 11. scalar Ex = r(sum)*increment 12. quietly sum f2xbar 13. scalar Eexpminexpx = r(sum)*increment 14. /* Print results */ . di "Evaluation points: " `neval' " over range: (" `a' "," `b' ") 15. di "Midpoint rule estimate of E[x] is: " Ex 16. di "Midpoint rule estimate of E[exp(-exp(x))] is: " Eexpminexpx 17. end . . midpointrule 20 -5 5 obs was 0, now 20 Evaluation points: 20 over range: (-5,5) Midpoint rule estimate of E[x] is: 0 Midpoint rule estimate of E[exp(-exp(x))] is: .38175625 . midpointrule 200 -5 5 obs was 0, now 200 Evaluation points: 200 over range: (-5,5) Midpoint rule estimate of E[x] is: 0 Midpoint rule estimate of E[exp(-exp(x))] is: .38175618 . midpointrule 2000 -5 5 obs was 0, now 2000 Evaluation points: 2000 over range: (-5,5) Midpoint rule estimate of E[x] is: 0 Midpoint rule estimate of E[exp(-exp(x))] is: .38175618 . . ********** (2) MONTE CARLO INTEGRATION USING DRAWS FROM DENSITY OF X ********** . . * To get E[g(x)] . * make draws from N[0,1], compute g(x), and average over draws . . program simintegration, rclass 1. version 8 2. /* define arguments. Here trueb2 = b2 in Phi(b1 + b2*x2) */ . args nsims 3. /* Generate the data: here x */ . drop _all 4. set obs `nsims' 5. set seed 10101 6. gen x = invnorm(uniform()) 7. /* Compute the function of interest */ . gen f1x = x /* For E[x] just need x */ 8. gen f2x = exp(-exp(x)) /* For E[exp(-exp(x))] */ 9. /* Compute the averages */ . quietly sum f1x 10. scalar Ex = r(mean) 11. quietly sum f2x 12. scalar Eexpminexpx = r(mean) 13. di "Number of simulations: " `nsims' 14. di "Monte Carlo estimate of E[x] is: " Ex 15. di "Monte Carlo estimate of E[exp(-exp(x))] is: " Eexpminexpx 16. end . . * Note a different program was used to obtain Table 12.1 on page 392 . * So results will differ somewhat from text, except for very high number of simulations . . simintegration 10 obs was 0, now 10 Number of simulations: 10 Monte Carlo estimate of E[x] is: -.10143571 Monte Carlo estimate of E[exp(-exp(x))] is: .42635197 . simintegration 25 obs was 0, now 25 Number of simulations: 25 Monte Carlo estimate of E[x] is: .17496346 Monte Carlo estimate of E[exp(-exp(x))] is: .35703296 . simintegration 50 obs was 0, now 50 Number of simulations: 50 Monte Carlo estimate of E[x] is: .0079132 Monte Carlo estimate of E[exp(-exp(x))] is: .37966293 . simintegration 100 obs was 0, now 100 Number of simulations: 100 Monte Carlo estimate of E[x] is: .11238423 Monte Carlo estimate of E[exp(-exp(x))] is: .3524417 . simintegration 500 obs was 0, now 500 Number of simulations: 500 Monte Carlo estimate of E[x] is: .06990338 Monte Carlo estimate of E[exp(-exp(x))] is: .36137551 . simintegration 1000 obs was 0, now 1000 Number of simulations: 1000 Monte Carlo estimate of E[x] is: .04309113 Monte Carlo estimate of E[exp(-exp(x))] is: .36945581 . simintegration 1000 obs was 0, now 1000 Number of simulations: 1000 Monte Carlo estimate of E[x] is: .04309113 Monte Carlo estimate of E[exp(-exp(x))] is: .36945581 . simintegration 100000 obs was 0, now 100000 Number of simulations: 100000 Monte Carlo estimate of E[x] is: -.00405425 Monte Carlo estimate of E[exp(-exp(x))] is: .38284684 . clear . set mem 20m (20480k) . simintegration 1000000 obs was 0, now 1000000 Number of simulations: 1000000 Monte Carlo estimate of E[x] is: -.00085186 Monte Carlo estimate of E[exp(-exp(x))] is: .38192861 . . ********** CLOSE OUTPUT ********** . log close log: c:\Imbook\bwebpage\Section3\mma12p1integration.txt log type: text closed on: 18 May 2005, 21:17:16 ----------------------------------------------------------------------------------------------------