Generalized Rank Dirichlet Distributions

TL;DR

This paper introduces the Generalized Rank Dirichlet (GRD) distributions, a new family of distributions on the ordered simplex that generalize the Dirichlet distribution, allowing for negative parameters and providing explicit moments and simulation methods.

Contribution

The paper defines GRD distributions on the ordered simplex, derives explicit moments, and develops exact and approximate simulation algorithms, expanding modeling capabilities for ranked data.

Findings

Explicit moments for GRD distributions across dimensions.

Series representations and simulation algorithms for the distributions.

Application potential in financial modeling and ranked statistics.

Abstract

We study a new parametric family of distributions on the ordered simplex $\nabla^{d-1} = \{y \in \mathbb{R}^d: y_1 \geq \dots \geq y_d \geq 0, \sum_{k=1}^d y_k = 1\}$, which we call Generalized Rank Dirichlet (GRD) distributions. Their density is proportional to $\prod_{k=1}^d y_k^{a_k-1}$ for a parameter $a = (a_1,\dots,a_d) \in \mathbb{R}^d$ satisfying $a_k + a_{k+1} + \dots + a_d > 0$ for $k=2,\dots,d$. The density is similar to the Dirichlet distribution, but is defined on $\nabla^{d-1}$, leading to different properties. In particular, certain components $a_k$ can be negative. Random variables $Y = (Y_1,\dots,Y_d)$ with GRD distributions have previously been used to model capital distribution in financial markets and more generally can be used to model ranked order statistics of weight vectors. We obtain for any dimension $d$ explicit expressions for moments of order $M \in \mathbb{N}$ for the $Y_k$'s and moments of all orders for the log gaps $Z_k = \log Y_{k-1} - \log Y_k$ when $a_1 + \dots + a_d = -M$. Additionally, we propose an algorithm to exactly simulate random variates in this case. In the general case $a_1 + \dots + a_d \in \mathbb{R}$ we obtain series representations for these quantities and provide an approximate simulation algorithm.

Equations

$${\nabla }^{d-1}=\{y\in {\mathbb{R}}^d:y_1\ge y_2\ge \cdots \ge y_d\ge 0\text{ and }{y}_1+\cdots +y_d=1\},$$

$$\prod _{k=1}^d{y}_k^{a_k-1}$$

$${\bar{a}}_k:=a_k+a_{k+1}+\cdots +a_d>0,\text{for }k=2,\dots ,d.$$

$$y^{a_1-1}(1-y{)}^{a_2-1},y\in [1/2,1].$$

$$Z_k=\log Y_{k-1}-\log Y_k,\text{for }k=2,\dots ,d$$

$$Q_a={\int }_{{\nabla }^{d-1}}\prod _{k=2}^d{\left(\frac{y_{k-1}}{y_k}\right)}^{-{\bar{a}}_k}\prod _{k=1}^d{y}_k^{-1}dy.$$

$$Q_a={\int }_{{\mathbb{R}}_{+}^{d-1}}\exp \left(-\sum _{k=2}^d{\bar{a}}_k{z}_k\right)dz=\prod _{k=2}^d{\int }_0^{\infty }{e}^{-{\bar{a}}_{k}z}dz.$$

$${\mathbb{P}}_a(A)=Q_a^{-1}{\int }_{{\nabla }^{d-1}}\prod _{k=1}^d{y}_k^{a_k-1}{1}_A(y)dy,A\in \mathcal{B}({\nabla }^{d-1})$$

$${\mathbb{E}}_a\left[\frac{{\prod }_{k=1}^d{Y}_k^{n_k}}{Y_1^M}\right]=\sum _{m\in {\mathbb{N}}_0^d(M-{\bar{n}}_1)}\left(\genfrac{}{}{0ex}{}{M-{\bar{n}}_1}{m_1,\dots ,m_d}\right)\prod _{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{m}}_k+{\bar{n}}_k}.$$

$${\mathbb{E}}_a\left[\frac{1}{Y_1^M}\right]=\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right)\prod _{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{m}}_k}.$$

$${\mathbb{E}}_a\left[\frac{{\prod }_{k=1}^d{Y}_k^{n_k}}{Y_1^M}\right]={\mathbb{E}}_a\left[\frac{{\prod }_{k=1}^d{Y}_k^{n_k}}{Y_1^{{\bar{n}}_1}}\right]=\prod _{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{n}}_k}$$

