{'fraction_non_alphanumeric': 0.043214578694480255, 'fraction_numerical': 0.011543419664830657, 'mean_word_length': 4.739532537413822, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulating bridges using confluent diffusions

Diffusions are a fundamental class of models in many fields, including finance, engineering, and biology. Simulating diffusions is challenging as their sample paths are infinite-dimensional and their transition functions are typically intractable. In statistical settings such as parameter inference for discretely observed diffusions, we require simulation techniques for diffusions conditioned on hitting a given endpoint, which introduces further complication. In this paper we introduce a Markov chain Monte Carlo algorithm for simulating bridges of ergodic diffusions which (i) is exact in the sense that there is no discretisation error, (ii) has computational cost that is linear in the duration of the bridges, and (iii) provides bounds on local maxima and minima of the simulated trajectory. Our approach works directly on diffusion path space, by constructing a proposal (which we term a confluence) that is then corrected with an accept/reject step in a pseudo-marginal algorithm. Our method requires only the simulation of unconditioned diffusion sample paths. We apply our approach to the simulation of Langevin diffusion bridges, a practical problem arising naturally in many situations, such as statistical inference in distributed settings.


Introduction
Diffusions are popular within a variety of application areas as models for stochastic dynamical systems. These applications include, amongst others, the physical sciences [30], the life sciences [19,20], and finance [2,8,18,25,26]. Given a model for a dynamical system of interest it is natural to want to conduct inference on the model parameters. The complexity of inference for diffusions often necessitates the use of computational techniques such as Markov chain Monte Carlo (MCMC). These require sampling diffusion bridges with different parameter values, while ensuring the bridges sampled are coherent with the observed data.
To introduce diffusion bridges, first consider the following one-dimensional stochastic differential equation (SDE): We term an unconditioned diffusion as being a weak solution to (1), and if X is conditioned on a terminal value, X | X T = x T , we term it a conditioned diffusion or a diffusion bridge.
Simulating diffusions is challenging. For very simple classes of diffusions in which the transition function is known, efficient schemes for simulating trajectories of unconditioned diffusions over long (temporal) durations (T 1) are available-one can simply apply the strong Markov property and iteratively simulate partial trajectories. As such, the computational cost in this setting scales linearly in the number of iterations. However, in general the transition function for the majority of interesting classes of diffusions is intractable. For diffusion bridges the conditioning precludes the iterative chaining of partial trajectories, and so the scaling with the number of iterations is more complicated.
Existing approaches for diffusion bridge simulation fall into two broad categories: methods based on simulating from time-discretisations (and thus approximations) of (1), and those which rely on constructions directly on diffusion path-space for (1) (which we term path-space samplers). The latter eliminates any discretisation error, but usually at a cost of introducing technical restrictions on b and σ. In both cases various embedded Monte Carlo algorithms can be used (for instance, importance and rejection sampling or MCMC), with repeated simulation from an appropriately chosen proposal distribution.
A key consideration is to ensure that the proposal and target bridge are well matched, while simultaneously accounting for the computational complexities of the proposal. One of the earliest breakthroughs was made by Pedersen [29], who used a time-discretised realisation of the unconditioned diffusion (1) as a proposal. However, this approach suffers from degeneracy between the target and proposal distribution as increasingly fine time-discretisations are considered. Further work by [16,17,20,28] attempt to directly resolve this, with other approaches such as 'pulling' drifts [16], guided and residual proposals [36,38,39], and the use of stochastic partial differential equations [22], being similarly motivated. However, these approaches become increasingly computationally expensive with the fine resolutions required in settings where controlling or removing the bias is important, for instance in inference for partially observed diffusions [35].
A conceptually novel alternative for simulating bridges of ergodic diffusions using a discretised approach was recently proposed by Bladt & Sørensen and subsequently extended [10,11,9]. The proposal is formed by simulating only unconditioned diffusions-one whose start points coincides with that of the desired bridge, another whose start point coincides with the end point of the bridge-these are spliced together at a well specified crossing time to form a single proposal path, which is then either accepted or rejected by means of a stochastic number of auxiliary unconditioned diffusions. The main strength of this approach is that only unconditioned diffusions are simulated, and so the desirable linear-in-time computational cost of such simulation is retained. This is in contrast to previously proposed approaches where the computational cost can grow exponentially in T . Termed the simple diffusion bridge approach, this recent innovation forms one key ingredient of the work we introduce in this paper.
Alternatively, path-space rejection samplers allow sample paths of both unconditioned and conditioned diffusions to be simulated at finite collections of time points, without any approximation error. Following the initial work of [3], there has been considerable interest in developing unbiased estimators for classes of stochastic differential equations (including [5,12,14,23,32,33]). Methodologically these approaches follow a common form: sample paths are drawn from a target measure by means of drawing sample paths from an equivalent proposal measure in which access to the conditioned transition density is available, and are accepted with probability proportional to the Radon-Nikodým derivative of the two measures. Path-space rejection sampling approaches have been applied to Bayesian inference for partiallyobserved diffusions [4,6,37]. However, a limitation is that the computational cost increases exponentially with bridge duration, and so are unsuitable in many applications.
It would be appealing to incorporate within the simple diffusion bridge framework of Bladt & Sørensen the discretisation free path-space rejection samplers for unconditioned ergodic diffusions-this would reduce the computational cost of path-space rejection sampling approaches from exponential to linear with respect to bridge duration. However, the difficulty with a naïve embedding of these methodologies is it requires the identification of the joint crossing times of a collection of diffusions-both for the spliced proposal, and the auxiliary trajectories used in determining acceptance. Instead, we resolve this issue by adapting ideas from -strong simulation within path-space rejection sampling [7,32]. In particular, diffusion trajectories are simulated together with upper and lower convergent bounding processes which enfold (almost surely) the trajectory. We further develop methodology within this framework allowing us to unbiasedly estimate the probability that two such diffusion paths intersect. The resulting methodology, which we term the confluent diffusion bridge approach, then allows us to overcome the challenge of embedding path-space rejection sampling within the simple diffusion bridge framework, and so attains the linear computational cost of the latter. We find practical application of our work motivated by Monte Carlo Fusion [15], by considering the simulation of Langevin diffusion bridges in which a mixture of Gaussian distributions is the invariant measure of the unconditioned diffusion.
We begin in Section 2 by outlining the simple diffusion bridge sampler of Bladt & Sørensen [9], and then in Section 3 the pertinent aspects of path-space rejection sampling and ε-strong simulation. These methodologies together with our approach in Section 4 form the basis of the confluent diffusion bridge (CDB) sampler introduced in this paper. In Section 5 we apply the CDB to numerical examples: in Section 5.1, we investigate the simulation of Ornstein-Uhlenbeck diffusion bridges with increasing duration, showing that the CDB sampler attains a linear computational scaling as opposed to the exponential scaling of competing methodologies; in Section 5.2 we consider the more complicated example of the Langevin diffusion bridge, which has a number of practical statistical applications. Finally, in Section 6 we provide a brief discussion on applications and extensions of our work.

Simple diffusion bridges
The simple diffusion bridge sampler introduced by Bladt & Sørensen [9], obtains approximate draws from the conditioned diffusion bridge of (1), with law which we denote P (x 0 ,x T ,T ) . For clarity of exposition, we denote P (x 0 ,T ) as the unconditioned law. Broadly, the approach is an MCMC scheme in which a proposal bridge is obtained via rejection sampling, and then a Metropolis-Hastings correction is made to account for the discrepancy between the law of the proposal and the law of the target. The novelty of the method is the construction of the proposal [9,11], which we will now summarise.
← − X (2) is further defined as T −t , t ∈ [0, T ], and referred to as the time-reversed diffusion. Let τ (Z) := inf{0 ≤ t ≤ T |X t } and: Then it can be shown that the following equivalence in distribution holds [11,Theorem 1]: where X is a solution to (1) independent of X (3) . The law of Z|{τ (Z) ≤ T } is denoted by Z (x 0 ,x T ,T ) , and is the law of the simple diffusion bridge proposal. Equation (3) shows that Z (x 0 ,x T ,T ) is equal to the law of the diffusion bridge conditioned on being hit by an independent auxiliary diffusion started at the stationary distribution ν. The methodology for this is summarised in Algorithm 1.

Algorithm 1 Simple Diffusion Bridge sampler -Proposal
Input: Moreover, let A represent the set of pairs of intersecting functions: Let X denote the target diffusion bridge (independent from X (3) ) and define: Then as established by [9, Theorem 2.1] which gives an explicit connection between Z (x 0 ,x T ,T ) and P (x 0 ,x T ,T ) , and provides the correction term in an independence sampler. More precisely, denoting by Z (n) the sample path kept at the n th iteration, we sample a new path Z ∼ Z (x 0 ,x T ,T ) and set Z (n+1) ← Z with probability: where we suppress (x 0 ,x T ,T ) from the superscripts of the laws for the sake of clarity. Otherwise set Z (n+1) ← Z (n) . Although 1/π T (Z) cannot be computed in closed form, it is possible to construct a positive, unbiased estimator. This is sufficient for implementing the algorithm [1]: by substituting the likelihood for a positive, unbiased estimator of it, the invariant distribution of the Markov Chain is unaltered. For a proposal sample Z ∈ C([0, T ]), one such estimator is with X (3,i) , i ∈ N + , i.i.d. copies of the auxiliary diffusion X (3) . Algorithm 2 provides a summary of the full pseudo-marginal independence sampler (including the embedded proposal mechanism detailed in Algorithm 1). We caution that Algorithm 2 is not geometrically ergodic [1, Theorem 8]. T ← T + 1

4:
Sample v ∼ ν 5: ← − X (2) , and X (3) is linear in T , then the overall algorithmic complexity of Algorithm 2 is also linear in T [9, Section 2.2]; [11]. In practice, and rather surprisingly, the contribution to the cost from the geometric number of trials of the event {τ (aux) < ∞} will typically decrease with increasing T , as the intersection of X (3) and Z becomes more likely over bridges of longer duration.
For reasons of computational tractability, the implementation of the simple diffusion bridge sampler presented by Bladt & Sørensen [9] employed a number of approximations. In particular, to simulate X (i) , i = 1, 2, 3 a discretisation approach based on a stochastic Taylor approximation [24] was taken. In addition, a further pragmatic heuristic was used to define the required intersection times of the unconditioned diffusions. The need for similar approximations arises in the independence sampler correction. Consequently these approximations in the simple diffusion bridge sampler propagate, and the desired law P (x 0 ,x T ,T ) is not recovered exactly. Furthermore, the accuracy of the achieved sampler is difficult to understand. It is these approximations which (for a class of diffusions we introduce fully later) we want to remove from this scheme. Full details are deferred to Section 4.

Exact and epsilon-strong simulation of unconditioned diffusions
Two sources of error arising in the simple diffusion bridges approach outlined in Section 2 are the discretisation errors associated with the simulation of the required diffusion sample paths (X (i) , i ∈ {1, 2, 3}), and the determination of the crossing times of the diffusion sample paths. The approach we outline to resolve these errors in Section 4 exploits the availability of methodology to simulate diffusions without time-discretisation error. One class of such methods, which we term ε-strong simulation, is appropriate here because it provides auxiliary information about the diffusion including a random collection of compact time-space regions over which the diffusion locally evolves. The flexibility of being able to simulate a trajectory with enfolding local bounds, and to any desired resolution, is useful when determining whether multiple such trajectories intersect, as indicated by the methodology in Section 2 and more fully detailed in Section 4. For simplicity in this section we will consider a single trajectory, X, and denote by (X ↑ (n) : n = 1, . . . ) and (X ↓ (n) : n = 1, . . . ) sequences of upper and lower bounding processes respectively, which converge to X as the 'resolution index' n → ∞. In particular, we require we have almost surely for all t ∈ [0, T ], It is possible to do this by means of Pollock et al. [32,Sections 3.2 & 5], which we will outline in the remainder of this section. We first describe how to simulate exactly (in the sense that it is without time-discretisation error) the trajectory of a class of diffusions at any finite collection of time points, by means of rejection and using a Brownian motion proposal. We then consider how to find appropriate upper and lower bounding processes of the type described in (6).
For Brownian motion to be an appropriate proposal, we first consider a transformation of (1) to one with unit volatility; this is always possible for a scalar diffusion using the Lamperti transformation [24,§IV.4]. In particular, if then Y t := η(X t ) solves the following SDE: where As such we will assume without loss of generality that the diffusion coefficient σ in (1) is identically equal to 1. We impose a number of other standard regularity conditions on α as given in [4, page 1080], which includes that α is continuously differentiable. Under these conditions, by application of [32,Algorithm 4] it is possible to simulate finite dimensional representations of diffusion sample paths, such as those solving (8), exactly (without discretisation error). In particular, the output of the algorithm is a skeleton, where κ is a Poisson distributed random variable determining the number of uniformly distributed time points, χ 1 , . . . , χ κ , between the end-points (χ 0 := 0 and χ κ+1 := T ), at which the trajectory is revealed. R and [U ↓ X[χ j ,χ j+1 ] , U ↑ X[χ j ,χ j+1 ] ] respectively). For the purposes of the methodology developed in this paper we require only L ↓ X[χ j ,χ j+1 ] and U ↑ X[χ j ,χ j+1 ] (which determines a compact interval which constrains the trajectory).
By means of [32,Algorithm 22] it is possible to retrospectively reveal the realised sample path at any additional time point, conditional on the skeleton (9), noting that the unrealised trajectory has the law of a Brownian bridge between each pair of skeleton points conditioned on the associated layer. It is then possible to augment the skeleton with additional layer information to account for the fact that an existing layer has been bisected by the new time point. In the later sections of this paper we will consider revealing the path, and augmenting the skeleton, at multiple additional time points. This can be done in an iterative manner, at a resolution sufficient for the user (for instance, some -tolerance). An illustrative example of a realised skeleton of a trajectory following from the entire -strong simulation approach can be found in Figure 1.

