2. Background

This section provides the necessary insights on the structure and main properties of the coupled Wright–Fisher family of diffusions and fixes some notation. It also provides a brief overview of exact rejection algorithms for diffusions, which constitute the basis of our work.

2.1. Coupled Wright–Fisher diffusions

The family of Wright–Fisher models, and more specifically their diffusion approximations, have been widely used in population genetics; see, for instance, [Reference Kimura36] and [Reference Griffiths24]. In its simplest form, the Wright–Fisher model describes the evolution of the frequency of two allele types in a single locus that have the same fitness, and whose configuration at each new generation of individuals is chosen uniformly and with replacement from that of the current generation in a haploid population of constant size. Extensions of the model include consideration of more than two allele types that might be located at different loci, and can incorporate other evolutionary forces such as mutation, selection, and recombination. A comprehensive overview of the family of Wright–Fisher models can be found, for example, in [Reference Crow and Kimura12] or [Reference Ewens18].

With the proliferation of genome-wide association studies, questions arise about how genetic variants associated to numerous diseases co-evolve or interact over time. Moreover, the increasing availability of allele frequency time series data is fostering the study of evolutionary forces such as mutation or selection; see [Reference Steinrücken, Bhaskar and Song44], [Reference Skwark43], [Reference Tataru, Simonsen, Bataillon and Hobolth45], and [Reference Nené, Mustonen and Illingworth39]. Within this framework, the coupled Wright–Fisher model [Reference Aurell, Ekeberg and Koski1] tracks the evolution of frequencies of allele types located at different loci, and, as well as locus-wise mutation and selection, describes possible selective pairwise allele interactions between loci in a haploid population. Its diffusion approximation can be derived as the weak limit of a sequence of discrete coupled Wright–Fisher models characterized by the assumption that the evolution of the population at one locus is conditionally independent of the other loci given the state of the previous generation. The coupled Wright–Fisher diffusion can be expressed as a system of stochastic differential equations of the form

(2.1) \begin{equation}dX_t=[\alpha(X_t)+G(X_t)]dt+ D^{\frac{1}{2}}(X_t)dB_t, \ X_0=x_0,\ t\in[0,T],\end{equation}

where $X_t$ is a vector of frequencies of allele types, $\alpha$ governs their mutations, and G contains the single and pairwise selective locus interactions.

Let L denote the total number of loci and $d_i\geq 2$ the number of different allele types in each of them. For $n\,:\!=\,\sum_{i=1}^L [d_i-1]$, let us index each element of an n-dimensional vector x by $i\in\{1,\ldots , L\}$, referring to a specific locus, and $j\in\{1,\ldots, d_i-1\}$, referring to an allele type, so that $x=\{x^i\}_{i=1}^L$ where each $x^i=\{x^{ij}\}_{j=1}^{d_i-1}$. If $x\in\mathbb{R}^n$ refers to the vector of allele type frequencies $X_t$, the elements of the drift $\alpha(x)\in\mathbb{R}^n$ take the form

(2.2) \begin{equation}\alpha^{ij}(x^{ij})=\frac{1}{2}\big(\theta^{i}_{j}-|\theta|x^{ij}\big),\end{equation}

where $|\theta|=\sum_{k=1}^{d_i}\theta^{i}_{k}$ and $\theta^{i}_{k}>0$ denote the parent-independent mutation rates to allele type $k\in\{1,\ldots, d_i\}$ at locus i, so that mutations occur at each locus separately. Wright–Fisher diffusions with drift $\alpha(x)$ correspond to the reversible neutral mutations allele model. The coupled Wright–Fisher model also admits parent-dependent mutations, but we will not consider these here.

The coupling term $G(x)\in\mathbb{R}^n$ has general form (Svirezhev–Shahshahani gradient)

\begin{equation*}G(x)=D(x)\nabla_{x} (\overline{V}\circ f)(x),\end{equation*}

where the square of the diffusion matrix $D(x)=\text{diag}(D^i(x^i))\in\mathbb{R}^{n\times n}$ is an L-block diagonal matrix with entries

\begin{equation*}D^i_{jk}=\begin{cases} x^{ij}(1-x^{ij}), \qquad j=k, \\[3pt] -x^{ij}x^{ik}, \qquad j\neq k,\end{cases} j,k \in\{1,\ldots, d_i-1\}, i\in\{1,\ldots , L\};\end{equation*}

$\nabla_{x}$ is the gradient operator with respect to each component of x; f transforms x into the augmented $(n+L)$-dimensional vector $\overline{x}$ that reflects the dependency between allele frequencies at locus i, i.e.,

