------------------------------------------------------------------------------------------------------
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
----------------------------------------------------------------------------------------------------