------------------------------------------------------------------------------------------------------------------------------- name: log: c:\acdbookrevision\stata_final_programs_2013\racd12.txt log type: text opened on: 19 Jan 2013, 13:00:02 . . ********** OVERVIEW OF racd12.do ********** . . * STATA Program . * copyright C 2013 by A. Colin Cameron and Pravin K. Trivedi . * used for "Regression Analyis of Count Data" SECOND EDITION . * by A. Colin Cameron and Pravin K. Trivedi (2013) . * Cambridge University Press . . * This STATA program does Bayesian MCMC . * 12.4 MARKOV CHAIN MONTE CARLO METHODS . . * To run you need no data (the data are generated) . * in your directory . . ********** SETUP ********** . . set more off . version 12 . clear all . * set linesize 82 . set scheme s1mono /* Graphics scheme */ . * set maxvar 100 width 1000 . . ************ . . * Bayesian Random Walk Metropolis-Hastings with uniformative prior . * for Poisson regression using generated data . * Based on robit example of Koop (2003) chapter 9.3 "Bayesian Econometrics" . * Koop's Matlab program ch9artdatb.m rewritten for Poisson not probit . * plus adaptation to use random walk chain MH rather than data augmentation . . ********** SETUP ********** . . set mem 5m set memory ignored. Memory no longer needs to be set in modern Statas; memory adjustments are performed on the fly automatically. . set more off . set scheme s1mono . * version 11 . . ********** CREATE DATA AND SUMMARIZE ********** . . * Generate artificial data set for Poisson . * - explanatory variable x ~ N[0,1]sum . * - dependent variable y = Poisson(exp(-0.5*x) . * for x ~ N[0,1] . . set obs 100 obs was 0, now 100 . set seed 10101 . gen x = rnormal() . gen mu = exp(-0.5*x) . gen y = rpoisson(mu) . gen cons = 1 . sum Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- x | 100 .0624741 .9061909 -2.104658 3.170709 mu | 100 1.068019 .4783192 .2048751 2.864315 y | 100 .93 1.084836 0 4 cons | 100 1 0 1 1 . . * Alternative seed . * set seed 12345 . . ********** POISSON MAXIMUM LIKELIHOOD ESTIMATION (for Table 12.1) . . poisson y x Iteration 0: log likelihood = -126.0206 Iteration 1: log likelihood = -126.0206 Poisson regression Number of obs = 100 LR chi2(1) = 14.60 Prob > chi2 = 0.0001 Log likelihood = -126.0206 Pseudo R2 = 0.0548 ------------------------------------------------------------------------------ y | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- x | -.4574116 .1216507 -3.76 0.000 -.6958425 -.2189806 _cons | -.1254831 .1094268 -1.15 0.251 -.3399557 .0889895 ------------------------------------------------------------------------------ . estimates store POISSON . . ********** POISSON BAYESIAN WITH RANDOM WALK CHAIN METROPOLIS HASTNGS ********** . . * Globals for number of reps and a key tuning parameter . . global s1 10000 // number of retained reps . global s0 10000 // number of burnin reps . global sdscale 0.25 // use random walk b + $sdscale * N(0, I) . . * Mata to obtain the posterior draws of b . . mata ------------------------------------------------- mata (type end to exit) ----------------------------------------------------- : // Create y vector and X matrix from Stat data set using st_view command : st_view(y=., ., "y") // dependent : st_view(X=., ., ("cons", "x")) // regressors : Xnames = ("cons", "x") // used to label output : : // Calculate a few quantities outside the loop for later use : n = rows(X) : k = cols(X) : ones = J(n,1,1) : : // Specify the number of replications : s0 = $s0 // number of burnin reps : s1 = $s1 // number of retained reps : s = s0+s1 // total reps : : // Store all draws and MH acceptance ratein the following matrices : b_all = J(k,s1,0) : accept_all = J(1,s1,0) : : // Initialization : bdraw = J(k,1,0) // starting b value is vector of zeroes : lpostdraw = -1*10^10 // starting value of ln(posterior) is small : // so accept initial MH draw : : // Now do Metropolis-Hastings loop : : for (irep=1; irep<=s; irep++) { > > // Draw new candidate value of b from MH random walk chain > bcandidate = bdraw + $sdscale*rnormal(k,1,0,1) > > // Note: For different data you may need to change the global sdscale > // And best is bcandidate = bdraw + z where z ~ N(0, posterior variance of b) > > // Compute the log posterior at the candidate value of b > // The assumed prior for b is uninformative > // so the posterior is proportional to the usual Poisson likelihood > // and for comparsions of lnL we can drop term ln(y_i!) > Xb = X*bcandidate > lpostcandidate = ones'(-exp(Xb) + y:*Xb) > > // Accept the candidate draw on basis of posterior probability ratio > // if uniform > (posterior(bcandidate) / posterior(bdraw)) > // where bcandidate is current b and bdraw is previous b > // Taking logs the rule is the same as > // if ln(uniform) > (lpostcandidate - lpostdraw) > laccprobability = lpostcandidate - lpostdraw > accept = 0 > if ( ln(runiform(1,1)) < laccprobability ) { > lpostdraw = lpostcandidate > bdraw = bcandidate > accept = 1 > } > > // Store the draws after burn-in of b and whether accept draw > if (irep>s0) { > // after discarding burnin, store all draws > j = irep-s0 > b_all[.,j] = bdraw // These are the posterior draws > accept_all[.,j] = accept // These are one if new draw accepted > } > > } : : // End MH loop : : // Pass results back to Stata : // This bit works only for k = 2 (intercept plus one slope) : beta = b_all' : accept = accept_all' : st_addvar("float", ("beta1", "beta2", "accept")) 1 2 3 +-------------+ 1 | 7 8 9 | +-------------+ : // The following line is needed for conformability : // It assumes that s1 (number of draws) exceeds the original sample size : stata("set obs $s1") obs was 100, now 10000 : st_store(., ("beta1", "beta2"), beta) : st_store(., ("accept"), accept) : end ------------------------------------------------------------------------------------------------------------------------------- . . ********** ANALYZE THE RESULTS ********** . . *** TABLE 12.1: BAYESIAN POSTERIOR SUMMARY . . * Posterior means and standard deviations . summarize beta1 beta2 Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- beta1 | 10000 -.1349377 .1102719 -.551262 .2277995 beta2 | 10000 -.4638151 .1206984 -.8691139 -.0125802 . . * Posterior credible regions . centile beta1 beta2, centile(2.5 97.5) -- Binom. Interp. -- Variable | Obs Percentile Centile [95% Conf. Interval] -------------+------------------------------------------------------------- beta1 | 10000 2.5 -.3560977 -.3717561 -.35214 | 97.5 .0718728 .0671539 .0743966 beta2 | 10000 2.5 -.7058128 -.7121991 -.7007356 | 97.5 -.2282387 -.2330735 -.2245314 . . summarize Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- x | 100 .0624741 .9061909 -2.104658 3.170709 mu | 100 1.068019 .4783192 .2048751 2.864315 y | 100 .93 1.084836 0 4 cons | 100 1 0 1 1 _est_POISSON | 100 1 0 1 1 -------------+-------------------------------------------------------- beta1 | 10000 -.1349377 .1102719 -.551262 .2277995 beta2 | 10000 -.4638151 .1206984 -.8691139 -.0125802 accept | 10000 .2565 .4367228 0 1 . * save bayespoisson.dta, replace . rename beta2 b . . * MLE of Poisson . estimates table POISSON, b(%7.4f) se(%7.3f) stats(N ll) stfmt(%9.1f) -------------------------- Variable | POISSON -------------+------------ x | -0.4574 | 0.122 _cons | -0.1255 | 0.109 -------------+------------ N | 100 ll | -126.0 -------------------------- legend: b/se . . *** FIGURE 12.1: POSTERIOR DRAWS AND ACF . . * Plot the posterior draws of b2 and acf . generate s= _n . tsset s time variable: s, 1 to 10000 delta: 1 unit . line b s if s < 100, name(first100, replace) ytitle("{&beta}{sub:2}") /// > xtitle("Draw s") scale(1.9) . line b s if s < 2000, name(all, replace) ytitle("{&beta}{sub:2}") /// > xtitle("Draw s") scale(1.9) . graph combine first100 all, iscale(1.0) ysize(3) xsize(6) ycommon . graph export racd12fig1.wmf, replace (note: file racd12fig1.wmf not found) (file c:\acdbookrevision\stata_final_programs_2013\racd12fig1.wmf written in Windows Metafile format) . graph export racd12fig1.eps, replace (note: file racd12fig1.eps not found) (file racd12fig1.eps written in EPS format) . . *** FIGURE 12.2: POSTERIOR DENSITY . . kdensity b, note("") legend(position(0)) xtitle("{&beta}{sub:2}") /// > ytitle("Posterior density") title("") scale(1.9) . graph export racd12fig2.wmf, replace (note: file racd12fig2.wmf not found) (file c:\acdbookrevision\stata_final_programs_2013\racd12fig2.wmf written in Windows Metafile format) . graph export racd12fig2.eps, replace (note: file racd12fig2.eps not found) (file racd12fig2.eps written in EPS format) . . *** OTHER DETAILS IN TEXT . . * Correlation of the draws . corrgram b, lags(20) -1 0 1 -1 0 1 LAG AC PAC Q Prob>Q [Autocorrelation] [Partial Autocor] ------------------------------------------------------------------------------- 1 0.8003 0.8003 6406.1 0.0000 |------ |------ 2 0.6362 -0.0121 10456 0.0000 |----- | 3 0.5069 0.0038 13027 0.0000 |---- | 4 0.4007 -0.0091 14633 0.0000 |--- | 5 0.3173 0.0021 15641 0.0000 |-- | 6 0.2474 -0.0109 16253 0.0000 |- | 7 0.1902 -0.0067 16615 0.0000 |- | 8 0.1418 -0.0120 16816 0.0000 |- | 9 0.1029 -0.0061 16922 0.0000 | | 10 0.0793 0.0147 16985 0.0000 | | 11 0.0640 0.0081 17026 0.0000 | | 12 0.0529 0.0036 17054 0.0000 | | 13 0.0452 0.0039 17075 0.0000 | | 14 0.0417 0.0090 17092 0.0000 | | 15 0.0370 -0.0023 17106 0.0000 | | 16 0.0383 0.0157 17120 0.0000 | | 17 0.0422 0.0122 17138 0.0000 | | 18 0.0417 -0.0039 17156 0.0000 | | 19 0.0378 -0.0051 17170 0.0000 | | 20 0.0292 -0.0122 17179 0.0000 | | . . * Give the acceptance rate in the random walk chain MH . quietly summarize accept . display "MH acceptance rate = " r(mean) " MH acceptance rate = .2565 . . * Check the efficiency loss due to correlation of draws . regress b Source | SS df MS Number of obs = 10000 -------------+------------------------------ F( 0, 9999) = 0.00 Model | 0 0 . Prob > F = . Residual | 145.666431 9999 .0145681 R-squared = 0.0000 -------------+------------------------------ Adj R-squared = 0.0000 Total | 145.666431 9999 .0145681 Root MSE = .1207 ------------------------------------------------------------------------------ b | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _cons | -.4638151 .001207 -384.28 0.000 -.4661811 -.4614492 ------------------------------------------------------------------------------ . newey b, lag(40) Regression with Newey-West standard errors Number of obs = 10000 maximum lag: 40 F( 0, 9999) = . Prob > F = . ------------------------------------------------------------------------------ | Newey-West b | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _cons | -.4638151 .0033588 -138.09 0.000 -.4703991 -.4572312 ------------------------------------------------------------------------------ . . ********** CLOSE OUTPUT ********** . . * log close . * clear . * exit . end of do-file . exit, clear