* MMA12P2MSLMSM.DO March 2005 for Stata Version 8.0 log using mma12p2msmmsl.txt, text replace ********** OVERVIEW OF MMA12P2MSLMSM.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.4.5 pages 397-8 and 12.5.5 pages 402-4 * Computes integral numerically and by simulation * (1) Maximum Simulated likelihood Table 12.2 * (2) Method of Simulated Moments Table 12.3 * with application to generated data * The application is only illustrative. * This is not a template program for MSL or MSM. * Different number of simulations S lead to different estimators. * This program gives entries in Tables 12.2 and 12.3 for S = 100 * For other values of S change the value of simreps * from the current global simreps 100 ********** SETUP ********** set more off version 8 ********** DATA DESCRIPTION ********** * Model is y = theta + u + e * where theta is a scalar parameter equal to 1 * u is extreme value type 1 * e is N(0,1) * n is set in global numobs ********** DEFINE GLOBALS ********** global simreps 100 /* change this to change the number of simulations */ global numobs 100 /* change this to change the number of observations */ ********** (1) MAXIMUM SIMULATED LIKELIHOOD (Table 12.2 p.398) ********** * This MSL program is inefficiently written computer code * as it requires drawing the same random variates at each iteration * Generate data clear set obs $numobs set seed 10101 gen u = -log(-log(uniform())) gen e = invnorm(uniform()) gen y = 1 + u + e summarize u e y * Write data to a text (ascii) file so can use with programs other than Stata outfile u e y using mma12p2mslmsm.asc, replace * Use the variant ml d0 as this gives the entire likelihood, not just one observation. * I want this so that seed is only reset for the entire data. * My program is inefficient as variates needs to be redrawn at each iteration program define msl version 6.0 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 */ tempvar theta1 /* create as needed to calculate lf, g, ... */ mleval `theta1' = `b', eq(1) /* theta1 is theta1_i = x_i'b */ local y "$ML_y1" /* create to make program more readable */ set seed 10101 tempvar denssim global isim=1 quietly gen `denssim' = exp(-0.5*(`y'-`theta1'+log(-log(uniform())))^2)/sqrt(2*_pi) while $isim < $simreps { quietly replace `denssim' = `denssim' + exp(-0.5*(`y'-`theta1'+log(-log(uniform())))^2)/sqrt(2*_pi) global isim=$isim+1 } mlsum `lnf' = ln(`denssim'/$isim) end gen one = 1 ml model d0 msl (y = one, nocons ) ml maximize *** Display MSL results in one column of Table 12.2 p.398 di "For number of simulations S = " $simreps di "MSL estimator: " _b[one] di "Standard error: " _se[one] ********** (2) METHOD OF SIMULATED MOMENTS (Table 12.3 p.404) ********** clear set obs $numobs set seed 10101 gen u = -log(-log(uniform())) gen e = invnorm(uniform()) gen y = 1 + u + e summarize u e y global isim=1 gen usim = -log(-log(uniform())) gen esim = invnorm(uniform()) while $isim < $simreps { quietly replace usim = usim-log(-log(uniform())) quietly replace esim = esim+invnorm(uniform()) global isim=$isim+1 } gen usimbar = usim/$isim gen esimbar = esim/$isim gen theta = y - usimbar - esimbar summarize * Results for Table 12.3 on page 404 * Here the st.eror of theta_MSM is approximated by the st. dev. of theta * divided by the square root of S (the number of simulations) quietly sum theta scalar theta_MSM = r(mean) scalar approx_sterror = r(sd)/sqrt($simreps) * Display MSM results in one column of Table 12.3 p.404 di "For number of simulations S = " $simreps di "MSM estimator: " theta_MSM di "Approximate standard error: " approx_sterror * As written this will not give the correct standard errors (see p.403). * Can get this by also computing the squared rv to get E[y^2] ********** CLOSE OUTPUT ********** log close clear exit