Confluent diffusion bridges
In this section we outline our confluent diffusion bridge approach, which incorporates a discretisation free path-space rejection sampler for an unconditioned ergodic diffusion using -strong simulation as described in Section 3, within the simple diffusion bridge framework of Bladt & Sørensen as outlined in Section 2. Here we carry over the assumptions from both of these areas: in particular, from Section 2 we assume that the unconditioned diffusion is ergodic with unique stationary distribution ν, and from Section 3 we assume the regularity conditions on α as described in [4, page 1080].

Sampling proposals
Our aim is to sample from Z (x 0 ,x T ,T ) , avoiding the discretisation errors in the simple diffusion bridge approach, by using the exact and ε-strong simulation approach outlined in Section 3. Recalling the construction of the proposal diffusion bridge Z as given in (2) in the simple diffusion bridge approach, it is itself composed of two unconditioned diffusions: A forward diffusion, X (1) , initialised at the left hand end-point x 0 , to some time τ (Z) ; and a time-reversed diffusion, ← − X (2) , from τ (Z) to the right hand end-point t } is simply the first intersection time of X (1) and ← − X (2) . A rejection sampling scheme to sample from Z (x 0 ,x T ,T ) would simply sample pairs of candidate paths X (1) and ← − X (2) until a pair is obtained that intersects (i.e. τ (Z) < ∞). Figure 2 illustrates the procedure graphically.
We begin by simulating skeletons for X (1) and X (2) (denoted by S(X (1) ) and S(X (2) ) respectively), which immediately provides the skeleton S( T −t (as shown in Figure 2(a)). These skeletons are then 'cross-populated', by revealing the trajectories (including layer information) at additional time points, as described in Section 3. In particular, S(X (1) ) is cross-populated by the temporal points {T − χ (2) j } κ (2) j=1 , and S( ← − X (2) ) is cross-populated by the temporal points {χ (1) j } κ (1) j=1 . Each resulting skeleton is thus revealed on the common grid (where κ := κ (1) + κ (2) ) as summarised by Figure 2(b). Next, we must determine τ (Z) for the two diffusions X (1) and ← − X (2) , in order to determine whether the candidate for Z (x 0 ,x T ,T ) is to be accepted. Cross-populating the skeletons as described above so they share a common time grid, which for simplicity we assume to be sorted in increasing order, we can simply sequentially consider the intervals [χ j , χ j+1 ], j = 0, . . . , κ,   (Figure 1(a)) and at increasing resolutions (Figures 1(b) -1(d)). The skeleton at a given resolution is given by the trajectory's temporal location (black dots) at a collection of times (black dashed vertical lines) and almost sure bounds between consecutive times (black solid horizontal lines).
querying each interval with the aid of Table 1. Throughout the remainder of this section we will assume, without loss of generality, that X and U (1) ] . Considering the current interval (say, [χ j , χ j+1 ]) in our sequential querying, and referring to Table 1, if it falls into case (i) then we simply proceed to the next interval ([χ j+1 , χ j+2 ])). In case (ii) we halt our sequential evaluation, having identified that the current interval is one in which the first crossing has occurred (so τ (Z) ∈ [χ j , χ j+1 ]). In case (iii) further information is required to determined whether crossing has occurred. We propose to exploit the fact that it is possible Time X (a) Initial skeletons for X (1) and ← − X (2) .  ← − X (2) to determine if they cross, and if so the interval in which they do cross (given by blue dashed vertical lines).  to retrospectively reveal additional points from the paths underlying the existing skeletons (as described in Section 3). In particular, we begin by revealing both paths at the mid-point (χ j + χ j+1 )/2 of the current interval, and then augmenting both skeletons with corresponding layer information for interval [χ j , (χ j + χ j+1 )/2] and [(χ j + χ j+1 )/2, χ j+1 ]. We then resume our sequential evaluation (continuing with an abuse of notation to write {χ j } κ+1 j=0 for the set of time points including any additionally revealed ones), starting at the interval [χ j , (χ j + χ j+1 )/2]. If upon evaluating all intervals in the sequential manner indicated we find that all intervals fall Scenario Interpretation (i) U (1) The two paths cannot cross The two paths must cross (iii) Otherwise Undetermined Table 1: Possible cases in sequential querying of S(X (1) ) and S( ← − X (2) ) to determine whether crossing has occurred, when dealing with the case X then crossing can be determined in the same manner by swapping the superscripts.
into case (i), then we can determine that no crossing of X (1) and ← − X (2) occurs (τ (Z) = ∞). Having simulated a skeleton for X (1) and ← − X (2) and determined whether they cross, we can determine whether the candidate for Z (x 0 ,x T ,T ) is to be accepted. If they do not cross, then we repeat the procedure above until a pair X (1) and ← − X (2) are obtained that do cross. If they do cross (as illustrated in Figure 2(c)) then we can simply return a skeleton for Z as suggested by (2) (as shown in Figure 2(d)). S(Z) is a more complicated object, and some care is required in its construction as we have only identified the interval in which the crossing has occurred and not τ (Z) itself. As a consequence, we need to retain both candidate skeletons for the interval in which crossing has occurred in order to conditionally reveal Z if required in this interval (as we describe later). For clarity, we provide an expression for S(Z) which, following the procedure evaluated above, will be resolved more finely than either of the original skeletons S(X (1) ) or S( ← − X (2) ). For simplicity, we denote {0 =: ξ (0) , ξ (1) , . . . , ξ (m) , ξ (m+1) := T } as the new, finer, ordered common temporal grid associated with each skeleton following the cross-population procedure used in determining crossing. We further denote m (1) := |{ξ (0) , . . . } ∩ [0, τ (Z) ]|. Note from our notation that τ (Z) ∈ [ξ (m (1) ) , ξ (m (1) +1) ], and we retain within S(Z) the layer information for both X (1) and ← − X (2) corresponding to that interval. Then We summarise the construction of the skeleton of the proposal within our CDB methodology, which we have described above, in Algorithm 3.

