bayesmh with skewed normal distribution of prior

[Sarah Frohnweiler]

Dear all,

I am running a Bayesian Analysis using the bayesmh command in Stata 16.1. The command runs perfectly fine, but I would like the customized prior for my main explanatory variable "treatment" to have a skewed normal distribution instead of its current normal distribution.

bayesmh y treatment x i.trade##i.region, likelihood(normal({sigma2}))   ///
       prior({y:treatment}, normal(1.920304, 0.3538458))                ///
       prior({y:x}, normal(0, 2.5))                                    ///
       prior({y:i.trade}, normal(0, 2.5))                              ///
       prior({y:i.region}, normal(0, 2.5))                            ///
       prior({y:i.trade#i.region}, normal(0, 2.5))                    ///
       prior({y:_cons}, normal(0.8740444, 3.875811))                 ///
       prior({sigma2}, gamma(1, 1))                                  ///
       nchains(4) mcmcsize(60000) burnin(20000) rseed(78643)       ///
       block({y:x _cons})                                      ///
       block({y:i.trade i.region i.trade#i.region})          ///
       block({sigma2})

The help file lists many prior distributions to choose from, but as far as I can tell no skewed normal distribution. For the variable of interest "treatment" I currently define a prior that is normally distributed with mean = 1.920304 and var = 0.3538458. I would like to adjust the prior such that it is skewed normally distributed with mean = 1.920304, var = 0.3538458, and shape = -0.1574642.

Do you know how to achieve this? Any suggestions are greatly appreciated!

Thanks in advance and best,
Sarah

[Joseph Coveney]

For the variable of interest "treatment" I currently define a prior that is normally distributed with mean = 1.920304 and var = 0.3538458. I would like to adjust the prior such that it is skewed normally distributed with mean = 1.920304, var = 0.3538458, and shape = -0.1574642.

Do you know how to achieve this?

I assume that what you refer to as shape is the parameter alpha in the Wikipedia article on the skew normal distribution. If so, then you could specify the skew-normal prior density using a substitutable expression. You'd compute the requisite parameters using the method of moments formulas shown under the Estimation header in the WIkipedia article. It would be something like the following. The parameters are named following the Wikipedia article.

scalar define mu = 1.920304
scalar define sigma = sqrt(0.3538458)
scalar define alpha = -0.1574642

scalar define delta = alpha / sqrt(1 + alpha * alpha)
scalar define omega = sigma / sqrt(1 - 2 / _pi * delta * delta)
scalar define ksi = mu - omega * delta * sqrt(2 / _pi)

And then your prior would be specified in a substitutable expression as follows.

prior({y:treatment}, ///
    density(2 / omega * normalden(({x} - ksi) / omega) * normal(alpha * ({x} - ksi) / omega)))

I show a worked example below with a few arbitrarily chosen parameters for mean (-1), standard deviation (2) and coefficient of skew (0.75) as the first example at the top. I then show yours as the second Bayesian estimation, beginning at the comment "With mean, variance and 'shape' given by OP:". I've left the likelihood flat so as to show the naked prior distribution as the posterior in order to confirm that the code for specifying the skew-normal prior density is correct. (The complete do-file and corresponding log file are attached to this post.)

version 16.1

clear *

// seedem
set seed 874018287

program define summary
    version 16.1
    syntax

    tempname Posterior
    frame create `Posterior'
    cwf `Posterior'

    use posterior
    
    quietly summarize eq0_p1 [fweight=_frequency], detail
    display in smcl as text _newline(1) "         Mean: " as result %04.2f r(mean)
    display in smcl as text "     Variance: " as result %04.2f r(Var)
    display in smcl as text "Coef. of Skew: " as result %04.2f r(skewness)

    cwf default
    erase posterior.dta
end

*
* Setting up parameters (method of moments)
*
// Desired moments
scalar define mu = -1 // mean
scalar define sigma = 2 // standard deviation
scalar define gamma1 = 0.75 // coefficient of skew

// Computed parameters
scalar define two_thirds = 2 / 3
scalar define fmpot = (4 - _pi) / 2
scalar define gamma23 = abs(gamma1)^two_thirds
scalar define delta = sign(gamma1) * ///
    sqrt(_pi / 2 * gamma23 / (gamma23 + fmpot^two_thirds))

scalar define alpha = delta / sqrt(1 - delta * delta)
scalar define omega = sigma / sqrt(1 - 2 / _pi * delta * delta)
scalar define ksi = mu - omega * delta * sqrt(2 / _pi)

*
* Draws from prior distribution
*
quietly set obs 1

generate double y = 0

bayesmh y, likelihood(llf(0)) ///
    prior({x}, density(2 / omega * normalden(({x} - ksi) / omega) * ///
        normal(alpha * ({x} - ksi) / omega))) ///
    mcmcsize(10000) nomodelsummary saving(posterior)

summary

*
* With mean, variance and "shape" given by OP:
*
scalar define mu = 1.920304
scalar define sigma = sqrt(0.3538458)
scalar define alpha = -0.1574642

// Skewness
scalar define delta = alpha / sqrt(1 + alpha * alpha)
scalar define gamma1 = (4 - _pi) / 2 * ((delta * sqrt(2 / _pi))^3) / ///
    ((1 - 2 * delta^2 / _pi)^(3/2))
display in smcl as text "Coeffient of skew: " as result %04.2f gamma1

scalar define omega = sigma / sqrt(1 - 2 / _pi * delta * delta)
scalar define ksi = mu - omega * delta * sqrt(2 / _pi)

bayesmh y, likelihood(llf(0)) ///
    prior({x}, density(2 / omega * normalden(({x} - ksi) / omega) * ///
        normal(alpha * ({x} - ksi) / omega))) ///
    mcmcsize(10000) nomodelsummary saving(posterior)

summary

exit

Please be aware that with your shape parameter at ca. -0.16, the skewness is essentially zero, given how bouncy that moment is. You might as well just go for the built-in normal density as you won't see much practical difference.

(Release 16 had frames, right?)

[Joseph Coveney]

And then your prior would be specified in a substitutable expression as follows.

prior({y:treatment}, ///
density(2 / omega * normalden(({x} - ksi) / omega) * normal(alpha * ({x} - ksi) / omega)))

Copy-paste error, sorry. The prior would be specified more like the following.

prior({y:treatment}, ///
    density(2 / omega * normalden(({y:treatment} - ksi) / omega) * normal(alpha * ({y:treatment} - ksi) / omega)))

[Sarah Frohnweiler]

Dear Joseph,

Thanks a lot for your detailed response and help! I was able to adjust the prior distribution accordingly and the code runs smoothly.

Best,
Sarah