\begin{equation*}\overline{x}^{ik}\,:\!=\,f^{ik}(x)=\bigg\{\begin{array}{cl} x^{ij},& j=k\in\{1,\ldots, d_i-1\},\\[4pt] 1-\sum_{j=1}^{d_i-1}x^{ij},& k=d_i;\\\end{array}\end{equation*}

and $\overline{V}(\overline{x})\in \mathbb{R}$,

\begin{equation*}\overline{V}(\overline{x})=(\overline{x})^Ts+\frac{1}{2}(\overline{x})^TH\overline{x},\end{equation*}

where $s\in\mathbb{R}^{n+L}$ is the vector of within-locus selection parameters and $H\in\mathbb{R}^{(n+L)\times(n+L)}$ is a symmetric between-loci pairwise interactions matrix. The matrix H is in fact built by L blocks of zeros of size $d_i\times d_i$ in the main diagonal (denoted by $0^{ii}$), and off-diagonal blocks of the form $H^{il}=(H^{li})^T\in\mathbb{R}^{d_i\times d_l}$, $i\neq l$, $i,l\in\{1,\ldots, L\}$,

\begin{equation*}H=\begin{pmatrix} 0^{11}& \quad H^{12}& \quad \ldots& \quad \ldots& \quad \ldots& \quad \ldots& \quad H^{1L}\\[1pt] \vdots& \quad \ddots& \quad \vdots& \quad \vdots& \quad \vdots & \quad \vdots & \quad \vdots\\[4pt] H^{i1}& \quad \ldots& \quad H^{i(i-1)}& \quad 0^{ii}& \quad H^{i(i+1)}& \quad \ldots& \quad H^{iL} \\[4pt] \vdots& \quad \ddots& \quad \vdots& \quad \vdots& \quad \ddots & \quad \vdots & \quad \vdots\\[1pt] H^{L1}& \quad \ldots& \quad \ldots& \quad \ldots& \quad \ldots& \quad \ldots& \quad 0^{LL} \\ \end{pmatrix},\end{equation*}

so that interactions of each locus with itself are not permitted. Note that if one removes the coupling term G, (2.1) simply becomes the usual multivariate neutral mutations Wright–Fisher diffusion, where $X_t$ describes the evolution of allele frequencies that evolve at each locus without interacting. The explicit form of $V(x)\,:\!=\,\nabla_{x} (\overline{V}\circ f)(x)$ in terms of s and H is specified in the following proposition, whose proof can be found in Section 7.

Proposition 2.1. Let $V(x)\in\mathbb{R}^n$ be the n-dimensional vector such that $V(x)\,:\!=\,\nabla_{x} (\overline{V}\circ f)(x)$. Then $\forall i \in\{1,\ldots, L\}$, $j\in\{1,\ldots, d_i-1\}$,

\begin{equation*}V^{ij}(x)=K_s^{ij}+\sum_{\underset{l\neq i}{l=1}}^{L}\bigg(K_l^{ij}+ \sum_{k=1}^{d_l-1} K_{lk}^{ij} x^{lk}\bigg),\end{equation*}

with

\begin{equation*}K_s^{ij}\,:\!=\,s^{ij}- s^{id_i}, \quad K_l^{ij}\,:\!=\,h^{il}_{jd_l}-h^{il}_{d_id_l}, \text{ and }K_{lk}^{ij}\,:\!=\,h^{il}_{jk}-h^{il}_{jd_l}-h^{il}_{d_ik}+h^{il}_{d_id_l},\end{equation*}

where $h^{il}_{jk}$ denotes the entry in row j and column k of the block $H^{il}$ of H.

Following Kimura’s formulation [Reference Kimura35], the explicit stationary density of (2.1) can also be obtained by solving the corresponding Fokker–Planck equation, see [Reference Aurell, Ekeberg and Koski1], whose solution takes the form

(2.3) \begin{equation} P(x)=\dfrac{1}{Z}\pi(x)e^{2(\overline{V}\circ f)(x)},\end{equation}

where

\begin{equation*}\pi(x)=\prod_{i=1}^L \pi^i(x^i)=\prod_{i=1}^L \left( (1-\sum_{j=1}^{d_i-1}x^{ij})^{\theta^i_{d_i}-1}\prod_{j=1}^{d_i-1}(x^{ij})^{\theta^i_j-1}\right)\end{equation*}

and Z is the normalizing constant

\begin{equation*}Z=\int_{x\in\mathcal{X}} \pi(x)e^{2(\overline{V}\circ f)(x)}dx.\end{equation*}

Here $\mathcal{X}=\{x\in\mathbb{R}^n\mid x^{ij}\geq 0, \sum_{j=1}^{d_i-1}x^{ij}\leq 1\}$.

The representation (2.3) resembles the stationary density of the haploid version of the model studied by Fearnhead (see [Reference Fearnhead20, Theorem 2]), which in the two-locus case agrees with the coupled Wright–Fisher diffusion.

2.2. Overview of exact rejection algorithms for simulation of diffusions

This subsection provides an overview of the exact rejection algorithm presented in [Reference Beskos and Roberts5], and presents the same sampling scheme followed for simulating from the coupled Wright–Fisher diffusion, as detailed further in Algorithm 4. Let

(2.4) \begin{equation} dX_t= \mu(X_t) dt+dB_t,\ X_0=x_0,\ t\in [0,T],\end{equation}

where $\mu({\cdot})$ is such that (2.4) admits a weakly unique solution $(X_t)_{t\in[0,T]}$. Let $\mathbb{Q}_{x_0}$ be the law of the process $(X_t)_{t\in[0,T]}$, and let $\mathbb{P}_{x_0}$ denote the law of a Brownian motion $(B_t)_{t\in[0,T]}$ starting at $B_0=x_0$. By means of a Girsanov transformation of measures one can write the Radon–Nikodym derivative of $\mathbb{Q}_{x_0}$ with respect to $\mathbb{P}_{x_0}$ as

(2.5) \begin{equation} \dfrac{d\mathbb{Q}_{x_0}}{d\mathbb{P}_{x_0}}=\exp\left\{\int_0^T \mu(B_t) dX_t-\dfrac{1}{2}\int_0^T \mu^2(B_t)dt\right\}.\end{equation}

Assuming that $\mu({\cdot})$ is differentiable everywhere and using Itô’s lemma, (2.5) can be rewritten as

\begin{equation*} \dfrac{d\mathbb{Q}_{x_0}}{d\mathbb{P}_{x_0}}=\exp\left\{\tilde{A}(B_T)-\tilde{A}(x_0)\right\}\exp\left\{-\dfrac{1}{2}\int_0^T (\mu^2(B_t)+\mu'(B_t)) dt\right\},\end{equation*}

where $\tilde{A}(x)\,:\!=\,\int_0^x\mu(u)du$. Imposing the further conditions that $\tilde{A}(x)$ be bounded above by a constant $K^A$, and that $(\mu^2+\mu')/2$ be bounded between constants $K^-$ and $K^+$, we find that

(2.6) \begin{equation} \dfrac{d\mathbb{Q}_{x_0}}{d\mathbb{P}_{x_0}}\propto\exp\left\{\tilde{A}(B_T)-K^A\right\} \exp\left\{-\int_0^T (\tilde{\phi}(B_t) -K^- )dt\right\},\end{equation}

with $\tilde{A}(x)\leq K^A$ and $K^-\leq \tilde{\phi}(x)\,:\!=\,\dfrac{1}{2}[\mu^2(x)+\mu'(x)]\leq K^+$, is a suitable acceptance–rejection probability.

Note that an exact evaluation of the integral in (2.6) is not possible without any approximation error, because this would require the storage of infinitely many points of the sample candidate path $B=(B_t)_{t\in[0,T]}$. However, the key point of the exact algorithms proposed in [Reference Beskos and Roberts5] and references therein is that an event occurring with probability (2.6) can be evaluated by sampling B only at a finite number of time points. This follows because the last term in (2.6) can be interpreted as the probability of the event $\omega_{\tilde{\phi}}$ that no points from a homogeneous spatial Poisson process $\Phi=\{(t_j,\psi_j): j=1,\ldots,J\}$ with unit intensity on $[0,T]\times[0,K^+-K^-]$ lie below the graph of $t\mapsto\tilde{\phi}(B_t)$. A formal statement and proof of this observation can be found in Theorem 1 of [Reference Beskos, Papaspiliopoulos and Roberts3].

Therefore, the exact sampling procedure (detailed in Algorithm 1) starts by drawing a sample from $\Phi$ that will determine the (random) time points at which the candidate B will be drawn, and then provides a skeleton of $(X_t)_{t\in[0,T]}$ at such time points upon acceptance of the sampled candidate (that is, if all the evaluated $\tilde{\phi}(B_t)$ lie above the sampled Poisson points). Note that the last point on the candidate path, $B_T$, serves to evaluate an event that occurs with probability $\exp\{\tilde{A}(B_T)-K^A\}$ and that is independent of $\omega_{\tilde{\phi}}$. In previously presented versions of the algorithm, $B_T$ is considered to follow a specific probability distribution, thus slightly modifying the candidate to be a certain biased Brownian motion, which is shown to improve the algorithm’s efficiency. Such an option is not needed for our purposes and is therefore not fully described here. In brief, this modification permits the boundedness condition on $\tilde{A}(x)$ to be relaxed, but the equivalent function in our proposed exact algorithm is already bounded.

3. Simulation of the candidate processes

This section is devoted to describing existing simulation strategies for the candidate processes in the exact rejection algorithm that will be presented later in Section 4. Suitable candidate processes in our setting will be L independent $(d_i-1)$-dimensional neutral Wright–Fisher diffusions $(X_t)_{t\in [0,T]}$, each one a weakly unique solution of

(3.1) \begin{equation}dX_t= \alpha(X_t) dt+D^{\frac{1}{2}}dB_t,\quad X_0=x_0,\ t\in [0,T],\end{equation}

with $\alpha(X_t)$ a $(d_i-1)$-dimensional vector with $\alpha^{ij}(x)=\frac{1}{2}(\theta^i_j-|\theta| x^{ij})$.

3.1. Transition density function expansions

Exact simulation of each neutral Wright–Fisher diffusion is possible by exploiting an available eigenfunction expansion of its transition density function that allows a probabilistic representation; see, for example, [Reference Griffiths, Spanò, Bingham and Goldie29].

For a fixed locus $i\in\{1,\ldots, L\}$, let $x=(x_1,\ldots, x_{d-1})$ be a vector of initial frequencies and $\theta+l$ a d-dimensional vector with entries $\theta_j+l_j, \ j\in\{1,\ldots, d\}$. Then the probabilistic representation of the transition density function of $(X_t)_{t\in [0,T]}$ in (3.1) is given by

(3.2) \begin{equation}g(x,y;t)=\sum_{m=0}^{\infty} q_m^{\theta}(t) \sum_{\underset{|l|=m}{l}} \mathcal{M}_{m,x}(l)\mathcal{D}_{\theta+l}(y) ,\end{equation}

where the $q_m^{\theta}(t)$ are transition functions of a pure death process $A_{\infty}^{\theta}(t)$ with an entrance boundary at $\infty$, $\mathcal{M}_{m,x}({\cdot})$ is the probability mass function (PMF) of a multinomial random variable, and $\mathcal{D}_{\theta+l}({\cdot})$ is the probability density function of a Dirichlet random variable; that is,

\begin{equation*} \mathcal{M}_{m,x}(l)=\frac{m!}{\prod_{j=1}^{d}l_j!}(1-\sum_{j=1}^{d-1}x_j)^{l_{d}} \prod_{j=1}^{d-1}x_j^{l_j} \end{equation*}

and

\begin{equation*} \mathcal{D}_{\theta+l}(y)=\frac{\Gamma( |\theta+l|)}{\prod_{j=1}^{d}\Gamma(\theta_j+l_j)} (1-\sum_{j=1}^{d-1}y_j)^{\theta_{d}+l_{d}-1}\prod_{j=1}^{d-1}y_j^{\theta_j+l_j-1}, \end{equation*}

with $l_d=m-\sum_{j=1}^{d-1} l_j$. A more detailed description of the process $A_{\infty}^{\theta}(t)$ and an exact sampling technique are provided later on.

In the case of one-dimensional ($d=2$) Wright–Fisher diffusions, the multivariate components of the mixture in (3.2) reduce to their one-dimensional counterparts, i.e., a binomial and a beta random variable, respectively [Reference Griffiths, Spanò, Bingham and Goldie29].

A sampling strategy for $g(x,\cdot;t)$ is summarized in Algorithm 2; see [Reference Griffiths and Li28] or [Reference Jenkins and Spanò32] for an analogous version in the one-dimensional case (the latter also includes a modification for the infinite-dimensional case, that is, for Fleming–Viot diffusions). Once expressed in probabilistic terms and given the simplicity of Algorithm 2, recovering samples from $g(x,\cdot;t)$ seems straightforward. However, sampling exactly from $q_m^{\theta}(t)$ poses some difficulty because it is only known in infinite series form. Previous approaches for simulating from approximated versions of $q_m^{\theta}(t)$ can be found in [Reference Griffiths and Li28], [Reference Griffiths25], [Reference Griffiths26], and [Reference Jewett and Rosenberg33]. In the next section, we review the exact simulation procedure presented in [Reference Jenkins and Spanò32], which is the one used here.

3.2. Exact simulation of the ancestral process $A_{\infty}^{\theta}$

We describe here the sampling procedure for recovering exact samples of $q_m^{\theta}(t)$, the transition functions of the aforementioned death process $A_{\infty}^{\theta}$. In more detail, let $\{A_n^{\theta}(t): t\geq 0\}$ be a pure death process on $\mathbb{N}$ such that $A_n^{\theta}(0)=n$ almost surely and with its only possible transition $m\to m-1$ occurring at rate $m(m+|\theta|-1)/2$ for each $m=1,\ldots,n$; that is, it represents the number of non-mutant lineages that coalesce backwards in time in the coalescent process with mutation. Then, let $q_m^{\theta}(t)=\lim_{n \to \infty}\mathbb{P}(A_n^{\theta}(t)=m)$.

For a more thorough interpretation of the transition density $g(x,\cdot, t)$ and its one-dimensional counterpart, it is worth noting that the expansion in (3.2) is derived via a duality principle for Markov processes [Reference Ethier and Kurtz16], that is, from the moment dual process of the Wright–Fisher diffusion, which is also a pure death process representing lineages backwards in time; see for example [Reference Etheridge and Griffiths15] or [Reference Griffiths, Spanò, Bingham and Goldie29] for a complete derivation and details. The corresponding dual (coalescent) process for coupled Wright–Fisher diffusions is derived in [Reference Favero, Hult and Koski19].

An expression for $q_m^{\theta}(t)$ starting from the entrance boundary at infinity is as follows (see [Reference Griffiths24]):

(3.3) \begin{equation} q_m^{\theta}(t)=\sum^\infty_{i=0}(-1)^{i}b_{m+i}^{(t,\theta)}(m),\end{equation}

where

(3.4) \begin{equation}b_{m+i}^{(t,\theta)}(m)=\frac{(|\theta|+2(m+i)-1)}{m!i!} \frac{\Gamma(|\theta|+2m+i-1)}{\Gamma(|\theta|+m)}e^{(m+i)(m+i+|\theta|-1)t/2}.\end{equation}

As shown in [Reference Jenkins and Spanò32], samples from $q_m^{\theta}(t)$ can be recovered exactly by means of a variant of the alternating series method, described in [Reference Devroye13, Chapter 4]. In brief, the alternating series method would require the sequence of coefficients $b_{m+i}^{(t,\theta)}(m)$ to be decreasing in i for each m, a condition that is not always met here. Nonetheless, one can exploit the fact that there exists a finite $C_m^{(t,\theta)}$ such that for all m and $i\geq C_m^{(t,\theta)}$ the sequence of coefficients $b^{(t,\theta)}_{m+i}(m)$ decreases monotonically as i tends to $\infty$. More explicitly, there exists

(3.5) \begin{equation}C_m^{(t,\theta)}=\inf\left\{i\geq 0 :\ (b^{(t,\theta)}_{m+i+1}(m)/b^{(t,\theta)}_{m+i}(m))<1 \right\}<\infty.\end{equation}

Then, once $C_m^{(t,\theta)}$ is available, the remaining of the sequence of coefficients is ensured to be decreasing and the alternating series method can be applied. The following proposition summarizes the main properties of the bound $C_m^{(t,\theta)}$.

Proposition 3.1. ([Reference Jenkins and Spanò32, Proposition 1]) Let $b_{m+i}^{(t,\theta)}(m)$ be the coefficients defined in (3.4) and let $C_m^{(t,\theta)}$ be as in (3.5). Then

1. (i) $C_m^{(t,\theta)}<\infty$ for all m;

2. (ii) $b_{m+i}^{(t,\theta)}(m)\downarrow 0$ as $i\to\infty$ for all $i\geq C_m^{(t,\theta)}$;

3. (iii) $C_m^{(t,\theta)}=0$ for all $m>D_{\epsilon}^{(t,\theta)}$, where

   (3.6) \begin{equation} D_{\epsilon}^{(t,\theta)}= \inf\left\{u\geq \left(\frac{1}{t}-\frac{|\theta|+1}{2}\right)\vee 0 : (|\theta|+2u-1)e^{u(u+|\theta|-1)t/2}<1-\epsilon \right\}, \end{equation}
   for $\epsilon\in[0,1)$.

Part (iii) of Proposition 3.1 will be of interest later in proposing the exact sampling algorithm for $(d-1)$-dimensional Wright–Fisher bridges (Section 6), where an explicit bound on m is needed.

Following Proposition 3.1, one can then recover exact samples from $q_m^{\theta}(t)$ because the terms $b^{(t,\theta)}_{m+i}(m)$ become monotonically smaller with increasing i, and for each m there exists $k_m$, an element of $\tilde{k}\in\mathbb{R}^{M+1}$ (i.e., $\tilde{k}=\{k_m\}_{m=0}^M$), such that

\begin{equation*}S^-_{\tilde{k}}(M)\,:\!=\,\sum_{m=0}^M\sum^{2k_m+1}_{i=0}(-1)^{i}b_{m+i}^{(t,\theta)}(m)\leq \sum_{m=0}^M q_m^{\theta}(t)\leq \sum_{m=0}^M\sum^{2k_m}_{i=0}(-1)^{i}b_{m+i}^{(t,\theta)}(m)=:S^+_{\tilde{k}}(M).\end{equation*}

Because

\begin{equation*}\limsup_{\tilde{k} \to (\infty, \ldots, \infty)}S^-_{\tilde{k}}(M)= \mathbb{P}(A_{\infty}^{\theta}(t)\leq M)\quad \text{and} \quad \liminf_{\tilde{k} \to (\infty, \ldots, \infty)}S^+_{\tilde{k}}(M)= \mathbb{P}(A_{\infty}^{\theta}(t)\leq M),\end{equation*}

and because both $S^-_{\tilde{k}}(M)$ and $S^+_{\tilde{k}}(M)$ can be computed from finitely many terms, given $U\sim \text{Uniform}(0,1)$ we can find $\tilde{k}^0\in\mathbb{R}^{M+1}$ with elements $k^0_m$ such that

\begin{equation*}k^0_m=\inf\left\{ k_m\in\mathbb{N}:\ S^-_{\tilde{k}}(M) > U \text{ or }S^+_{\tilde{k}}(M) < U\right\}\end{equation*}

for each $m\in\{0,\ldots, M\}$.

Now, if $\tilde{k}^0$ is such that $S^-_{\tilde{k}^0}(M) > U$, standard inverse sampling provides

(3.7) \begin{equation}\inf\left\{ M\in\mathbb{N}: \sum_{m=0}^M q_m^{\theta}(t) \geq S^-_{\tilde{k}^0}(M)> U\right\},\end{equation}

with M exactly distributed following $q_m^{\theta}(t)$. The sampling strategy will consist in exploring the summands in $S^-_{\tilde{k}}(M)$ and $S^+_{\tilde{k}}(M)$ through their respective indexes m and $k_m$, until, for a given realization of U, the condition (3.7) is satisfied.

A complete simulation procedure is presented in Algorithm 3, where several improvements mentioned in [Reference Jenkins and Spanò32] have been incorporated. Most notably, the variable M is initialized at the nearest integer around the mean of a certain distribution (not necessarily at 0) that serves as an estimate of the mode $\hat{q}_{\text{mod}}$ of $q_m^{\theta}(t)$, a modification that decreases computation times substantially. Such initialization originates from an asymptotic approximation of the transition functions $q_m^{\theta}(t)$, see [Reference Griffiths25], which states that as $t\to 0$, $A_{\infty}^{\theta}(t)$ converges to a normal distribution, that is,

(3.8) \begin{equation}\frac{A_{\infty}^{\theta}(t)-\mu^{(t, \theta)}}{\sigma^{(t, \theta)}}\stackrel{\mathcal{D}}{\longrightarrow} \mathcal{N}(0,1)\text{ as } t\to 0,\end{equation}

where $\mu^{(t, \theta)}=2\eta/t$ and

\begin{equation*}(\sigma^{(t, \theta)})^2=\begin{cases} \frac{2\eta}{t\beta^2}(\eta+\beta)^2\left(1+\frac{\eta}{\eta+\beta}-2\eta\right), \quad \beta\neq 0,\\[5pt] \frac{2}{3t}, \quad \beta= 0,\end{cases}\end{equation*}

with $\eta=\beta/e^{\beta}-1$ for $\beta\neq 0$ or $\eta=1$ otherwise, $\beta=\frac{1}{2}(|\theta|-1)t$, and where $\stackrel{\mathcal{D}}{\longrightarrow}$ denotes convergence in distribution; see [Reference Jenkins and Spanò32, Theorem 1].

Note that this initialization is possible because the exploration of the summands in $S^-_{\tilde{k}}(M)$ and $S^+_{\tilde{k}}(M)$ does not need to occur in any specific order. It is also worth mentioning that, when M is initialized at 0, the vector $\tilde{k}$ is updated increasingly, i.e. a new element of the vector is added at every new iteration where M is increased by one unit. In Algorithm 3, however, M is updated telescopically; that is, at each new iteration j, M moves farther from $\hat{q}_{\text{mod}}$ by one unit alternatingly above or below. This in turn entails updating the corresponding $M+1$ elements of $\tilde{k}$ accordingly; i.e., the number of elements might either increase or decrease at each iteration. In addition, because M is updated telescopically, when an M such that $S^-_{\tilde{k}}(M)>U$ is found, one needs to ensure whether this is the infimum M for which this condition is satisfied. For precise results on the complexity of Algorithm 3 and simulation performance we refer the reader to [Reference Jenkins and Spanò32].

At this stage, Algorithm 3 can be used in Step 1 of Algorithm 2, and an exact sampling procedure for $(d-1)$-dimensional neutral Wright–Fisher diffusions is completed.

4. Exact rejection algorithm for simulating coupled Wright–Fisher diffusions

Let $X_t$ be the n-dimensional vector of allele frequencies satisfying (2.1). Following the same scheme as Algorithm 1 in Section 2, Algorithm 4 simulates exact skeletons of paths of coupled Wright–Fisher diffusions with L loci and $d_i$ allele types each, $i\in\{1,\ldots,L\}$. Candidate processes in this case are L independent $(d_i-1)$-dimensional neutral Wright–Fisher diffusions, each one a weakly unique solution of (3.1) and sampled exactly following Algorithm 2.

The exact rejection algorithm proposed in this paper relies on the existence and characterization of the following acceptance–rejection probability, which is detailed in Theorem 4.1 and whose proof is deferred to Section 7.

Theorem 4.1. Let $\mathbb{CWF}_{\alpha,G,x_0}$ be the law of X, solution of (2.1), and let $\mathbb{WF}L_{\alpha,x_0}$ be the joint law of L independent $(d_i-1)$-dimensional neutral Wright–Fisher diffusions weakly unique solutions of (3.1), $i\in\{1,\ldots, L\}$. Then the Radon–Nikodym derivative of $\mathbb{CWF}_{\alpha,G,x_0}$ with respect to $\mathbb{WF}L_{\alpha,x_0}$ is of the form

\begin{equation*} \dfrac{d\mathbb{CWF}_{\alpha,G,x_0}}{d\mathbb{WF}L_{\alpha,x_0}} =\exp\left\{A(X_0,X_T)\right\}\exp\left\{-\int_0^T \phi(X_t) dt\right\},\end{equation*}

and there exist constants $A^-, A^+, C^-$, and $C^+$ such that $A(X_0,X_T)$ is bounded on $[0,1]^n\times[0,1]^n$ by $A^-\leq A(X_0,X_T)\leq A^+$ and $\phi(X_t)$ is bounded on $[0,1]^n$ by $C^-\leq\phi(X_t)\leq C^+$, with

\begin{multline*} A(X_0,X_T)\,:\!=\,\displaystyle\int_0^T V(X_t)\cdot dX_t= \sum_{i=1}^{L}\sum_{j=1}^{d_i-1} \bigg(K^{ij}_s (X^{ij}_T-X^{ij}_0) +\sum_{l=1}^{L}K^{ij}_l (X^{ij}_T-X^{ij}_0)\\ +\displaystyle\sum_{l=i+1}^{L}\sum_{k=1}^{d_l-1} K^{ij}_{lk}\ (X^{ij}_TX^{lk}_T-X^{ij}_0X^{lk}_0)\bigg)\end{multline*}

and

\begin{equation*}\phi(X_t)\,:\!=\,\dfrac{1}{2}\left[ (V(X_t))^T D(X_t) V(X_t) +2 (V(X_t))^T \alpha(X_t)\right].\end{equation*}

Using Algorithm 2 in Step 3 (or the corresponding modification for one-dimensional diffusions, if $d_i=2$ for some i), Algorithm 4 returns an exactly simulated skeleton of $(X_t)_{t\in[0,T]}$ solution of (2.1).

The algorithm’s computational complexity can also be established; this is made precise in Proposition 4.1, whose proof can be found in Section 7.

Proposition 4.1. Let L be the number of loci; let M(t) denote the total number of coefficients that must be computed in the implementation of Algorithm 3, where $t\in(0,T)$ is the time distance between two sampled skeleton points; and let N(T) denote the number of Poisson points required until the first skeleton in Algorithm 4 is accepted. Then $\text{E}[LM(t)]<\infty$ and $\text{E}[N(T)]<\infty$, and more specifically, there exists $\kappa>0$ such that

\begin{align*} \text{E}[LM(t)] &=o(t^{-(1+\kappa)}) \quad \text{ as } t\to 0, \text{ and } \\ \text{E}[N(T)] &\leq T(C^+-C^-)e^{T(C^+-C^-)+A^+-A^-}. \end{align*}

In summary, the complexity of Algorithm 4 increases either as $t\to 0$, when the average number of coefficients to be computed in Algorithm 3 increases as $1/t$, or with increasing T, when the average number of Poisson points needed until acceptance increases exponentially.

The latter is easily solvable, simply by considering shorter intervals $[t_{k-1}, t_k]$ such that

\begin{equation*}\bigcup_{k=1}^K[t_{k-1}, t_k]=[0,T], \quad \text{ with } t_0=0, \ t_K=T, \text{ and }t_{k-1}<t_{k}, \ \forall k\in\{1, \ldots, K\},\end{equation*}

and then concatenating the accepted skeletons in each of them. To solve the problem when $t\to 0$, we follow the recommendation in [Reference Jenkins and Spanò32], and whenever $t<0.05$, resort to the approximation in (3.8).

While asymptotically the algorithm’s growth rate does not depend on the number of loci, it is worth mentioning that with increasing L the acceptance probability decreases, as there are a larger number of skeletons that need to be accepted simultaneously, which naturally affects the algorithm’s feasibility. This will be clearer in the next section, where simulation results for examples with $L=2$ and $L=4$ are provided. As expected, the acceptance probability also decreases whenever the target diffusion differs more from the neutral Wright–Fisher candidate. This is exemplified in the next section, where results are shown for two coupled Wright–Fisher models with the same number of loci and mutation parameters, but different selective pairwise interactions.

5. Numerical experiments

In the following, several implementations of Algorithm 4 are shown along with their simulation results. The examples below represent plausible network structures present in biological applications, and whose interaction parameters are of interest. Examples of these include frequency-dependent selection in population dynamics with applications to vaccine interventions [Reference Corander11], genome-wide discovery of interdependent loci affecting antibiotic resistance [Reference Schubert, Maddamsetti, Nyman, Farhat and Marks42], co-evolution of interacting human gamete-recognition genes [Reference Rohlfs, Swanson and Weir41], and analysis of sequence data [Reference Gao, Zhou and Aurell23].

Consider the case of two loci with two allele types each, i.e., $L=2$ and $d_1=d_2=2$. A particular example with mutation, one-type allele interaction between loci, and no within-locus selection reduces (2.1) to

(5.1) \begin{equation} \begin{cases} dX_t^{11}= \alpha^{11}(X_t^{11}) dt+X_t^{11}(1-X_t^{11})hX_t^{21} dt+\sqrt{X_t^{11}(1-X_t^{11})}dB_t^1,\\[8pt] dX_t^{21}= \alpha^{21}(X_t^{21}) dt+ X_t^{21}(1-X_t^{21})hX_t^{11} dt+\sqrt{X_t^{21}(1-X_t^{21})}dB_t^2, \end{cases} \end{equation}

where the $\alpha^{i1}({\cdot})$ are as in (2.2), $B=(B_t^1,B_t^2)$ is a vector of independent Brownian motions, $H^{12}=H^{21}=\begin{pmatrix} h&\quad 0 \\ 0&\quad 0\end{pmatrix}$, and H is 0 elsewhere.

In this case, $A(X_0,X_T)\,:\!=\,h(X_T^{11}X_T^{21}-X_0^{11}X_0^{21})$ and

\begin{align*} \phi(X_t) = & \frac{1}{2}\left(h^2 \left[ (X_t^{21})^2 X_t^{11}(1-X_t^{11}) + (X_t^{11})^2 X_t^{21}(1-X_t^{21})\right]\right.{}\\[5pt] & \left.{}+ 2h \left[X_t^{21}\alpha^{11}(X_t^{11}) + X_t^{11}\alpha^{21}(X_t^{21})\right]\right), \end{align*}

and the bound constants were set to $A^+=h(1-x_0^{11}x_0^{21})$, $C^-=-\frac{h}{2}(|\theta^{1}|+|\theta^{2}|)$, and $C^+=\frac{h}{2}(\frac{h}{2}+\theta^{1}_{1}+\theta^{2}_{1})$, where $|\theta^i|=\theta_{1}^{i}+\theta_{2}^{i}$.

Table 1.

Table for 10,000 sampled paths satisfying (5.1) with $h=0.1$ and $\theta^{1}_{1}=\theta^{1}_{2}=\theta^{2}_{1}=\theta^{2}_{2}=0.01$.

The results of several simulation scenarios for sampled skeletons of the diffusion solution of (5.1) are shown in Table 1 and Table 2, which report the total length T of the considered interval [0,T], the initialization of the path $(x_0^{11}, x_0^{21})$, the average number of attempts (drawn skeletons) until acceptance, the average number of Poisson points needed until acceptance, the average number of coefficients computed in Algorithm 3, the acceptance probability, the number of approximations needed due to small ts in between drawn points of the skeleton, and the average time in seconds per accepted path.

Table 2.

Table for 10,000 sampled paths satisfying (5.1) with $h=1$ and $\theta^{1}_{1}=\theta^{1}_{2}=\theta^{2}_{1}=\theta^{2}_{2}=0.01$.

Figure 1.

Stationary distribution (left) and histogram of 10,000 samples (right) from $(X_t^{11},X_t^{21})$ satisfying (5.1) at $T=5$ with $(x_0^{11},x_0^{21})=(0.5, 0.5)$, $h=0.1$, and $\theta^{i}_j=0.01$.

As shown in Table 1 and Table 2, the average number of coefficients needed is larger for shorter intervals, where the sampled Poisson points are more likely to be close to each other ($t\to 0$), than for longer intervals, as expected from the results presented in Proposition 4.1. Also, the acceptance probabilities when $h=1$ drop to around half compared with the simulations when $h=0.1$. This is also expected, as the model with larger pairwise interaction parameter differs more from the candidate paths, so acceptance of the candidate becomes harder. This is also reflected in the average number of attempts, needed Poisson points, and coefficients, which increase consistently in the case $h=1$. Running time also increases with increasing T. In the case $h=1$ and $T=5$ the total running times became prohibitive.

In order to establish the correctness of Algorithm 4 and with the aim of providing a qualitative comparison, sampled paths of the solution of (5.1) at large T are compared with the corresponding stationary density, and this is done for different mutation parameters (see Figure 1 and Figure 2) with satisfactory results. Note that the stationary density for (5.1) can be explicitly computed, and its normalizing constant reads as follows, see [Reference Aurell, Ekeberg and Koski1]:

\begin{equation*}Z=\frac{\Gamma(\theta^2_1)\Gamma(\theta^2_2)}{\Gamma(\theta^2_1+ \theta^2_2)}\sum_{n=0}^{\infty}\frac{(\theta^2_1)^{(n)}(2h)^n}{|\theta^2|^{(n)}n!} \frac{\Gamma(\theta^1_2)}{\Gamma(|\theta^1|+n)}\Gamma(\theta^1_1+n),\end{equation*}

where $a^{(n)}=a(a+1)\ldots(a+n-1)$. For the qualitative comparisons in Figure 1 and Figure 2, the infinite sum is truncated at a sufficiently high n so that the remainder is negligible. For the parameters in Table 1 and 2 the truncation level is set to $n=70$.

Figure 2.

Stationary distribution (left) and histogram of 10,000 samples (right) from $(X_t^{11},X_t^{21})$ satisfying (5.1) at $T=5$ with $(x_0^{11},x_0^{21})=(0.5, 0.5)$, $h=0.1$, and $\theta^{1}_j=1.2,\theta^{2}_j=0.8$.

Consider now the case of four loci with two allele types each, i.e., $L=4$ and $d_i=2$ for $i\in\{1,2,3,4\}$. A particular example with mutation, one-type allele interactions between loci, and no within-locus selection reduces (2.1) to

(5.2) \begin{equation} \begin{cases} dX_t^{11}=\alpha^{11}(X_t^{11}) dt+X_t^{11}(1-X_t^{11})(h_1X_t^{21}+h_2X_t^{31}h_3X_t^{41}) dt+{}\\[3pt] \qquad\quad +\sqrt{X_t^{11}(1-X_t^{11})}dB_t^1, \\[3pt] dX_t^{21}=\alpha^{21}(X_t^{21}) dt+ X_t^{21}(1-X_t^{21})h_1X_t^{11} dt+\sqrt{X_t^{21}(1-X_t^{21})}dB_t^2,\\[3pt] dX_t^{31}=\alpha^{31}(X_t^{31}) dt+ X_t^{31}(1-X_t^{31})h_2X_t^{11} dt+\sqrt{X_t^{31}(1-X_t^{31})}dB_t^3,\\[3pt] dX_t^{41}= \alpha^{41}(X_t^{41}) dt+ X_t^{41}(1-X_t^{41})h_3X_t^{11} dt+\sqrt{X_t^{41}(1-X_t^{41})}dB_t^4,\\ \end{cases}\end{equation}

where the $\alpha^{i1}({\cdot})$ are as in (2.2), $B=(B_t^1,B_t^2,B_t^3,B_t^4)$ is a vector of independent Brownian motions,

\begin{equation*}H^{12}=H^{21}=\begin{pmatrix} h_1&\quad 0 \\[3pt] 0&\quad 0\end{pmatrix},\, H^{13}=H^{31}=\begin{pmatrix} h_2&\quad 0 \\[3pt] 0&\quad 0\end{pmatrix}, \text{ and }H^{14}=H^{41}=\begin{pmatrix} h_3&\quad 0 \\[3pt] 0&\quad 0\end{pmatrix},\end{equation*}

and H is 0 elsewhere. In this case,

\begin{equation*} A(X_0,X_T)=h_1(X_T^{11}X_T^{21}-X_0^{11}X_0^{21}]+h_2[X_T^{11}X_T^{31}-X_0^{11}X_0^{31}] +h_3[X_T^{11}X_T^{41}-X_0^{11}X_0^{41})\end{equation*}

and

\begin{align*} \phi(X_t) = & \frac{1}{2}([(h_1)^2X_t^{21}(1-X_t^{21}) +(h_2)^2X_t^{31}(1-X_t^{31})+(h_3)^2X_t^{41}(1-X_t^{41})](X_t^{11})^2\\[2pt] & +(h_1X_t^{21}+h_2X_t^{31}+h_3X_t^{41})^2X_t^{11}(1-X_t^{11})+ 2[(h_1X_t^{21}+h_2X_t^{31}+h_3X_t^{41})\alpha^{11}(X_t^{11})\\[2pt] & +(h_1\alpha^{21}(X_t^{21})+h_2\alpha^{31}(X_t^{31})+h_3\alpha^{41}(X_t^{41}))X_t^{11}]), \end{align*}

and the bound constants were set to

\begin{align*} &A^+=h_1(1-x_0^{11}x_0^{21})+h_2(1-x_0^{11}x_0^{31})+h_3(1-x_0^{11}x_0^{41}),\\& C^-=-\frac{1}{2}(|\theta^1|(h_1+h_2+h_3)+|\theta^{2}|h_1+|\theta^{3}|h_2+|\theta^{4}|h_3), \text{ and }\\& C^+=\frac{1}{2}(\theta^{1}_{1}(h_1+h_2+h_3)+\frac{(h_1+h_2+h_3)^2+(h_1)^2+(h_2)^2+(h_3)^2}{4}+h_1\theta^{1}_{2}+h_2\theta^{1}_{3}+h_3\theta^{1}_{4}).\end{align*}

Similarly to the previous example with $L=2$ and interaction parameter $h=1$, the model in (5.2) differs more from the candidate process than, say, a model with only one pairwise interaction parameter (that is, a model with, for example, $h_2=h_3=0$). This is clearly reflected in the low average acceptance probabilities, or, equivalently, in the average number of attempts or simulated Poisson points needed until acceptance; see Table 3. Moreover, simulating from model (5.2) implicitly requires the simultaneous acceptance of four candidate paths, which makes acceptance more difficult. Nonetheless, it is still feasible to use Algorithm 4 in these scenarios, as other approximate simulation strategies would be affected by similar problems.

Table 3.

Table for 10,000 sampled paths with $h_1=0.1$, $h_2=0.15$, $h_3=0.2$, and $\theta^{i}_{j}=0.2$.

6. Simulation of multidimensional neutral Wright–Fisher bridges

To complete our simulation scheme, this section presents an exact simulation technique for sampling from neutral $(d-1)$-dimensional Wright–Fisher bridges. As mentioned before, once Algorithm 4 recovers a skeleton of the desired coupled Wright–Fisher diffusion, the remainder of the path can be filled by sampling from the corresponding neutral Wright–Fisher bridges, with no further reference to the target distribution needed; see, for example, [Reference Beskos, Papaspiliopoulos and Roberts3].

Consider a $(d-1)$-dimensional Wright–Fisher bridge, between x at time 0 and z at time t. Its transition density is given by

(6.1) \begin{equation}g_{z,t}(x,y;s)=\frac{g(x,y;s)g(y,z;t-s)}{g(x,z;t)}, \quad 0<s<t,\end{equation}

where $g(\cdot,\cdot;\cdot)$ is as in (3.2), see [Reference Fitzsimmons, Pitman and Yor22]: The precise eigenfunction expansion for $g_{z,t}(x,\cdot;s)$ is provided in the following proposition, whose proof can be found in Section 7.

Proposition 6.1. Let $g_{z,t}(x,\cdot;s)$ be as in (6.1), the transition density function of a $(d-1)$-dimensional Wright–Fisher bridge. Then its eigenfunction expansion reads

\begin{equation*} g_{z,t}(x,y;s)=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty} \sum_{\underset{|l|=m}{l}}\sum_{\underset{|r|=n}{r}} p_{m,n,l,r}^{(x,z,s,t,\theta)}\mathcal{D}_{\theta+l+r}(y), \end{equation*}

with

\begin{equation*} p_{m,n,l,r}^{(x,z,s,t,\theta)}=\frac{q_m^{\theta}(s)q_n^{\theta}(t-s)}{g(x,z;t)} \mathcal{M}_{m,x}(l)\mathcal{D}_{\theta+r}(z) \mathcal{DM}_{\theta+l;n}(r), \end{equation*}

where $\mathcal{DM}_{\theta+l;n}({\cdot})$ denotes the PMF of a Dirichlet–multinomial random variable, with

\begin{equation*} \mathcal{DM}_{\theta+l;n}(r)=\frac{n!\Gamma(|\theta+l|)} {\Gamma(|\theta+l+r|)}\prod_{j=1}^{d}\frac{\Gamma(\theta_j+l_j+r_j)} {l_j!\Gamma(\theta_j+l_j)}. \end{equation*}

Following the result in Proposition 6.1, a sampling scheme for $g_{z,t}(x,y;s)$ is provided in Algorithm 5.

Similarly to Step 1 in Algorithm 2, sampling exactly from the discrete random variable with PMF $p_{m,n,l,r}^{(x,z,s,t,\theta)}$ is not straightforward. However, given the results obtained so far, it only remains to find how to evaluate $g(x,z;t)$ without approximation error, and the sampling strategy will be complete. As pointed out in [Reference Jenkins and Spanò32] for the one-dimensional case, note that the problem of evaluating $g(x,z;t)$ at x and z is different from that of sampling from it.

By (3.2) and (3.3), one obtains

(6.2) \begin{equation} g(x,z;t)= \sum_{m=0}^{\infty} \sum^\infty_{i=0}(-1)^{i} c_{m+i,m}^{(x,z,t,\theta)}(m)= \sum_{m=0}^{\infty} \sum^\infty_{i=0}(-1)^{i}b_{m+i}^{(t,\theta)}(m)\mathbb{E}[\mathcal{D}_{\theta+L_m}(z)], \end{equation}

where $L_m\sim\text{Multinomial}(m,x)$.

As in Section 3, the aim is to find monotonically converging bounds on $p_{m,n,l,r}^{(x,z,s,t,\theta)}$ so that the alternating series method can be applied.

Let

\begin{equation*}d_{2m}=\sum_{i=0}^{m}c_{m+i,m-i}^{(x,z,t,\theta)}(m), \qquad d_{2m+1}=\sum_{i=0}^{m}c_{m+1+i,m-i}^{(x,z,t,\theta)}(m),\end{equation*}

for $m=0,1,\ldots,$ which, rearranging the terms in (6.2), gives the alternating series

(6.3) \begin{equation}g(x,z;t)=\sum_{m=0}^{\infty} (d_{2m} - d_{2m+1})= d_0-d_1+d_2-\ldots.\end{equation}

Indeed, it can be proved that the terms $(d_i)_{i\geq 0}$ are monotonically decreasing from a certain threshold that is characterized in the following results, which is a condition required to apply the alternating series method.

First, note that the strategy presented in [Reference Jenkins and Spanò32] for the one-dimensional case applies here almost without change, with the exception of the terms involving $\mathbb{E}[\mathcal{D}_{\theta+L_m}(z)]$ which by an analogous (generalized) argument can be shown to decrease in m, as shown in the following lemma, proved in Section 7.

Lemma 6.1. Let $L_m\sim{Multinomial}(m,x)$. Then for all $m\in\mathbb{N}$,

\begin{equation*} \mathbb{E}[\mathcal{D}_{\theta+L_{m+1}}(z)]\leq \tilde{K}^{(\theta,x,z)} \mathbb{E}[\mathcal{D}_{\theta+L_m}(z)], \end{equation*}

where

(6.4) \begin{align} \tilde{K}^{(\theta,x,z)} = & \left(\dfrac{|\theta|}{\theta_d}\left(1-\sum_{j=1}^{d-1}z_j\right)\vee \dfrac{2(1+|\theta|)}{1-\sum_{j=1}^{d-1}z_j}\right)\left(1-\sum_{j=1}^{d-1}x_j\right)\nonumber\\[4pt] & +\sum_{j=1}^{d-1}\left(\dfrac{|\theta|}{\theta_j}z_j\vee \dfrac{2(1+|\theta|)}{z_j}\right)x_j.\end{align}

Now, similarly to (3.5), the following bound is defined:

(6.5) \begin{equation}E^{(t,\theta)}=\inf\left\{m\geq 0 :\ 2j\geq C_{m-j}^{(t,\theta)} \text{ for all } j=0,\ldots, m\right\}.\end{equation}

This is used in Proposition 6.2 below, which fully characterizes the bound on m. The proof is again deferred to Section 7.

Proposition 6.2. Let the sequence $(d_i)_{i\geq 0}$ be as defined in (6.3), and consider the bounds $E^{(t,\theta)}$, $D_{\epsilon}^{(t,\theta)}$, and $\tilde{K}^{(\theta,x,z)}$ as in (6.5), (3.6), and (6.4) respectively. Then, for $\epsilon\in(0,1)$,

\begin{equation*} d_{2m+2}<d_{2m+1}<d_{2m} \end{equation*}

for all $m\geq E^{(t,\theta)}\vee D_{\epsilon}^{(t,\theta)}\vee 2\tilde{K}^{(\theta,x,z)}/\epsilon$.

Once the bound in Proposition 6.2 is established, exact simulation of $(d-1)$-dimensional Wright–Fisher bridges ($d>2$) is possible by setting

\begin{equation*}\tilde{F}_{m,n,l,r}^{(s,t,\theta)}\,:\!=\, C_{m}^{(s,\theta)} \vee C_{n}^{(t-s,\theta)} \vee E^{(t,\theta)}\vee D_{\epsilon}^{(t,\theta)}\vee 2\tilde{K}^{(\theta,x,z)}/\epsilon,\end{equation*}

which for $2u\geq \tilde{F}_{m,n,l,r}^{(s,t,\theta)}$ provides the monotonically converging bounds

\begin{equation*}e_{m,n,l,r}(2u+1)<e_{m,n,l,r}(2u+3)<p_{m,n,l,r}^{(x,z,s,t,\theta)}< e_{m,n,l,r}(2u+2)<e_{m,n,l,r}(2u),\end{equation*}

where

\begin{equation*}e_{m,n,l,r}(u)= \frac{\sum^{u}_{i=0}(-1)^{i}b_{m+i}^{(s,\theta)}(m)\sum^{u}_{i=0}(-1)^{i}b_{n+i}^{(t-s,\theta)}}{\sum^{u+1}_{i=0}(-1)^{i} d_i}\mathcal{M}_{m,x}(l)\mathcal{D}_{\theta+r}(z)\mathcal{DM}_{\theta+l;n}(r);\end{equation*}

see [Reference Jenkins and Spanò32, Proposition 4].

To recover exact samples from $p_{m,n,l,r}^{(x,z,s,t,\theta)}$, consider the pairing bijective function $\Sigma:\mathbb{N}\rightarrow \mathbb{N}\times \mathbb{N}\times \mathbb{N}^d \times \mathbb{N}^d$ such that $\Sigma(j)=(m,n,l,r)$. Now, for each j there exists $v_j$, an element of $\tilde{v}\in\mathbb{R}^{J+1}$ (i.e., $\tilde{v}=\{v_j\}_{j=0}^J$), such that

\begin{equation*}R^-_{\tilde{v}}(J)\,:\!=\,\sum_{j=0}^J e_{\Sigma(j)}(2v_j+1)\leq \sum_{j=0}^J p_{\Sigma(j)}^{(x,z,s,t,\theta)}\leq \sum_{j=0}^J e_{\Sigma(j)}(2v_j)\,:\!=\,R^+_{\tilde{v}}(J),\end{equation*}

providing a setting analogous to the one presented in Section 3. The proposed exact sampling scheme can be found in [Reference Jenkins and Spanò32, Algorithm 5], which we reproduce here as Algorithm 6 for completeness.

Other approaches to exact simulation of one-dimensional Wright–Fisher bridges include that recently proposed in [Reference Griffiths, Jenkins and Spanò27], which restricts to the case where either $\theta_1$ or $\theta_2$ is 0, and one or both of x and z are 0, which is not applicable here.

7. Proofs

Proof of Proposition 2.1. Let $V:\mathbb{R}^n \to \mathbb{R}^{n}$ be such that $V\,:\!=\,\nabla_{x}(\overline{V}\circ f)(x)$, where we recall that

\begin{equation*}f^{ik}(x)=\overline{x}^{ik}=\bigg\{\begin{array}{cl} x^{ij},& j,k=\{1,\ldots, d_i-1\},\\[5pt] 1-\sum_{j=1}^{d_i-1}x^{ij},& k=d_i.\\\end{array}\end{equation*}

Then for all $i\in\{1,\ldots, L\}$, $j\in\{1,\ldots,d_i-1\}$, we have

\begin{equation*}\begin{array}{ll}V^{ij}(x)&=\displaystyle\frac{\partial (\overline{V}\circ f)(x)}{\partial x^{ij}}=\displaystyle\frac{\partial}{\partial \overline{x}^{ij}}\overline{V}(\overline{x})-\frac{\partial}{\partial \overline{x}^{id_i}}\overline{V}(\overline{x})\\[8pt] &=\displaystyle s^{ij} + \sum_{l=1}^{L} \sum_{k=1}^{d_l} h^{il}_{jk} \overline{x}^{lk}- s^{id_i}-\sum_{l=1}^{L}\sum_{k=1}^{d_l} h^{il}_{d_ik}\overline{x}^{lk},\\[8pt] &=\displaystyle s^{ij}- s^{id_i} + \sum_{l=1}^{L} \bigg(h^{il}_{jd_l}-h^{il}_{d_id_l} +\sum_{k=1}^{d_l-1} (h^{il}_{jk}-h^{il}_{jd_l}-h^{il}_{d_ik}+h^{il}_{d_id_l}) x^{lk}\bigg)\\[8pt] &=\displaystyle K_s^{ij}+\sum_{l=1}^{L}\bigg(K_l^{ij} +\sum_{k=1}^{d_l-1} K_{lk}^{ij} x^{lk}\bigg)=\displaystyle K_s^{ij}+\sum_{\underset{l\neq i}{l=1}}^{L}\bigg(K_l^{ij} +\sum_{k=1}^{d_l-1} K_{lk}^{ij} x^{lk}\bigg),\end{array}\end{equation*}

where the last equality holds because whenever $i=l$, all entries of the blocks $H^{ii}$ are 0.

Proof of Theorem 4.1. By definition of $A({\cdot})$, $\phi({\cdot})$, and $X_t$,

(7.1) \begin{align} &\exp\left\{A(X_0,X_T)\right\}\exp\left\{-\int_0^T \phi(X_t) dt\right\}=\nonumber\\[4pt] &=\exp\left\{\int_0^T V(X_t)\cdot dX_t -\int_0^T\dfrac{1}{2}\left[ (V(X_t))^T D(X_t) V(X_t) +2 (V(X_t))^T \alpha(X_t)\right]dt\right\}\nonumber\\[4pt] &=\exp\left\{\int_0^T D^{\frac{1}{2}}(X_t)V(X_t) \cdot dB_t -\int_0^T\dfrac{1}{2}\left[ (V(X_t))^T D(X_t) V(X_t) \right]dt\right\}.\end{align}

Because $V({\cdot})$ and $D({\cdot})$ are continuous on $[0,1]^n$, there exist constants $C^-$ and $C^+$ such that

\begin{equation*}C^-\leq\dfrac{1}{2} \displaystyle V(X_t)^T D(X_t)V(X_t) \leq C^+, \quad \text{almost surely.}\end{equation*}

Consequently,

\begin{equation*}\displaystyle\int_0^T\dfrac{1}{2} (V(X_t))^T D(X_t) V(X_t) dt<\infty, \quad \text{almost surely,}\end{equation*}

so Novikov’s condition is fulfilled and (7.1) can be identified as a Girsanov transformation [Reference Karatzas and Shreve34] with Girsanov kernel $(V(X_t))^T D^{\frac{1}{2}}(X_t)$.

Let $\mathbb{Q}$ be the probability measure with

\begin{equation*} \frac{d\mathbb{Q}}{d\mathbb{WF}L_{\alpha,x_0}} =\exp\left\{A(X_0,X_T)\right\}\exp\left\{-\int_0^T \phi(X_t) dt\right\}.\end{equation*}

It follows that the law of X under $\mathbb{Q}$ coincides with $\mathbb{CWF}_{\alpha,G,x_0}$. Indeed, by Girsanov’s theorem,

\begin{equation*} \tilde{B}_t = B_t - \int_0^t D^{\frac{1}{2}}(X_s)V(X_s) ds\end{equation*}

is a $\mathbb{Q}$-Brownian motion and

\begin{equation*} dX_t = [\alpha(X_t) + D(X_t)V(X_t)]dt+ D^{\frac{1}{2}}(X_t)d\tilde{B}_t \\ =[\alpha(X_t) + G(X_t)] dt+ D^{\frac{1}{2}}(X_t)d\tilde{B}_t.\end{equation*}

Now, by Proposition 2.1, and for $K_s,K_l, K_{lk}\in\mathbb{R}^n$,

\begin{align*} &\displaystyle\int_0^T V(X_t)\cdot dX_t =\displaystyle\int_0^T \displaystyle \bigg( K_s+\sum_{l=1}^{L}\bigg(K_l +\sum_{k=1}^{d_l-1} K_{lk}X_t^{lk}\bigg)\bigg) \cdot dX_t\\ &\quad =\displaystyle \int_0^T K_s\cdot dX_t +\int_0^T\sum_{l=1}^{L} K_l\cdot dX_t +\int_0^T\sum_{l=1}^{L}\sum_{k=1}^{d_l-1} K_{lk} X_t^{lk} \cdot dX_t\\ &\quad =\displaystyle \sum_{i=1}^{L}\sum_{j=1}^{d_i-1}\Bigg(\int_0^T K^{ij}_s dX^{ij}_t +\int_0^T \sum_{l=1}^{L}K^{ij}_l dX^{ij}_t +\displaystyle\int_0^T\sum_{l=i+1}^{L}\sum_{k=1}^{d_l-1} K^{ij}_{lk}\ d(X_t^{lk}X^{ij}_t)\Bigg),\\ \end{align*}

where the rightmost term comes from pairing terms of the form

\begin{equation*}K^{ij}_{lk}X_t^{lk}dX^{ij}_t+K^{lk}_{ij}X_t^{ij}dX^{lk}_t=K^{ij}_{lk} d(X_t^{lk}X^{ij}_t),\end{equation*}

and recalling that $K^{ij}_{lk}=K^{lk}_{ij}$, and $i\neq l$ prevents squared terms. Hence,

\begin{align*} \int_0^T V(X_t)\cdot dX_t = & \displaystyle \sum_{i=1}^{L}\sum_{j=1}^{d_i-1} \bigg(K^{ij}_s (X^{ij}_T-X^{ij}_0) +\sum_{l=1}^{L}K^{ij}_l (X^{ij}_T-X^{ij}_0)\nonumber\\[4pt] & +\displaystyle\sum_{l=i+1}^{L}\sum_{k=1}^{d_l-1} K^{ij}_{lk}\ (X^{lk}_TX^{ij}_T-X^{lk}_0X^{ij}_0)\bigg). \end{align*}

The fact that

\begin{align*}A(X_0,X_T)\,:\!=\,\sum_{i=1}^{L}\sum_{j=1}^{d_i-1} \bigg(K^{ij}_s (X^{ij}_T-X^{ij}_0) +\sum_{l=1}^{L}K^{ij}_l (X^{ij}_T-X^{ij}_0) +\displaystyle\sum_{l=i+1}^{L}\sum_{k=1}^{d_l-1} K^{ij}_{lk}\ (X^{ij}_TX^{lk}_T-X^{ij}_0X^{lk}_0)\bigg) \end{align*}

is bounded follows immediately, concluding the proof.

Proof of Proposition 4.1. Let M(t) be the total number of coefficients that must be computed in Algorithm 3, with $t\in(0,T)$ the time distance between two sampled skeleton points. By [Reference Jenkins and Spanò32, Proposition 5(iv)], there exists a $\kappa>0$ such that $\mathbb{E}[M(t)]=o(t^{-(1+\kappa)})$ as $t\to 0$, and further random coefficients needed in Algorithm 2 do not add to the algorithms’ complexity. Similarly, although our rejection scheme uses Algorithm 3 L times,

\begin{equation*}\mathbb{E}[LM(t)]=L\mathbb{E}[M(t)]=o(t^{-(1+\kappa)}),\end{equation*}

so the algorithm’s complexity is proportional to L, but its growth rate as $t\to 0$ remains the same as with $L=1$.

Now let $\epsilon\,:\!=\,d\mathbb{CWF}_{\alpha,G,x_0}/d\mathbb{WF}L_{\alpha,x_0}$ be the acceptance–rejection probability in Algorithm 4. Because $A^-\leq A(X_T)\leq A^+$ and $C^-\leq\phi(X_t)\leq C^+$,

\begin{align*}\epsilon\propto\exp\left\{A(X_T)-A^+\right\}\exp\left\{-\int_0^T (\phi(X_t)- C^-)dt\right\} \geq \exp\{-T(C^ +-C^-) + A^--A^+\}.\end{align*}

Let D refer to the number of Poisson points needed to decide upon acceptance or rejection of a proposed path. Then, following [Reference Beskos, Papaspiliopoulos and Roberts3, Proposition 3],

\begin{equation*}\mathbb{E}[N(T)]=\mathbb{E}[D]/\epsilon=T(C^+-C^-)/\epsilon\leq T(C^+-C^-)e^{T(C^+-C^-)+A^+-A^-},\end{equation*}

where the first equality follows from considering the expectation of the sum of all drawn Poisson points $\sum_{i=1}^I D_i$ over I iterations of the algorithm until the first path is accepted. First conditioning on I and applying the law of iterated expectations, and then applying the law of total expectation, concludes the proof.

Proof of Proposition 6.1. Let $g_{z,t}(x,y;s)$ be the transition density function of a $(d-1)$-dimensional Wright–Fisher bridge, between x at time 0 and z at time t. By (6.1) and (3.2),

\begin{align*}g_{z,t}(x,y;s)= & \frac{g(x,y;s)g(y,z;t-s)}{g(x,z;t)}\\[3pt] & =\frac{1}{g(x,z;t)}\sum_{m=0}^{\infty} q_m^{\theta}(s)\sum_{\underset{|l|=m}{l}} \mathcal{M}_{m,x}(l)\mathcal{D}_{\theta+l}(y)\sum_{n=0}^{\infty} q_n^{\theta}(t-s) \sum_{\underset{|r|=n}{r}} \mathcal{M}_{n,y}(r)\mathcal{D}_{\theta+r}(z)\\[3pt] &=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\sum_{\underset{|l|=m}{l}}\sum_{\underset{|r|=n}{r}}\frac{q_m^{\theta}(s)q_n^{\theta}(t-s)}{g(x,z;t)}\mathcal{M}_{m,x}(l)\mathcal{D}_{\theta+r}(z)\mathcal{M}_{n,y}(r)\mathcal{D}_{\theta+l}(y).\end{align*}

By definition,

\begin{align*}\mathcal{M}_{n,y}(r)\mathcal{D}_{\theta+l}(y)= & \frac{n!}{\prod_{j=1}^{d}l_j!}(1-\sum_{j=1}^{d-1}y_j)^{l_{d}}\prod_{j=1}^{d-1}y_j^{r_j}\\[4pt] & \times\frac{\Gamma(|\theta+l|)}{\prod_{j=1}^{d}\Gamma(\theta_j+l_j)} (1-\sum_{j=1}^{d-1}y_j)^{\theta_{d}+l_{d}-1}\prod_{j=1}^{d-1}y_j^{\theta_j+l_j-1}. \end{align*}

Multiplying and dividing by $\dfrac{\Gamma(|\theta+l+r|)}{\prod_{j=1}^d \Gamma(\theta_j+l_j+r_j)}$ and rearranging shows that

\begin{multline*} \frac{n!\Gamma(|\theta+l|)}{\Gamma(|\theta+l+r|)}\prod_{j=1}^{d}\frac{\Gamma(\theta_j+l_j+r_j)}{\Gamma(\theta_j+l_j)r_j!}\dfrac{\Gamma(|\theta+l+r|)}{\prod_{j=1}^d \Gamma(\theta_j+l_j+r_j)}(1-\sum_{j=1}^{d-1}y_j)^{\theta_{d}+l_{d}+r_d-1}\\\times\prod_{j=1}^{d-1}y_j^{\theta_j+l_j+r_j-1}=\mathcal{DM}_{\theta+l;n}(r)\mathcal{D}_{\theta+l+r}(y). \end{multline*}

Now, identifying the coefficients of $p_{m,n,l,r}^{(x,z,s,t,\theta)}$, the proof is complete.

Proof of Lemma 6.1. For the sake of simplicity, we will first consider the case $d=3$ and then show that by analogous arguments the results can be extended to the general $(d-1)$-dimensional case. First, note that the indexes of the sum on the right-hand side of

\begin{equation*}\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)]=\sum_{\underset{|l|=m}{l}}\mathbb{P}(L_{m}=l)\mathcal{D}_{\theta+l}(z)\end{equation*}

can be seen as being placed on the $(d-1)$-face of the d-simplex $\mathcal{S}_d=\{l\in\mathbb{R}^d\mid l_{j}\geq 0,\sum_{j=1}^{d}l_j=m\}$. For example, if $d=3$, we only need to consider indexes $l_1$ and $l_2$ (see Figure 3). Thus,

(7.2) \begin{equation}\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)]=\sum_{l_1=0}^{m}\sum_{l_2=0}^{m-l_1}\mathbb{P}(L_{m}=l)\mathcal{D}_{\theta+l}(z).\end{equation}

Figure 3.

Graphical representation of the two-dimensional case ($d=3$). Each point in the figure represents an element of the sum in (7.2). For every and up to each m, all points depicted in the plane $l_1-l_2$ (left) correspond to points depicted in the mth $(l_1-l_2-l_3)$-face (right), which in turn represent each and every one of the summands in (7.2).

Let us define the quantities

\begin{equation*}\mathcal{Q}_{m}\,:\!=\,\mathbb{P}(L_{m}=l)\mathcal{D}_{\theta+l}(z), \qquad \mathcal{Q}_{m+1}^{(j)}\,:\!=\,\mathbb{P}(L_{m+1}=\tilde{l})\mathcal{D}_{\theta+\tilde{l}}(z),\end{equation*}

where $\tilde{l}\in\mathbb{R}^d$ is equal to l except in its jth component, for which $\tilde{l}_j=l_j+1$. As a special case, when $j=d$, one can write $\tilde{l}_d=m-\sum_{j=1}^{d-1}l_j+1$.

Then, for $l_1\leq\lfloor mz_1 \rfloor, l_2\leq\lfloor mz_2 \rfloor$,

(7.3) \begin{align}\mathcal{Q}_{m+1}^{(d)} =&\frac{m+1}{m+1-l_1-l_2}\frac{|\theta|+m}{\theta_3+m-l_1-l_2}(1-x_1-x_2)(1-z_1-z_2)\mathcal{Q}_{m}\\[4pt] &\leq \frac{m+1}{1+m(1-z_1-z_2)}\frac{|\theta|+m}{\theta_3+m(1-z_1-z_2)}(1-x_1-x_2)(1-z_1-z_2)\mathcal{Q}_{m}\nonumber\\[4pt] &\leq \left(\dfrac{|\theta|}{\theta_3}\vee \frac{m^2+m(1+|\theta|)+|\theta|}{m^2(1-z_1-z_2)^2}\right)(1-x_1-x_2)(1-z_1-z_2)\mathcal{Q}_{m}\nonumber\\[4pt] &\leq \left(\dfrac{|\theta|}{\theta_3}(1-z_1-z_2)\vee \dfrac{2(1+|\theta|)}{1-z_1-z_2}\right)(1-x_1-x_2)\mathcal{Q}_{m}. \nonumber\end{align}

Here, in the second inequality, note that the function

\begin{equation*}g(m)\,:\!=\,\frac{m^2+m(1+|\theta|)+|\theta|}{m^2(1-z_1-z_2)^2}\end{equation*}

is decreasing in m, and thus for $m\geq 1$, it attains its maximum value at g(1). Then, if

\begin{equation*}f(m)\,:\!=\, \frac{m+1}{1+m(1-z_1-z_2)}\frac{|\theta|+m}{\theta_3+m(1-z_1-z_2)},\end{equation*}

one obtains $f(m)\leq f(0)\vee g(m)\leq f(0)\vee g(1)$ yielding the desired result.

Similarly, for $l_1\leq\lfloor mz_1 \rfloor, \lfloor mz_2 \rfloor\leq l_2$,

(7.4) \begin{align}\mathcal{Q}_{m+1}^{(2)} &=\frac{m+1}{l_2+1}\frac{|\theta|+m}{\theta_2+l_2}x_2z_2\mathcal{Q}_{m}\leq \frac{m+1}{mz_2+1}\frac{|\theta|+m}{\theta_2+mz_2}x_2z_2\mathcal{Q}_{m}\nonumber\\[5pt] &\leq \left(\dfrac{|\theta|}{\theta_2}\vee \frac{m^2+m(1+|\theta|)+|\theta|}{m^2z_2^2}\right)x_2z_2\mathcal{Q}_{m}\leq \left(\dfrac{|\theta|}{\theta_2}z_2\vee \dfrac{2(1+|\theta|)}{z_2}\right)x_2\mathcal{Q}_{m},\end{align}

and for $\lfloor mz_1 \rfloor\leq l_1$ and all $l_2$,

(7.5) \begin{align}\mathcal{Q}_{m+1}^{(1)} &=\frac{m+1}{l_1+1}\frac{|\theta|+m}{\theta_1+l_1}x_1z_1\mathcal{Q}_{m}\leq \frac{m+1}{mz_1+1}\frac{|\theta|+m}{\theta_1+mz_1}x_1z_1\mathcal{Q}_{m}\nonumber\\[5pt] &\leq \left(\dfrac{|\theta|}{\theta_1}\vee \frac{m^2+m(1+|\theta|)+|\theta|}{m^2z_1^2}\right)x_1z_1\mathcal{Q}_{m}\leq \left(\dfrac{|\theta|}{\theta_1}z_1\vee \dfrac{2(1+|\theta|)}{z_1}\right)x_1\mathcal{Q}_{m}.\end{align}

Combining the inequalities in (7.3), (7.4), and (7.5),

\begin{align*}\mathbb{E}[\mathcal{D}_{\theta+L_{m+1}}(z)]=&\sum_{l_1=0}^{\lfloor mz_1 \rfloor}\sum_{l_2=0}^{\lfloor mz_2 \rfloor}\mathcal{Q}_{m+1}^{(d)} + \sum_{l_1=0}^{\lfloor mz_1 \rfloor}\sum_{l_2=\lfloor mz_2\rfloor+1}^{m+1-l_1} \mathcal{Q}_{m+1}^{(2)}+ \sum_{l_1=\lfloor mz_1 \rfloor+1}^{m+1}\!\!\sum_{l_2=0}^{m+1-l_1} \mathcal{Q}_{m+1}^{(1)}\\[8pt] &\leq \left(\dfrac{|\theta|}{\theta_3}(1-z_1-z_2)\vee \dfrac{2(1+|\theta|)}{1-z_1-z_2}\right)(1-x_1-x_2)\sum_{l_1=0}^{\lfloor mz_1 \rfloor}\sum_{l_2=0}^{\lfloor mz_2 \rfloor} \mathcal{Q}_{m}\\[8pt] &+\left(\dfrac{|\theta|}{\theta_2}z_2\vee \dfrac{2(1+|\theta|)}{z_2}\right)x_2\sum_{l_1=0}^{\lfloor mz_1 \rfloor}\sum_{l_2=\lfloor mz_2 \rfloor}^{m-l_1} \mathcal{Q}_{m}\\[8pt] &+\left(\dfrac{|\theta|}{\theta_1}z_1\vee \dfrac{2(1+|\theta|)}{z_1}\right)x_1\sum_{l_1=\lfloor mz_1 \rfloor}^{m}\sum_{l_2=0}^{m-l_1} \mathcal{Q}_{m}\\[8pt] &\leq\left[\left(\dfrac{|\theta|}{\theta_3}(1-z_1-z_2)\vee \dfrac{2(1+|\theta|)}{1-z_1-z_2}\right)(1-x_1-x_2)\right.\\[8pt] &+\left.\displaystyle\sum_{j=1}^2\left(\dfrac{|\theta|}{\theta_j}z_j\vee \dfrac{2(1+|\theta|)}{z_j}\right)x_j\right]\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)],\end{align*}