$${\mathbb{E}}_a\left[\frac{{\prod }_{k=1}^d{Y}_k^{n_k}}{Y_1^M}\right]$$

$$=\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M-{\bar{n}}_1}{m_1,\dots ,m_d}\right){\mathbb{E}}_a\left[\frac{{\prod }_{k=1}^d{Y}_k^{n_k+m_k}}{Y_1^M}\right]$$

$$=\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M-{\bar{n}}_1}{m_1,\dots ,m_d}\right)\prod _{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{m}}_k+{\bar{n}}_k}.$$

$${\mathbb{E}}_a[f(Y)]=\frac{{\mathbb{E}}_b[f(Y){\prod }_{k=1}^d{Y}_k^{a_k-b_k}]}{{\mathbb{E}}_b[{\prod }_{k=1}^d{Y}_k^{a_k-b_k}]}.$$

$${\mathbb{E}}_a[f(Y)]$$

$$=\frac{{\int }_{{\nabla }^{d-1}}f(y){\prod }_{k=1}^d{y}_k^{a_k-b_k}{\prod }_{k=1}^d{y}_k^{b_k-1}dy}{{\int }_{{\nabla }^{d-1}}{\prod }_{k=1}^d{y}_k^{b_k-1}dy}\times \frac{{\int }_{{\nabla }^{d-1}}{\prod }_{k=1}^d{y}_k^{b_k-1}dy}{{\int }_{{\nabla }^{d-1}}{\prod }_{k=1}^d{y}_k^{a_k-b_k}{\prod }_{k=1}^d{y}_k^{b_k-1}dy}$$

$$=\frac{{\mathbb{E}}_b[f(Y){\prod }_{k=1}^d{y}_k^{a_k-b_k}]}{{\mathbb{E}}_b[{\prod }_{k=1}^d{y}_k^{a_k-b_k}]},$$

$${\mathbb{E}}_a[f(Y)]=\frac{{\mathbb{E}}_{a-{\bar{a}}_1{e}_1}[f(Y)Y_1^{{\bar{a}}_1}]}{{\mathbb{E}}_{a-{\bar{a}}_1{e}_1}[Y_1^{{\bar{a}}_1}]}.$$

$${\mathbb{E}}_a\left[\prod _{k=1}^d{Y}_k^{n_k}\right]=\frac{\sum _{m\in {\mathbb{N}}_0^d(M-{\bar{n}}_1)}\left(\genfrac{}{}{0ex}{}{M-{\bar{n}}_1}{m_1,\dots ,m_d}\right)\prod _{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{m}}_k+{\bar{n}}_k}}{\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right)\prod _{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{m}}_k}}.$$

$${\mathbb{E}}_a[Y_k]=C^{-1}\prod _{j=2}^k\frac{{\bar{a}}_j}{{\bar{a}}_j+1},\text{where}C=1+\sum _{k=2}^d\prod _{j=2}^k\frac{{\bar{a}}_j}{{\bar{a}}_j+1}.$$

$$a_k=\{\begin{array}{ll}-1-\frac{y_2}{y_1-y_2},& k=1,\\ \frac{y_k}{y_{k-1}-y_k}-\frac{y_{k+1}}{y_k-y_{k+1}},& k=2,\dots ,d-1,\\ \frac{y_d}{y_{d-1}-y_d}& k=d.\end{array}$$

$$\begin{array}{rl}{\mathbb{E}}_{a+M{e}_1}\left[\frac{f(Y)}{Y_1^M}\right]& =\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right){\mathbb{E}}_{a+M{e}_1}\left[f(Y)\frac{{\prod }_{k=1}^d{Y}_k^{m_k}}{Y_1^M}\right]\\ & =\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right){\mathbb{E}}_{a+M{e}_1}\left[\frac{{\prod }_{k=1}^d{Y}_k^{m_k}}{Y_1^M}\right]{\mathbb{E}}_{a+m}[f(Y)]\\ & =\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right)\prod _{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{m}}_k}{\mathbb{E}}_{a+m}[f(Y)],\end{array}$$

$${\mathbb{E}}_a[f(Y)]=\sum _{m\in {\mathbb{N}}_0^d(M)}{w}_m{\mathbb{E}}_{a+m}[f(Y)]\text{where}{w}_m=\frac{\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right){\prod }_{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{m}}_k}}{{\sum }_{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right){\prod }_{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k+{\bar{m}}_k}}$$

