Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • bayesmh with skewed normal distribution of prior

    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.

    Code:
    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

  • #2
    Originally posted by Sarah Frohnweiler View Post
    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.
    Code:
    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.
    Code:
    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.)
    Code:
    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?)
    Attached Files

    Comment


    • #3
      Originally posted by Joseph Coveney View Post
      And then your prior would be specified in a substitutable expression as follows.
      Code:
      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.
      Code:
      prior({y:treatment}, ///
          density(2 / omega * normalden(({y:treatment} - ksi) / omega) * normal(alpha * ({y:treatment} - ksi) / omega)))

      Comment


      • #4
        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

        Comment

        Working...
        X