where, in the first inequality, terms starting at $l_j=\lfloor mz_j \rfloor+1$ are compared one-to-one to terms starting at $l_j=\lfloor mz_j \rfloor$, and the last inequality holds after taking common factors and noting that the terms for $l_j=\lfloor mz_j \rfloor$ are bounded by both

\begin{equation*}\left(\dfrac{|\theta|}{\theta_3}(1-z_1-z_2)\vee \dfrac{2(1+|\theta|)}{1-z_1-z_2}\right)(1-x_1-x_2)\end{equation*}

and

\begin{equation*}\sum_{j=1}^2\left(\dfrac{|\theta|}{\theta_j}z_j\vee \dfrac{2(1+|\theta|)}{z_j}\right)x_j.\end{equation*}

The proof for the general $(d-1)$-dimensional case follows analogously. Consider

\begin{equation*}\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)]=\sum_{l_1=0}^{m}\sum_{l_2=0}^{m-l_1}\ldots\sum_{l_{d-1}=0}^{m-|l|_{d-2}}\mathbb{P}(L_{m}=l)\mathcal{D}_{\theta+l}(z),\end{equation*}

where $|l|_{d-2}=\sum_{j=1}^{d-2} l_j$.

Following the strategy used above, the sums can be partitioned into terms such that either