$${\mathbb{E}}_a[g(Z)]=\sum _{m\in {\mathbb{N}}_0^d(M)}{w}_m{\mathbb{E}}_{a+m}[g(Z)],$$

$${\mathbb{E}}_a[e^{t_2{Z}_2+\cdots +t_d{Z}_d}]=C^{-1}\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right)\prod _{k=2}^d\frac{{\bar{a}}_k}{{\bar{a}}_k-t_k+{\bar{m}}_k};t_k<{\bar{a}}_k\text{for }k=2,\dots ,d,$$

$${\mathbb{E}}_a\left[\prod _{k=2}^d{Z}_k^{n_k}\right]=C^{-1}\sum _{m\in {\mathbb{N}}_0^d(M)}\left(\genfrac{}{}{0ex}{}{M}{m_1,\dots ,m_d}\right)\prod _{k=2}^d\frac{{\bar{a}}_k{n}_k!}{({\bar{a}}_k+{\bar{m}}_k{)}^{n_k+1}}.$$

$${\mathbb{E}}_a\left[\frac{1}{Y_1^r}\right]=\sum _{k=0}^{\infty }\left(\genfrac{}{}{0ex}{}{r}{k}\right)\sum _{j=0}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right)(-1{)}^{k-j}{d}^{r-j}\sum _{m\in {\mathbb{N}}_0^d(j)}(\genfrac{}{}{0ex}{}{j}{m_1,\dots ,m_d})\prod _{i=2}^d\frac{{\bar{a}}_i}{{\bar{a}}_i+{\bar{m}}_i}.$$

$${\mathbb{E}}_a\left[\frac{1}{Y_1^r}\right]=d^r\sum _{k=0}^{\infty }\left(\genfrac{}{}{0ex}{}{r}{k}\right){\mathbb{E}}_a\left[{\left(\frac{1}{d{Y}_1}-1\right)}^k\right].$$

$${\mathbb{E}}_a[f(Y)]=\sum _{k=0}^{\infty }\sum _{j=0}^k\sum _{m\in {\mathbb{N}}_0^d(j)}{w}_m^{r,j,k}{\mathbb{E}}_{a+m+(r-j)e_1}[f(Y)],$$

$$w_m^{r,j,k}=C^{-1}\left(\genfrac{}{}{0ex}{}{r}{k}\right)\left(\genfrac{}{}{0ex}{}{k}{j}\right)(-1{)}^{k-j}{d}^{r-j}(\genfrac{}{}{0ex}{}{j}{m_1,\dots ,m_d})\prod _{i=2}^d\frac{{\bar{a}}_i}{{\bar{a}}_i+{\bar{m}}_i}$$

$${\mathbb{E}}_a[f(Y)]=\frac{{\mathbb{E}}_{a+r{e}_1}[f(Y)Y_1^{-r}]}{{\mathbb{E}}_{a+r{e}_1}[Y_1^{-r}]}$$

Full text

Keywords:

Generalized Rank Dirichlet Distribution, Dirichlet Distribution, Poisson–Dirichlet Distribution, Exponential Distribution, Ordered Simplex, Ranked Weights.

MSC 2020 Classification:

Primary 60E05; Secondary 62G30

1 Introduction

For an integer $d\geq 2$ we study a parametric family of distributions defined on the ordered simplex

[TABLE]

whose density is proportional to

[TABLE]

for a parameter $a\in\mathbb{R}^{d}$. It was shown in [6] (and reproduced below in Proposition 1) that this density induces a probability measure on $\nabla^{d-1}$, when appropriately normalized, if and only if

[TABLE]

Notably, condition (2) allows for certain $a_{k}$’s to be negative as long as the tail sum $\bar{a}_{k}$ remains positive. In fact, even parameters satisfying $a_{k}<0$ for $k=1,\dots,d-1$ can be compatible with condition (2) (as long as $a_{d}$ is sufficiently positive).