Embedding the proposal within Markov chain Monte Carlo
To determine whether the proposal skeleton from Algorithm 3, S(Z), is accepted or rejected, then following the simple diffusion bridge approach outlined in Section 2 we repeatedly simulate skeletons of auxiliary diffusion paths X (3) until we can determine that one intersects with Z. This is illustrated in Figure 3. Recall that T denotes the number such paths needed, as given by (5). For each realisation of X (3) we again exploit the ε-strong simulation we introduced in Section 3. To do so we simply initialise the trajectory by simulating from ν, and then apply [32, Algorithm 4] as in Section 3 to obtain S(X (3) ). To determine whether the auxiliary path has crossed the proposal Z, we begin by cross-populating S(Z) and S(X (3) ) so that they share a common set of temporal points (and associated layers) by application of [32,Algorithm 22].
Following this cross-population it is important we do not lose track of the interval containing τ (Z) . If a new temporal point has been added in the time interval containing τ (Z) , we run Algorithm 3 in this interval to identify the subinterval to which τ (Z) belongs. Subsequently we populate S(X (3) ) with any additional temporal points in order to maintain a common resolution.
To determine whether the auxiliary path X (3) crosses the proposal path Z or not, that is whether for we have τ (aux) < ∞ or τ (aux) = ∞, we can now replicate the procedure outlined in Section 4.1. We proceed in almost the same fashion as we did for determining whether X (1) and ← − X (2) crossed: As presented in Section 3, ε-strong simulation can be used to determine whether two diffusion paths cross one another over a given interval. The only difference is should we need to determine whether the auxiliary path crosses the proposal path in the interval containing τ (Z) , then to apply Section 3 we first temporally bisect (iteratively if necessary) both S(Z) and S(X (3) ).
Having developed a procedure to determine whether X (3) and Z cross, we can now return to the simple diffusion bridge approach of Bladt & Sørensen [9], outlined in Section 2: a proposal Z is employed within a pseudo-marginal independence sampler for which 1/π T (Z) appearing in the Metropolis-Hastings ratio is estimated by T in (5). We summarise this in Algorithm 4. For presentation we resolve whether crossing has occurred in temporal order, but this is not necessary-any crossing will do-and therefore the search can often be accelerated by looking for "trivial" crossings first, before having to employ any refinements of layers.
Note that in determining whether the auxiliary path crossed the proposal, we have revealed the proposal path at additional time-points. Importantly, we never discard any additional information that has been sampled for the proposal path Z. If after sampling the auxiliary path it turns out that the two trajectories did not cross, we must keep all information about the path Z and carry it into the next experiment in which we sample the subsequent auxiliary path. It is the ε-strong simulation approach we employ that allows us to retain all pertinent information of the previously simulated auxiliary paths without a corresponding increase in computational cost at each iteration. T ← T + 1