1. (a) $l_j\leq \lfloor mz_j \rfloor$ for all j;

2. (b) $\lfloor mz_j \rfloor\leq l_j$ for all $j\neq 1$, and $l_i\leq \lfloor mz_i \rfloor$ for all $i\neq j$; or

3. (c) $\lfloor mz_1 \rfloor\leq l_1$, and $l_i$ is free for all $i\neq 1$.

Note that this partition covers all (non-exclusive) combinations and includes all the elements of the sum. Now, comparing $\mathbb{E}[\mathcal{D}_{\theta+L_{m+1}}(z)]$ with $\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)]$, the bounding constants for each case are as follows:

1. (a)

   \begin{equation*}\left(\dfrac{|\theta|}{\theta_d}\left(1-\sum_{j=1}^{d-1}z_j\right)\vee \dfrac{2(1+|\theta|)}{1-\sum_{j=1}^{d-1}z_j}\right)\left(1-\sum_{j=1}^{d-1}x_j\right),\end{equation*}

2. (b)

   \begin{equation*}\left(\dfrac{|\theta|}{\theta_j}z_j\vee \dfrac{2(1+|\theta|)}{z_j}\right)x_j,\end{equation*}

3. (c)

   \begin{equation*}\left(\dfrac{|\theta|}{\theta_1}z_1\vee \dfrac{2(1+|\theta|)}{z_1}\right)x_1.\end{equation*}