In the special case $a_{1}=a_{2}=\dots=a_{d}>0$, if $X\sim\mathrm{Dirichlet}(a)$ then the ranked vector of decreasing order statistics $Y=(X_{(1)},\dots,X_{(d)})$ has density proportional to (1). In the case that the components of $a$ are not all the same this relationship is no longer true. However, since the functional form of (1) is the same as for the Dirichlet density – just defined on the ordered simplex rather than the standard simplex – we call the induced probability distribution the generalized ranked Dirichlet distribution with parameter $a$, or GRD($a$) for short.

The GRD distribution can be used to model the distribution of ranked weight vectors that sum to one even for a general $a$ parameter. Indeed, if $X=(X_{1},\dots,X_{d})$ is a random (unordered) vector of nonnegative weights that sum to one with density proportional to $\prod_{k=1}^{d}x_{(k)}^{a_{k}-1}$ then the decreasing order statistics $Y=(X_{(1)},\dots,X_{(d)})$ follow a GRD($a$) distribution.

To the best of the author’s knowledge the general form of the GRD($a$) distribution under the condition (2) first appeared as the invariant density of a certain stochastic process, called a rank Jacobi process in [6]. Previously, the special case with $\bar{a}_{1}=\sum_{k=1}^{d}a_{k}=0$ had appeared in [1, 8, 4, 3], where it arose as the invariant measure to a class of processes known as Atlas or first-order models. In particular, in [1], a connection to independent exponential random variables via the log gaps (see equation (3) below) was established. The analysis in this paper heavily exploits this relationship to exponential random variables in the case $\bar{a}_{1}=0$ to study GRD($a$) distributions for more general parameters $a$.

Arguably, the most well-studied distribution that models ranked weights is the Poisson–Dirichlet (PD) distribution introduced by Kingman in [7]. Indeed, it has found applications in a large number of fields including population genetics, number theory, physics, finance and statistics (see [9, 2] for detailed accounts of the PD distribution). However, it is defined on the infinite dimensional Kingman simplex $\{y\in\mathbb{R}^{\infty}:y_{1}\geq y_{2}\geq\dots\geq 0,\ \sum_{k=1}^{\infty}y_{k}=1\}$ and as such is an infinite-dimensional distribution. In the author’s PhD thesis [5], it was shown that, under appropriate assumptions on the parameter vector, the GRD distribution converges as $d\to\infty$ to a distribution on the Kingman simplex which is absolutely continuous with respect to a PD distribution with an explicitly given density. As such, the GRD family can be viewed as a finite dimensional relative of the PD distribution.

Remarkably, even in the most basic case $d=2$, the GRD distribution does not in general seem to be a standard probability distribution with a previously recorded name. When $d=2$ we can write $Y_{2}=1-Y_{1}$ and reduce to a one-dimensional random variable $Y_{1}$, which has density proportional to

[TABLE]

When $a_{1}>0$ this coincides with a truncated Beta distribution, but the case $a_{1}\leq 0$ does not seem to have an established name.

Nevertheless, this distribution has remarkable structural properties. In Section 2 we formally define the GRD distribution. Under the condition $\bar{a}_{1}=0$ the aforementioned relationship to independent exponential distributions is explored in Section 3, which we use to obtain negative moments of all orders for the largest weight $Y_{1}$. In Section 4 we then obtain a change of measure identity which establishes a relationship between GRD distributions with different parameters. In Section 5 we explore the case $\bar{a}_{1}=-M$ for some positive integer $M$. In this case the change of measure formula can be leveraged to obtain explicit expressions for the positive moments of the $Y_{k}$’s up to order $M$, which are derived in Section 5.1. In particular, when $M=1$, the moment formula is invertible with respect to the parameter vector $a$ allowing for explicit first moment matching. Additionally, it is shown in Section 5.3 that the log gaps

[TABLE]

can be represented as a mixture of exponential random variables when $\bar{a}_{1}=-M$. This leads us to explicit formulas for the moment generating function and moments of all orders for the log gaps. Using the log gaps as an intermediary, in Section 5.4, we derive an algorithm to simulate exactly from the GRD($a$) distribution in the case $\bar{a}_{1}=-M$. The general case when $\bar{a}_{1}$ is not assumed to be a negative integer is studied in Section 6. In this case we obtain a series representation for moments of the log gaps and leverage this to propose an approximate simulation algorithm to generate GRD($a$) random variates.

Notation.