An Ornstein-Uhlenbeck process
We consider the computational scaling of the CDB sampler (Algorithm 4) as the bridge duration increases (i.e. as T increases), contrasted with simulating the bridge directly using path-space rejection sampling: we use the ε-strong simulation approach of [7,32]. Recall, the primary motivation of the simple and confluent diffusion bridge approaches is to reduce the simulation cost from exponential to linear in T . Our first goal is to illustrate this reduction using the wellstudied example of the Ornstein-Uhlenbeck process, defined as the solution to the following For the experiments below we set (θ, µ, σ) = (1, 0, 1) and measure the execution time of the algorithm (in seconds) when applied to the sampling of 0-0 bridges (with x 0 -x T bridge denoting a bridge that start in x 0 and terminates in x T ) over a range of different values of bridge duration. Figure 4 makes it clear that the computational cost of the confluent diffusion bridges approach scales linearly with T , as opposed to the exponential scaling of ε-strong simulation. The latter is more efficient only when the bridge duration is very short. This is in line with our expectation: if the bridge duration is short then path-space rejection sampling will have a large acceptance probability, whereas the confluent diffusion bridge algorithm requires at least three diffusion paths to be simulated. Furthermore, over a short bridge duration it is less likely that X (1) and ← − X (2) cross, and that Z and X (3) cross. However, the motivation of this paper is to circumvent the more challenging simulation of diffusion bridges of large duration; in this regime the favourable performance of the confluent diffusion bridge approach is stark.