This yields

\begin{equation*}\mathbb{E}[\mathcal{D}_{\theta+L_{m+1}}(z)]\leq \tilde{K}^{(\theta,x,z)} \mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)],\end{equation*}

with

\begin{align*} \qquad\qquad\qquad\tilde{K}^{(\theta,x,z)}= &\left(\dfrac{|\theta|}{\theta_d}\left(1-\sum_{j=1}^{d-1}z_j\right)\vee \dfrac{2(1+|\theta|)}{1-\sum_{j=1}^{d-1}z_j}\right)\left(1-\sum_{j=1}^{d-1}x_j\right)\\[5pt] & +\sum_{j=1}^{d-1}\left(\dfrac{|\theta|}{\theta_j}z_j\vee \dfrac{2(1+|\theta|)}{z_j}\right)x_j.\qquad\qquad\qquad\qquad\qquad\qquad\qquad\end{align*}

Proof of Proposition 6.2. The proof follows from that of [Reference Jenkins and Spanò32, Proposition 3], which is reproduced here for completeness.

The inequality $d_{2m+1}<d_{2m}$ follows because if $m\geq E^{(t,\theta)}$ then $2j\geq C_{m-j}^{(t,\theta)}$ for all $j=0,\ldots, m$, which, by Proposition 3.1, implies