The tail sum notation of $\bar{a}_{k}=a_{k}+a_{k+1}+\dots+a_{d}$, as in (2), is in force throughout the paper. We write $e_{1},\dots,e_{d}$ for the standard basis vectors in $\mathbb{R}^{d}$. We denote by $\mathbb{N}$ the natural numbers (starting from one) and $\mathbb{N}_{0}=\mathbb{N}\cup\{0\}$. For an integer $M>0$ we define $\mathbb{N}_{0}^{d}(M)=\{m\in\mathbb{N}^{d}_{0}:\bar{m}_{1}=M\}$. By convention, empty sums are taken to be zero, while empty products are taken to be one. Since $\nabla^{d-1}$ is a $(d-1)$-dimensional subset of $\mathbb{R}^{d}$, all integrals over $\nabla^{d-1}$ should be understood as the pushforward of Lebesgue measure on $\mathbb{R}^{d-1}$ under the map $(y_{1},\dots,y_{d-1})\mapsto(y_{1},\dots,y_{d-1},1-y_{1}-\dots-y_{d-1})$.

2 The GRD Distribution

Given $a\in\mathbb{R}^{d}$ we set $Q_{a}=\int_{\nabla^{d-1}}\prod_{k=1}^{d}y_{k}^{a_{k}-1}\,dy$. Then we have the following result already established in [6]. The proof is short and insightful so we reproduce it here.

Proposition 1** (Finite normalizing constant).**

$Q_{a}<\infty$* if and only if $\bar{a}_{k}>0$ for $k=2,\dots,d$.*

• Proof.

First note that the size or sign of $a_{1}$ does not effect integrability of $Q_{a}$ since $1/d\leq y_{1}\leq 1$. Hence we assume without loss of generality that $a_{1}=-\bar{a}_{2}$. Then we rewrite the integral as

[TABLE]

Next consider the change of variables $z_{k}=\log(y_{k-1})-\log(y_{k})$ for $k=2,\dots,d$. This transformation maps the ordered simplex onto $\mathbb{R}_{+}^{d-1}$ and its Jacobian is determined by $dz=\prod_{k=1}^{d}y_{k}^{-1}dy$. Thus we obtain

[TABLE]

This expression is finite if and only if $\bar{a}_{k}>0$ for every $k=2,\dots,d$ completing the proof. ∎

This leads us to the standing assumption mentioned in the introduction.

Assumption 2**.**

The parameter vector $a\in\mathbb{R}^{d}$ satisfies $\bar{a}_{k}>0$ for $k=2,\dots,d$.

We can now formally define the GRD distribution.

Definition 3** (Generalized Rank Dirichlet (GRD) Distribution).**

For a parameter $a\in\mathbb{R}^{d}$ satisfying Assumption 2 the probability measure

[TABLE]

is called a Generalized Rank Dirichlet (GRD) distribution with paremeter $a$. We will write $Y\sim\mathrm{GRD}(a)$ for a random variable $Y$ with law $\mathbb{P}_{a}$ and denote by $\mathbb{E}_{a}[\cdot]$ expectation under $\mathbb{P}_{a}$.

3 The case $\bar{a}_{1}=0$

An important special case of interest is when $\bar{a}_{1}=0$. In this case a similar calculation as in the proof of Proposition 1 shows that the log gaps $(Z_{2},\dots,Z_{d})$ given by (3) are distributed as independent exponentially distributed random variables whenever $Y\sim\mathrm{GRD}(a)$, and consequently, the weight ratios $Y_{k-1}/Y_{k}$ follow a Pareto distribution. Moreover, the normalizing constant $Q_{a}$ is explicitly computable in this case. To the best of the author’s knowledge the Pareto property was first observed in [3] and the relationship to independent exponential random variables was explored in [1]. We collect these results in the following proposition.

Proposition 4** (Section 4 in [1]).**

When $\bar{a}_{1}=0$ we have that $Q_{a}=\prod_{k=2}^{d}\bar{a}_{k}^{-1}$. Additionally the log gaps $(Z_{2},\dots,Z_{d})$ are independent and satisfy $Z_{k}\sim\mathrm{Exp}(\bar{a}_{k})$, while the ratios $Y_{k-1}/Y_{k}$ are independent and satisfy $Y_{k-1}/Y_{k}\sim\mathrm{Pareto}(1,\bar{a}_{k})$ for $k=2,\dots,d$.

These facts can be leveraged to compute certain expected ratios and negative moments of $Y_{1}$.