A Langevin diffusion
In this section we consider a far more complicated example inspired by the Monte Carlo Fusion approach of [15]. Monte Carlo Fusion is a methodology for unifying distributed analyses and inferences on shared parameters from multiple sources, into a single coherent inference. This is a commonly arising problem in areas such as distributed 'big data' problems and privacy settings. The underpinning theory requires the exact simulation of (multiple) Langevin diffusion bridge with known invariant distribution. It is natural to consider a confluent diffusion bridge as it does not introduce any approximation errors.
Consider the diffusion X, a solution to the following SDE: where we set: with g(·|µ, σ) denoting the probability density function of a Gaussian random variable with mean µ and standard deviation σ. The 9 constants in (12) were set to the following values: Equation (11) defines a Langevin diffusion with invariant density proportional to f , which for the particular choice (12) of f becomes a mixture of Gaussian densities. An example trajectory of a 0-0 diffusion bridge, simulated using the confluent diffusion bridges algorithm is given in Figure 5. Note that simulating diffusion bridges of long duration from (11) is a challenging problem: the process can exhibit multimodal behaviour and can yield trajectories that do not resemble paths of Brownian motion. The numerical simulations conducted had an acceptance rate above 97%, suggesting that the proposal bridges are very close to being independent samples from the target diffusion law, and again indicating the confluent diffusion bridge approach can offer substantial computational advantage over competing methodologies.

