Simulating Draws of Multivariate EV-Type Distribution

Question

Let $\varepsilon = [\varepsilon_1,...,\varepsilon_J]$ be a random vector that we can partition into $K$ disjoint subvectors. $\varepsilon$ has this cdf:

\begin{equation} F(\varepsilon) = \exp \bigg[-\sum_{k=1,...,K}\Big ( \sum_{j\in J_k} e^{\varepsilon_j / \gamma} \Big )^ \gamma \bigg].\end{equation}

This is the distribution of nested logit errors in a discrete choice model, where the elements belonging to the same nest $k$ are correlated according to $\gamma \in [0,1]$. I need to simulate draws from this distribution but cannot figure out how, without differentiating $J$ times and getting the pdf so I can do MCMC.

There is another question here that gives a solution that seems pretty cumbersome for $J$ large. And another question here that went unresolved, but maybe the "nice" properties of EV distributions may be helpful?

Answer 1

I asked and then answered this question on Stack Overflow. Here is an example for a particular nesting structure using the rmvevd command in the evd package for R:

D = 4 alternatives, alternative 1 is in nest 1; alternatives 2,3,4 are in nest 2.

lambda_2 = 0.5 is the nesting parameter.

B = the set of all possible sets that can be formed from the 4 alternatives: (B = {1}, {2}, {3}, {4}, {1,2} {1,3}, {1,4}, {2,3}, {2,4} {3,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}, {1,2,3,4})

mvevd requires the user to specify the arguments dep and asy. dep is a vector of nesting parameters on the sets in B of length 2 or greater. asy (in the nested logit example) is a list that specifies the nesting structure. In the example here asy will be vector of 1s in the 1st and 14th entries (corresponding to alternative 1 in its own group, and alternatives 2,3,4 in group 2). dep will have a 1 in all entries except for the 10th entry, which will equal lambda_2 (as {2,3,4} is the 10th element of B that is not of unit length). A code example to simulate 10 4-d nested logit deviates according to this example is as follows:

library('evd')
vdep <- c(rep(1,9), lambda, 1)
asy_list <- list(1, 0, 0, 0, c(0,0), c(0,0), c(0,0), c(0,0), c(0,0), c(0,0),
            c(0,0,0), c(0,0,0), c(0,0,0), c(1,1,1), c(0,0,0,0))
deviates <- rmvevd(10, dep = vdep, asy = asy_list, model = "alog", d = 4)