* This program takes random draws of the estimated parameters and recomputes WTP measures
*! Author PWJ
program wtpkrinsky, rclass
	version 9.2
     	syntax varlist, Bname(name) Vname(name) MODel(string) [MYMean(name) EXPOnential]
	tempname z
   	mata: takedraw()

	* Now use the new parameter vector z to compute WTP measures
 	wtpmeasure `varlist', bvec(`z') mod(`model') mym(`mymean') `exponential'
	if "`exponential'"=="" ret sca meanwtp=r(meanwtp)
	else {
		ret sca medianwtp=r(medianwtp)
		ret sca meanwtp=r(meanwtp)
	}
end
version 9.2
mata:
void takedraw()
{
	real matrix b, v, Z
	b = st_matrix(st_local("bname"))
	v = st_matrix(st_local("vname"))
	Z = (b' :+ cholesky(v)*invnormal(uniform(cols(b),1)))'  // Take a random draw from multivariate normal distribution
	st_matrix(st_local("z"),Z)
}
end