Conclusion
In this paper we proposed a novel approach to simulate bridges for ergodic diffusions, based on the combination of three algorithms: ε-strong simulation of Brownian motion, path space rejection sampling, and simple diffusion bridges. The result is a methodology which produces unbiased samples from the target diffusion bridge measure together with bounds on local minima and maxima of the sampled trajectory (as in ε-strong simulation). The resulting algorithm possesses the desirable property that the computational time scales linearly with bridge duration (just as in simple diffusion bridges), while remaining exact (just as in path-space rejection sampling). A particularly interesting feature is it is well-suited to bridging over distant, rather than near, points. That the independence sampler gets better with increasing T contrasts with existing approaches. Indeed, it is also clear that the path space rejection sampler or ε-strong simulation of [7] is still preferable for small T . A useful by-product of acquiring local maxima and minima is that it makes it possible to compute, without introducing any approximation errors, various quantities that depend on a path having crossed some pre-specified (random) barriers. Such problems are common for instance in finance, and more precisely in pricing barrier options (see [34]).
Our numerical studies in Section 5 verify empirically that we attain the linear-in-time computational scaling. As discussed in Section 5.2 this is particularly useful in practical approaches for unifying distributed analyses and inferences on shared parameters from multiple sources, into a single coherent inference. As such, one fruitful direction of this work would be to embed it within the fusion approach of [15].
Extending our methodology to a multi-dimensional diffusion bridge setting is a promising direction. Although this seems challenging as we rely on the crossing of unconditioned diffusions, the simple diffusion bridges methodology has in fact been extended to multi-dimensional ergodic diffusions by [10].
Another clear course for future work is to extend the class of models to include jump diffusion bridges. Again, there has been some work in the context of this for path-space rejection sampling: [13] consider the setting of unconditioned jump diffusions, and [31] and [21] consider conditioned jump diffusions. A confluent diffusion bridge approach would be of particular promise as it could potentially avoid a number of strong (and unnatural) assumptions that are needed in existing rejection sampling algorithms for jump diffusion bridges (in particular, assumptions on the jump intensity).
More broadly, having available diffusion simulation techniques which are linear-in-time between end points opens up new avenues in the developing area of Bayesian inference for diffusions. For instance, it is commonplace when considering inference for discretely observed diffusions to employ data augmentation schemes (or knots) [27], and so this approach could aid insight into the optimal placing of such knots. Another direction is to study whether the methodology could itself be adapted to construct a perfect simulation algorithm for diffusion bridges using coupling-from-the-past (noting that our methodology is itself underpinned by an independence sampler). detailed comments helped identify a methodological problem with an earlier version of this article. PJ, MP, and GOR were supported by The Alan Turing Institute under the EPSRC grant EP/N510129/1. GOR was additionally supported under the EPSRC grants EP/K034154/1, EP/K014463/1, and EP/R018561/1. Time X (a) Initial skeleton of (non-crossing) auxiliary path X (3) . Time X (b) Skeleton of (non-crossing) auxiliary path X (3) , at a resolution sufficient to show no crossing occurs.
Time X (c) Initial skeleton of (crossing) auxiliary path X (3) . Time X (d) Skeleton of (crossing) auxiliary path X (3) , at a resolution which shows crossing occurs at an interval which does not coincide with the intersection interval of S(Z).