\begin{equation*}b^{(t,\theta)}_{m+j+1}(m-j)<b^{(t,\theta)}_{m+j}(m-j)).\end{equation*}

Multiplying by $\mathbb{E}[\mathcal{D}_{\theta+L_{m-j}}(z)]$ and then summing over $j=0,\ldots,m$ gives

(7.6) \begin{equation}d_{2m+1}=\sum_{j=0}^m c_{m+j+1,m-j}^{(x,z,t,\theta)}<\sum_{j=0}^m c_{m+j,m-j}^{(x,z,t,\theta)}=d_{2m}.\end{equation}

Proving $d_{2m+2}<d_{2m+1}$ requires some extra steps. First, note that

\begin{equation*}d_{2m+2}=d_{2(m+1)}=\sum_{r=0}^{m+1} c_{m+1+r,m+1-r}^{(x,z,t,\theta)}=\sum_{j=-1}^m c_{m+2+j,m-j}^{(x,z,t,\theta)},\end{equation*}

where $j=r-1$. Noting now that $2j+1>2j\geq C_{m-j}^{(t,\theta)} \text{ for all } j=0,\ldots, m$, which implies

\begin{equation*}b^{(t,\theta)}_{m+j+2}(m-j)<b^{(t,\theta)}_{m+j+1}(m-j),\end{equation*}