Theorem 5**.**

Let $a\in\mathbb{R}^{d}$ satisfying Assumption 2 be given and suppose that $\bar{a}_{1}=0$.

1. (i)

(Moments of ratios) Let $n\in\mathbb{N}_{0}^{d}$ and $M\in\mathbb{N}$ such that $M\geq\bar{n}_{1}$ be given. Then

[TABLE] 2. (ii)

(Negative moments of $Y_{1}$) For any $M\in\mathbb{N}$,

[TABLE]

• Proof.

First we assume that $\bar{n}_{1}=M$. In this case note that the expectation on the left hand side of (4) is given by $Q_{a+n-Me_{1}}/Q_{a}$. Since $(\overline{a+n-Me_{1}})_{1}=0$ we obtain

[TABLE]

by Proposition 4, which proves (i) in this case.

To prove (i) in the general case we use the multinomial formula to obtain

[TABLE]

In the last equality we used (6), which is applicable since $\bar{n}_{1}+\bar{m}_{1}=M$. Finally (5) follows by taking $n=0$ in (4). ∎

4 A change of measure formula

We now derive a change of measure identity, which holds for any GRD distribution. This identity is the workhorse for the computations to come.

Theorem 6** (Change of measure).**

Fix $a,b\in\mathbb{R}^{d}$ satisfying Assumption 2. Let $f:\nabla^{d-1}\to\mathbb{R}$ be a function that is integrable under $\mathbb{P}_{a}$. Then

[TABLE]

• Proof.

We see that

[TABLE]

where in the intermediate equality we multiplied and divided by $Q_{b}=\int_{\nabla^{d-1}}\prod_{k=1}^{d}y_{k}^{b_{k}-1}\,dy.$ ∎

As we saw in Section 3, the case when the sum of the parameters is zero is particularly tractable. Thus a canonical choice for the vector $b$ in the change of measure formula is $b=a-\bar{a}_{1}e_{1}$, in which case $\bar{b}_{1}=0$. Under this choice (7) becomes

[TABLE]

5 The case $\bar{a}_{1}=-M$

5.1 Moments of the $Y_{k}$’s

Remarkably, the identities for the negative moments of $Y_{1}$ when $\bar{a}_{1}=0$ can be used to derive positive moments, up to order $M$, for a GRD($a$) distribution when $\bar{a}_{1}=-M$. This is the content of the next theorem.

Theorem 7** (Moment formulas for $\bar{a}_{1}=-M$).**

Suppose that $a\in\mathbb{R}^{d}$ satisfies Assumption 2 and that $\bar{a}_{1}=-M$ for some $M\in\mathbb{N}$. Then for any $n\in\mathbb{N}^{d}_{0}$ with $\bar{n}_{1}\leq M$ we have that

[TABLE]

• Proof.

This follows directly by taking $f(Y)=\prod_{k=1}^{d}Y_{k}^{n_{k}}$ in (8) and invoking Theorem 5 to compute the right hand side of (8). ∎

When $M=1$ this formula takes a particularly simple form

[TABLE]

In particular this formula is invertible, which allows for explicit first moment matching, which can be used to calibrate the parameters to data.

Corollary 8** (First moment matching).**

Let $y\in\nabla^{d-1}$ satisfying $y_{1}>y_{2}>\dots>y_{d}$ be given. Define $a\in\mathbb{R}^{d}$ via

[TABLE]

Then $a$ satisfies Assumption 2, $\bar{a}_{1}=-1$ and $\mathbb{E}_{a}[Y_{k}]=y_{k}$ for $k=1,\dots,d$.

• Proof.

This is readily verified by applying (9) to this choice of $a$. ∎

5.2 An improved change of measure formula

In the case that $\bar{a}_{1}=-M$ for some $M\in\mathbb{N}$, the denominator of (8) is explicitly computable courtesy of Theorem 5. By writing $1=(Y_{1}+\dots+Y_{d})^{M}$ we can also expand the numerator to obtain that

[TABLE]

where the intermediate equality followed from Theorem 6 (with $a$ taken to be $a+m$ and $b$ taken to be $a+Me_{1}$ in the notation of the theorem), while the final equality followed from Theorem 5(i) since $\overline{(a+Me_{1})}_{1}=0$. This leads us to the following improved change of measure formula.

