*! mkbilogn.ado version 1.1.1 spj revised 27 Jan 98 *! Syntax: mkbilogn var1 var2, r(#) m1(#) s1(#) m2(#) s2(#) *! [defaults #=0.5,0,1,0,1; respectively] prog def mkbilogn version 5.0 local varlist "req new min(2) max(2)" local options "Rho(real 0.5) m1(real 0) m2(real 0) s1(real 1) s2(real 1)" parse "`*'" if `rho' < -1 | `rho' >= 1 { di in red "Need rho s.t. -1 <= rho < 1" exit 198 /* need rho ~= 1 or Cholesky fails */ } if `s1' <=0 | `s2' <= 0 { di in red "Std. dev. must be positive" exit 198 } parse "`varlist'", parse(" ") di "Creating 2 r.v.s X1 X2 s.t. x1=log(X1), x2=log(X2) are bivariate" di " Normal with mean(x1) = `m1' ; mean(x2) = `m2' ; s.d.(x1) = `s1' ;" di " s.d.(x2) = `s2' ; corr(x1,x2) = `rho' " /* Method of creation based on Stata FAQ at http://www.stata.com/support/faqs/stat/mvnorm.html */ tempname a1 a2 A P var1 var2 c1 c2 lnc1 lnc2 matrix `P' = (1,`rho'\ `rho',1) matrix `A' = cholesky(`P') matrix colnames `A' = `lnc1' `lnc2' quietly { ge `c1' = exp(invnorm(uniform())) ge `c2' = exp(invnorm(uniform())) ge `lnc1' = ln(`c1') ge `lnc2' = ln(`c2') matrix `a1' = `A'[1,.] matrix score `var1' = `a1' matrix `a2' = `A'[2,.] matrix score `var2' = `a2' replace `1' = exp(`s1'*`var1' + `m1') replace `2' = exp(`s2'*`var2' + `m2') } end /* Note: mean,variances,corr refer to mean,variances,corr of the log of the vbles. Mean of lognormal distribution: exp[mu + (sig_sq/2)] Variance of lognormal distribution: [exp(2*mu)]*[exp(2*sig_sq) - exp(sig_sq)] */