and using the same argument as in (7.6) yields

\begin{equation*}\sum_{j=1}^m c_{m+j+2,m-j}^{(x,z,t,\theta)}<\sum_{j=1}^m c_{m+j+1,m-j}^{(x,z,t,\theta)},\end{equation*}

where the sum is taken only over $j=1,\ldots, m$ so that the remaining terms in $d_{2m+2}$ and $d_{2m+1}$ can be compared. Indeed, it only remains to prove that

\begin{equation*}c_{m+1,m+1}^{(x,z,t,\theta)}+c_{m+2,m}^{(x,z,t,\theta)}<c_{m+1,m}^{(x,z,t,\theta)}.\end{equation*}

Note now that

(7.7) \begin{equation}\frac{c_{k+1,m}^{(x,z,t,\theta)}}{c_{k,m}^{(x,z,t,\theta)}}=\frac{b^{(t,\theta)}_{k+1}(m)}{b^{(t,\theta)}_{k}(m)}=h_m(k)e^{(2k+|\theta|)t/2}\leq (|\theta|+2k+1)e^{(2k+|\theta|)t/2},\end{equation}

where

\begin{equation*}h_m(k)\,:\!=\,\frac{|\theta|+m+k-1}{k-m+1}\frac{|\theta|+2k+1}{|\theta|+2k-1},\end{equation*}