Theorem 9** (Change of measure v2).**

Let $a\in\mathbb{R}^{d}$ satisfying Assumption 2 be given and suppose that $\bar{a}_{1}=-M$ for some $M\in\mathbb{N}$. Then we have that

[TABLE]

for any $\mathbb{P}_{a}$-integrable function $f:\nabla^{d-1}\to\mathbb{R}$.

Since the $w_{m}$’s appearing in (11) are positive weights which sum to one, Theorem 9 establishes that $\mathbb{P}_{a}$ can be explicitly represented as a mixture of GRD distributions with parameters that sum to zero. This relationship can be leveraged to obtain certain moment formulas for the weights and log gaps, which are explored in the sections below. Additionally, marginal distributions for the weights under the GRD($a$) distribution can be studied with this change of measure identity as well, though we do not pursue this direction in detail here.

5.3 The log gaps as a mixture of exponential random variables

The change of measure formula of Theorem 9 is particularly insightful when we consider the log gaps $Z_{k}=\log Y_{k-1}-\log Y_{k}$ for $k=2,\dots,d$. Indeed, since $Z$ is a function of $Y$, we readily obtain the following corollary to Theorem 9.

Corollary 10** (Change of measure for log gaps).**

Let $a\in\mathbb{R}^{d}$ satisfying Assumption 2 be given and suppose that $\bar{a}_{1}=-M$ for some $M\in\mathbb{N}$. For any function $g:\mathbb{R}^{d-1}_{+}\to\mathbb{R}$ such that $g(Z)$ is $\mathbb{P}_{a}$-integrable we have

[TABLE]

where $w_{m}$ is defined in (11). In particular the the log gaps $(Z_{2},\dots,Z_{d})$ under $\mathbb{P}_{a}$ are a mixture of independent exponential random vectors.

• Proof.

The formula (12) is a direct consequence of Theorem 9, while the claim regarding the mixture of independent exponential distributions follows from Proposition 4 and the fact that $\bar{a}_{1}+\bar{m}_{1}=0$ for every $m\in\mathbb{N}^{d}_{0}(M)$. ∎

As an application of Corollary 10 we obtain the moment generating function and moments of the log gaps.

Corollary 11** (Log gap moments).**

Let $a\in\mathbb{R}^{d}$ satisfying Assumption 2 be given and suppose that $\bar{a}_{1}=-M$ for some $M\in\mathbb{N}$. Set $C=\mathbb{E}_{a+Me_{1}}[1/Y_{1}^{M}]$, which is explicitly given by (5) since $(\overline{a+Me_{1}})_{1}=0$. Then

1. (i)

the moment generating function of the log gaps $Z_{2},\dots,Z_{d}$ is given by

[TABLE] 2. (ii)

for any $n=(n_{2},\dots,n_{d})\in\mathbb{N}_{0}^{d-1}$ we have that

[TABLE]

• Proof.

This follows directly from Corollary 10 and known formulas for exponential random variables. ∎

5.4 Generation of random variates

We finish Section 5 by discussing a way to simulate a random vector $Y$ following a $\mathbb{P}_{a}$ distribution when $\bar{a}_{1}=-M$. This can be done by first simulating the log gap random vector $Z$ under $\mathbb{P}_{a}$ using the relationship in Corollary 10 and then inverting the maps $Y\mapsto(Z_{2},\dots,Z_{d})=(\log Y_{1}-\log Y_{2},\dots,\log Y_{d-1}-\log Y_{d})$. To carry this out we define a random variable $V$ on $\mathbb{N}^{d}_{0}(M)$ via $\mathbb{P}(V=m)=w_{m}$. The simulation steps are then as follows

This ensures that $Y\sim\mathbb{P}_{a}$. We note that the presentation of the algorithm above is simply pseudocode and the implementation can be made more efficient by vectorizing the operations.

6 The General Case

In the case that $\bar{a}_{1}\neq-M$ the change of measure formula can still be used to study the GRD distributions. Indeed, by applying Newton’s generalized binomial theorem we can obtain a series representation $\mathbb{E}_{a}[Y_{1}^{-r}]$ for arbitrary $r\in\mathbb{R}$ in the case $\bar{a}_{1}=0$.

Proposition 12** (Expected powers of $Y_{1}$).**

