Does this definite integral have a closed form?
$$ \mathbb{P}_\varepsilon(v_2 + \varepsilon_2 > v_1 + \varepsilon_1 > \varepsilon_0) = \int_{-\infty}^\infty \int_{-\infty}^{v_2+\varepsilon_2-v_1} \int_{-\infty}^{v_1+\varepsilon_1} f(\varepsilon) \,\mathrm{d}\varepsilon_0 \,\mathrm{d}\varepsilon_1 \,\mathrm{d}\varepsilon_2 $$
where $\varepsilon$ is distributed generalized extreme value with cdf
$$ F(\varepsilon) = \exp\left(-e^{-\varepsilon_0}-\left(e^{-\varepsilon_1/\lambda} + e^{-\varepsilon_2/\lambda}\right)^\lambda\right) $$
and pdf
$$ \begin{align} f(\varepsilon) &= \frac{\partial^3 F(\varepsilon)}{\partial\varepsilon_0\partial\varepsilon_1\partial\varepsilon_2} \\ &= F(\varepsilon) \cdot e^{-\varepsilon_0 - (\varepsilon_1 + \varepsilon_2) / \lambda} \cdot \left(\left(e^{-\varepsilon_1 / \lambda} + e^{-\varepsilon_2 / \lambda}\right)^{2(\lambda-1)} - \frac{\lambda-1}{\lambda}\cdot\left(e^{-\varepsilon_1 / \lambda} + e^{-\varepsilon_2 / \lambda}\right)^{\lambda-2}\right). \end{align} $$
Both $v_1$ and $v_2$ are real-valued and $0 < \lambda \leq 1$.
In a nested logit choice model (e.g., Train, 2009) with three alternatives $j \in \{0, 1, 2\}$ and indirect utility $u_j = v_j + \varepsilon_j$ where alternatives $j \in \{1, 2\}$ are in the same nest and we set $v_0 = 0$, this integral is of interest because it gives the probability of $j = 2$ being the best and $j = 1$ being the second best.
The simpler probability of just $j = 2$ being the best is nice and well-known:
$$ \mathbb{P}_\varepsilon(u_2 > \max\{u_1, u_0\}) = \underbrace{\frac{\left(e^{v_1/\lambda} + e^{v_2/\lambda}\right)^\lambda}{1 + \left(e^{v_1/\lambda} + e^{v_2/\lambda}\right)^\lambda}}_{\mathbb{P}_\varepsilon(\max\{u_1, u_2\} > u_0)} \cdot \underbrace{\frac{e^{v_2/\lambda}}{e^{v_1/\lambda} + e^{v_2/\lambda}}}_{\mathbb{P}_\varepsilon(u_2 > u_1)}. $$
If $\lambda = 1$, this collapses to a non-nested logit model, and we also get a nice expression:
$$ \begin{align} \lambda = 1 \implies \mathbb{P}_\varepsilon(u_2 > u_1 > u_0) = \underbrace{\frac{e^{v_2}}{1 + e^{v_1} + e^{v_2}}}_{\mathbb{P}_\varepsilon(u_2 > \max\{u_1, u_0\})} \cdot \underbrace{\frac{e^{v_1}}{1 + e^{v_1}}}_{\mathbb{P}_\varepsilon(u_1 > u_0)}. \end{align} $$
Mathematica manages to recover both of these known results, but it fails when I try to give it the integral of interest. So it may not have a closed form, but I'm hopeful!
Thanks so much for any thoughts or pointers.
Figured it out!
$$ \mathbb{P}_\varepsilon(u_2 > u_1 > u_0) = \mathbb{P}_\varepsilon(u_1 > u_0) - \mathbb{P}_\varepsilon(u_1 > \max\{u_2, u_0\}), $$
in which the two expression on the right hand side are given in my original question, with the latter term having $u_1$ and $u_2$ swapped.