the second equality follows from (3.4), and the last inequality holds because $h_m(k)<h_k(k)=|\theta|+2k+1$; see the proof of [Reference Jenkins and Spanò32, Proposition 1].

Because by hypothesis $m\geq D_{\epsilon}^{(t,\theta)}$, recalling the definition of $D_{\epsilon}^{(t,\theta)}$ in (3.6) and choosing $k=m+1$ in (7.7) yields

\begin{equation*}c_{m+2,m}< (|\theta|+2k+1)e^{(2k+|\theta|)t/2}c_{m+1,m}^{(x,z,t,\theta)}<(1-\epsilon)c_{m+1,m}^{(x,z,t,\theta)}.\end{equation*}

Finally,

\begin{align*}\frac{c_{m+1,m+1}^{(x,z,t,\theta)}}{c_{m+1,m}^{(x,z,t,\theta)}}&=\frac{b^{(t,\theta)}_{m+1}(m+1)}{b^{(t,\theta)}_{m+1}(m)}\frac{\mathbb{E}[\mathcal{D}_{\theta+L_{m+1}}(z)]}{\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)]}=\frac{|\theta|+2m}{(m+1)(|\theta|+m)}\frac{\mathbb{E}[\mathcal{D}_{\theta+L_{m+1}}(z)]}{\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)]}\\&=\frac{1}{(m+1)}\left(1+\frac{m}{|\theta|+m}\right)\frac{\mathbb{E}[\mathcal{D}_{\theta+L_{m+1}}(z)]}{\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)]}<\frac{2}{(m+1)}\frac{\mathbb{E}[\mathcal{D}_{\theta+L_{m+1}}(z)]}{\mathbb{E}[\mathcal{D}_{\theta+L_{m}}(z)]}<\epsilon,\end{align*}

where the last inequality follows because $m+1>m\geq 2\tilde{K}^{(\theta,x,z)}/\epsilon$ and using Lemma 6.1, yielding

\begin{equation*}c_{m+1,m+1}^{(x,z,t,\theta)}+c_{m+2,m}<\epsilon c_{m+1,m}^{(x,z,t,\theta)}+(1-\epsilon)c_{m+1,m}^{(x,z,t,\theta)}=c_{m+1,m}^{(x,z,t,\theta)},\end{equation*}

which concludes the proof.