Let $a\in\mathbb{R}^{d}$ satisfying Assumption 2 be given and suppose that $\bar{a}_{1}=0$. Then for any $r\in\mathbb{R}$ we have

[TABLE]

• Proof.

We write $1/Y_{1}=d(1+\frac{1-dY_{1}}{dY_{1}})$. Note that since $1/d\leq Y_{1}\leq 1$ we have that $|\frac{1-dY_{1}}{dY_{1}}|<1$. Hence, applying Newton’s binomial theorem and taking expectation yields

[TABLE]

Now applying the standard binomial theorem to the term inside the expectation and using the identity derived in Theorem 5(ii) completes the proof. ∎

We now combine this with the change of measure formula to obtain the following theorem.

Theorem 13** (Change of measure series representation).**

Let $a\in\mathbb{R}^{d}$ satisfying Assumption 2 be given and suppose that $\bar{a}_{1}=-r$ for some $r\in\mathbb{R}$. Then for any $\mathbb{P}_{a}$-integrable function $f:\nabla^{d-1}\to\mathbb{R}$ we have that

[TABLE]

where

[TABLE]

and $C=\mathbb{E}_{a+re_{1}}[1/Y_{1}^{r}]$ is given explicitly by (13).

• Proof.

From the change of measure identity (8) we have that

[TABLE]

The denominator has the series representation given by Proposition 12. To handle the numerator we use Newton’s binomial theorem to expand out $Y_{1}^{-r}=d(1+\frac{1-dY_{1}}{dY_{1}})$ as before, multiply both sides by $f(Y)$ and take expectation to obtain

[TABLE]

where we used the standard binomial theorem in the final equality. Proceeding as in (10) we obtain

[TABLE]

Plugging this into (15) completes the proof. ∎

The upshot of this theorem is that we can represent an arbitrary GRD($a$) distribution as a countable mixture of GRD distributions where the parameter vectors sum to zero. Applying this to the log gap process $Z$ as in Section 5.3 shows, in turn, that the log gaps under an arbitrary GRD($a$) distribution are a countable mixture of independent exponential random variables. This leads to series representation formulas for the log generating function and moments of the log gaps.

Corollary 14** (Log gap moments series representation).**

Let $a\in\mathbb{R}^{d}$ satisfying Assumption 2 be given. Then

1. (i)

the moment generating function of the log gaps $Z_{2},\dots,Z_{d}$ is given by

[TABLE] 2. (ii)

for any $n=(n_{2},\dots,n_{d})\in\mathbb{N}_{0}^{d-1}$ we have that

[TABLE]

where $w_{m}^{-{\bar{a}_{1}},j,k}$ is defined in the statement of Theorem 13.

Moreover, the representation of $Z$ as a countable mixture of independent exponential random variables suggests an approximate algorithm for generating random GRD($a$) variates for arbitrary parameter $a$ by truncating the series appearing in (14). If we keep the first $K+1\in\mathbb{N}$ terms in the series then by rearranging the terms in the sum we obtain from (14) that

[TABLE]

where

[TABLE]

Consequently, if we define the random variable $V^{K}$ on the discrete set $\{m\in\mathbb{N}_{0}^{d}:\bar{m}_{1}\leq K\}$ via

[TABLE]

then we obtain an algorithm to approximately sample from the GRD($a$) distribution for arbitrary parameter $a$.

7 Conclusion

We introduced the family GRD($a$) of distributions on the ordered simplex $\nabla^{d-1}$. We established change of measure formulas that relate GRD($a$) distributions with different parameters to each other. In the case that $\bar{a}_{1}=-M$ for some $M\in\mathbb{N}$ we exploited the change of measure identity to show that such a distribution is a (finite) mixture of GRD distributions with parameters that sum to zero. This, together with the fact that the log gaps $Z$ are independent exponential random variables when the parameters sum to zero, was used to establish moment formulas, up to order $M$, for the weights as well as for moments of all orders for the log gaps. This led to an algorithm which allows one to exactly sample the weights $Y$. In the case $M=1$, the first moment formula is invertible allowing for explicit moment matching which can be used for calibration to data. In the general case when $\bar{a}_{1}\in\mathbb{R}$, we were able to recover many of the same properties, but under series representations rather than finite sums. This led us to an algorithm for approximately sampling the weights $Y$ in this case.

Acknowledgements.

I am grateful to Martin Larsson for helpful discussions.