Parisian Time of Reﬂected Brownian Motion with Drift on Rays and Its Application in Banking

: In this paper, we study the Parisian time of a reﬂected Brownian motion with drift on a ﬁnite collection of rays. We derive the Laplace transform of the Parisian time using a recursive method, and provide an exact simulation algorithm to sample from the distribution of the Parisian time. The paper was motivated by the settlement delay in the real-time gross settlement (RTGS) system. Both the central bank and the participating banks in the system are concerned about the liquidity risk, and are interested in the ﬁrst time that the duration of settlement delay exceeds a predeﬁned limit. We reduce this problem to the calculation of the Parisian time. The Parisian time is also crucial in the pricing of Parisian type options; to this end, we will compare our results to the existing literature.


Introduction
Suppose we have a system of rays emanating from a common origin and a particle moving on this system. On each ray, the particle behaves as a reflected Brownian motion with drift; and once at the origin, it instantaneously chooses a ray for its next excursion randomly according to a given distribution. We are interested in the time length the particle spends on each ray, and the first time that the excursion time length on a ray exceeds a predefined threshold. We call this first exceeding time of threshold a Parisian time, as it generalizes the concept of Parisian time of standard Brownian motion in literature.
The study of excursion time length of Brownian motion goes back to Chung (1976). Other aspects of Brownian excursion have also been considered. Durrett et al. (1977) developed the relationships between the Brownian excursions, meanders and bridges using the limit processes of conditional functionals of Brownian motion. Imhof (1984) derived the joint density concerning the maximum of Brownian motion and 3-dimensional Bessel process. Kennedy (1976) derived the distribution of the maximum of excursion via the limiting processes and relates it to the standard Brownian motion. Getoor and Sharpe (1979) obtained a limiting result on the distribution of additive functionals over Brownian excursions. A literature review can be found in Zhang (2014).
More recently, Chesney et al. (1997) studied the Parisian time of Brownian motion, and used the result to price the Parisian type options. They are path-dependent options whose payoff depends not only on the final value of the underlying asset, but also on the path trajectory of the underlying above or below a predetermined barrier for a length of time. The two-sided Parisian option was considered in Dassios and Wu (2010), its pricing depends on the Parisian time of a drifted Brownian motion with a two-sided excursion time threshold. It turns out that the Parisian times derived in Chesney et al. (1997) and Dassios and Wu (2010) can be viewed as the special cases of our result. Moreover, the results in the current paper can be used to price more complicated Parisian type options. For more details about Parisian options, see Schröder (2003), Anderluh and van der Weide (2009) and Labart and Lelong (2009). This paper is motivated by the real-time gross settlement system (RTGS, and known as CHAPS in the UK, see McDonough 1997;Padoa-Schioppa 2005). The participating banks in the RTGS system are concerned about liquidity risk and wish to prevent the excessive liquidity exposure between two banks. There is evidence suggesting that in CHAPS, banks usually set bilateral or multilateral limits on the exposed position with others (see Becher et al. 2008), this mechanism was studied by Che and Dassios (2013) using a Markov model. For a single bank, namely bank A, let a reflected Brownian motion be the net balance between bank A and bank i, and let u i be the bilateral limit set up by bank A for bank i, Che and Dassios (2013) calculated the probability that the limit is exceeded in a finite time.
We consider another source of liquidity risk, the time-lag between the execution of the transaction and its final completion. As it is explained in McDonough (1997) and Padoa-Schioppa (2005), if a counterparty does not settle an obligation for the full value when due but at some unspecified time thereafter, the expected liquidity position of the payee could be affected. The settlement delay may force the payee to cover its cash-flow shortage by funding at short notice from other sources, which may result in a financial loss due to higher financing costs or to damage to its reputation. In more extreme cases, it may be unable to cover its cash-flow shortage at any price, in which case it may be unable to meet its obligation to others. As the settlement delay is the major source of liquidity risk in the RTGS system, both the central bank and the participating banks are interested in the length of the delay. Previous research in Che and Dassios (2013) has shown that the Markov-type models are adequate for CHAPS, we will extend this model here to study the settlement delay. For bank A and bank i in CHAPS, we view the net balance between them as a reflected Brownian motion with drift. Assume that bank A has set a time limit d i on the duration of settlement delay for bank i, and they are interested in the first time that the limit is exceeded. In practice, an individual bank could set multiple limits or even remove the limit on different types of counterparties. We reduce this problem to the calculation of the Parisian time of a reflected Brownian motion with drift on rays. For more details about the CHAPS, see Che (2011) and Soramäki et al. (2007).
We construct the reflected Brownian motion with drift on rays in Section 2, then calculate the Laplace transform of the Parisian time in Section 3. An exact simulation algorithm to sample from the distribution of the Parisian time is provided in Section 4. We discuss the application of these results in Section 5.

Construction of the Underlying Process and the Parisian Time
In this section, we construct the reflected Brownian motion with drift on a finite collection of rays, and define the Parisian time we are interested in. Let n be a finite positive integer, we denote by S a system containing n rays emanating from the common origin, i.e., S := {S 1 , . . . , S n }, and fix a distribution P := {P i } i=1,...,n , so that ∑ n i=1 P i = 1. We also define the functions µ(S i ) := µ i and σ(S i ) := σ i for i = 1, . . . , n, where µ i ∈ R and σ i ∈ R + are constants (see Figure 1).
Consider a planar process X(t) on the system of rays S. We represent the position of X(t) by (|X(t)|, Θ(t)), where |X(t)| denotes the distance between X(t) and the origin, and Θ(t) ∈ {S 1 , ..., S n } indicates the current ray of the process. Let U(t) := µ(Θ(t))t + σ(Θ(t))W t be the "driving process", and |X(t)| be the Skorokhod reflection of U(t), i.e., Then, |X(t)| has the same distribution as a reflected Brownian motion with drift µ(Θ(t)) and dispersion σ(Θ(t)). A proof of this can be seen in Jeanblanc et al. (2009) Section 4.1, Peskir (2006 and Graversen and Shiryaev (2000). ray S 2 , excursion time threshold d 2 µ(S 2 ) = µ 2 , σ(S 2 ) = σ 2 ray S 1 , excursion time threshold d 1 µ(S 1 ) = µ 1 , σ(S 1 ) = σ 1 ray S 5 , excursion time threshold We expect Θ(t) to be constant during each excursion of X(t) away from the origin and has the same distribution as P when X(t) returns to the origin. To this end, we initialize Θ(t) with P(Θ(0) = S i ) = P i , i = 1, . . . , n, and let Θ(t) remain constant whenever |X(t)| = 0. Once |X(t)| = 0, Θ(t) is randomized according to P, i.e., This means the coefficients of U(t) remain constant whenever |X(t)| = 0, and the Skorokhod reflection of U(t) has the same distribution as a reflected Brownian motion with drift µ i and dispersion σ i on each ray S i . Therefore, we summarize the behaviour of X(t) as follows. The initial state of X(t) is distributed as P(X(0) = (0, S i )) = P i , i = 1, . . . , n. Then it behaves as a Brownian motion with drift µ i and dispersion σ i on ray S i , as long as it does not return to the origin. Once at the origin, it instantaneously chooses a new ray according to P, independently of the past behaviour; that is, P(X(t) = (0, S i ) | |X(t − )| = 0) = P i , i = 1, . . . , n.
Next, we define the last zero time and excursion time length of X(t) as g(t) := sup{s ≤ t | |X(s)| = 0} and U(t) := t − g(t). Then U(t) represents the time length X(t) has spent in the current ray since last time at the origin. On each ray S i , there is a threshold d i > 0 for the excursion time length, our target is to find the first time that the threshold is exceeded by U(t). Thus, we are interested in the Parisian time τ defined as Note that X(t) may make an excursion with infinite time length on a ray S i if the drift µ i on this ray is positive. Since our target is to study the Parisian time τ, we are only interested in the excursion time length up to d i , even if the total length is infinite.
We need to calculate the excursion time length of X(t), but the problem is there is no first excursion from zero; before any t > 0, the process has made an infinite number of small excursions away from the origin. To approximate the dynamic of a Brownian motion, Dassios and Wu (2010) introduced the "perturbed Brownian motion", we will extend this idea here.
For every > 0, we define a perturbed process X (t) = (|X (t)|, Θ (t)) on the system of rays S. On each ray S i , |X (t)| behaves as a reflected Brownian motion with drift µ i , dispersion σ i and starting from , as long as it does not return to the origin. Once at the origin, X (t) not only chooses a new ray according to P, but also jumps to on the new ray. In other words, X (t) has a perturbation of size at the origin which can be described as Hence, we describe the behaviour of X (t) as follows. The initial state of X (t) is distributed as P(X (0) = ( , S i )) = P i , i = 1, . . . , n. Then it behaves as a Brownian motion with drift µ i , dispersion σ i and starting from on ray S i , as long as it does not return to the origin. Once at the origin, it instantaneously chooses a new ray according to P and jumps to on the new ray.
We define the Parisian time of X (t) similarly as before. Let g (t) := sup{s ≤ t | |X (s)| = 0} and U (t) := t − g (t). We are interested in the Parisian time τ defined as As → 0, the perturbation at origin vanishes, and X (t) → X(t) in a pathwise sense, then τ → τ in distribution. Hence, we will first derive the Laplace transform of τ , then take the limit lim →0 E(e −βτ ) to calculate the Laplace transform of the Parisian time τ.

Laplace Transform of τ
We present the main result of the paper in this section. For simplicity, we denote the symmetric function where Φ(.) is the cumulative distribution function of standard normal distribution, and the constant where µ i , σ i , P i and d i are defined in Section 2. For µ i ∈ R, σ i ∈ R + , P i ∈ (0, 1] and d i > 0, we deduce from the definition that C i > 0.
Theorem 1. Let X(t) be a reflected Brownian motion with drift on a system of rays S, where µ i ∈ R, σ i ∈ R + , P i ∈ (0, 1] and d i > 0 are the drift, dispersion, entering probability and excursion time threshold of ray S i , i = 1, . . . , n. For β ≥ 0, the Laplace transform of the Parisian time τ is and the expectation of τ is Proof. We prepare some preliminary formulas for the proof. From Section 2, we know X (t) starts from on ray S i , and behaves as a Brownian motion with drift µ i and dispersion σ i as long as it does not return to the origin. Let g i ( , t) be the density of the first hitting time at 0 of such a Brownian motion, then We define the following functions in , and call them L i and U i for convenience. In their Laplace transforms respectively, L i represents the excursion time length of X (t) on ray S i if it is shorter than the threshold d i , and U i represents the excursion time length if it is longer than d i . In the latter case, we set the excursion time length to be d Moreover, we calculate the limits of their derivatives to be the last equation can be checked using Ψ(x) = 1 + 1 0 (1 − e −x 2 v ) 1 2v 3/2 dv, which is obtained by a direct calculation from the definition of Ψ(x). Now we study the Parisian time τ . Define the sequence of random times recursively, and the mutually exclusive events Then, C m denotes the event that the exceeding of threshold occurs during the (m + 1)-th excursion of X (t) away from the origin. Next, we set {X (0) = ( , S i )} for an arbitrary but fixed i, and calculate the Laplace transforms E(e −βτ 1 {C m } | X (0) = ( , S i )) for m ∈ N 0 .
For m = 0, we interpret C 0 as follows. Starting from on ray S i , X (t) spends more than d i time before hitting the origin, hence the exceeding occurs during the first excursion. This is equivalent to the event that a Brownian motion with drift µ i and dispersion σ i spends more than d i time to travel from to 0, which has probability Next, we consider the event C 1 . In this case, the duration of the first excursion of X (t) is shorter than d i , and the Laplace transform of the duration is d i 0 e −βs g i ( , t)dt. After the first excursion, X (t) returns to the origin and jumps to ( , S k ) with probability P k , then exceeds the excursion time threshold d k before returning to the origin. The behaviour of X (t) during the second excursion is similar to what we described for C 0 , with the index i replaced by k. Thus, we have In the same way, we consider the event C 2 . In this case, the duration of the first excursion of X (t) is shorter than d i , with the Laplace transform d i 0 e −βs g i ( , t)dt. After the first excursion, X (t) returns to the origin and jumps to ( , S k ) with probability P k . Restarting from ( , S k ), X (t) will exceed the excursion time threshold exactly during the second excursion (hence the third in total). The behaviour of X (t) during the second and third excursions is similar to what we described for C 1 , with the index i replaced by k. Hence, The same explanation applies to C m for any positive integer m, i.e., the duration of the first excursion of X (t) is shorter than d i , after that X (t) restarts from ( , S k ) and exceeds the threshold exactly during the m-th excursion. Hence, we deduce that This implies a recursive structure between the Laplace transforms of τ conditioned on C m and C m−1 , we solve for Since the exceeding of threshold may occur during any excursion of X (t), we need to sum the result over m ∈ N 0 , this gives the last equation holds because for each k, g k ( , t) is a probability density function on (0, ∞), so for any Equation (6) boils down the Laplace transform of τ to the initial state of X (t), which is distributed as P(X (0) = ( , S i )) = P i . Then we calculate the Laplace transform of τ to be As → 0, both numerator and denominator of the right hand side of (7) tend to 0, then we can calculate the limit lim →0 E(e −βτ ) using (4) and (5), and this gives E(e −βτ ). The expectation of τ is obtained by applying the moment generating function.
As in Section 2, X(t) can be reduced to a Brownian motion with drift or a standard Brownian motion by choosing the parameters accordingly, then we can compare Theorem 1 with the results in the existing literature.

Exact Simulation Algorithm of the Parisian Time
In this section, we provide an exact simulation algorithm to sample from the distribution of the Parisian time τ. Our algorithm is based on the exact simulation schemes of the truncated Lévy subordinator developed in Dassios et al. (2020). We refer to Algorithms 4.3 and 4.4 of Dassios et al. (2020) as AlgorithmI(.) and AlgorithmII(. , .),their full steps are attached in Appendix A.
Theorem 2. Exact simulation algorithm of the Parisian time τ.
1. Initialize µ i , σ i , P i , d i and calculate C i for i = 1, . . . , n. Set λ = ∑ n i=1 C i . 2. Generate a multinomial random variable I whose probability function is P(I = i) = C i ∑ n j=1 C j for i = 1, ..., n, via the following steps: (a) Generate an uniform random variable U 1 ∼ U[0, 1].
Since C i > 0 for i = 1, . . . , n, we know λ > 0, and the denominator of E(e −βτ ) is positive. This enables us to rewrite the Laplace transform in an integration format using the exponential function Equation (8) can be understood as a product of the Laplace transforms of two independent random variables, hence we can generate them separately, and view the Parisian time τ as their summation.
Denote by I a multinomial random variable with the probability function P(I = i) = C i ∑ n j=1 C j for i = 1, ..., n, then we can generate I using the strip method, this becomes Step 2. Note that the random variable d I = {d 1 , . . . , d n } has the Laplace transform Next, we denote by τ * the random variable whose Laplace transform is For each i, we interpret the expression as the Laplace transform of the random variable Comparing (10) with (A1), we know X i (.) can be generated via Algorithms 4.3 and
Finally, since E(e −βτ ) = E(e −βd I )E(e −βτ * ), we have the representation τ law = d I + τ * , where d I and τ * are independent, then τ can be generated via Step 4.
Next, we illustrate the accuracy and performance of the exact simulation algorithm with a numerical example. We set n = 7, and µ 1 = 0, µ 2 = 0.5, µ 3 = −0.3, µ 4 = 0, µ 5 = 0.2, µ 6 = 0, µ 7 = −0.1; σ 1 = 1.5, σ 2 = 2, σ 3 = 1.3, σ 4 = 1, σ 5 = 2, σ 6 = 1, σ 7 = 1; P 1 = 0.1, P 2 = 0.2, P 3 = 0.1, P 4 = 0.2, P 5 = 0.2, P 6 = 0.1, P 7 = 0.1; Using the exact simulation algorithm, we generate samples from the Parisian time and calculate their average. On the other hand, we use Equation (3) to calculate the true expectation of τ to be 3.0534. Then we consider the following two standard measures for the associated error of the algorithm, 1. difference = sample average − true expectation 2. standard error = sample standard deviation √ number of samples Table 1 reports the results, we see that the algorithm can achieve a high accuracy, and one has to generate more samples to decrease the standard error. In addition, we estimate the distribution function of the Parisian time. Using the exact simulation algorithm and the smoothing techniques (see Bowman and Azzalini 1997), we get the estimated curve for the distribution function. On the other hand, we apply the Gaver-Stehfest method (see Cohen 2007) to invert the Laplace transform E(e −βτ ) β numerically and obtain the inverted curve for the distribution function. These curves are provided in Figure 2, they show that the exact simulation algorithm provides a good approximation for the distribution of the Parisian time.
We also illustrate the performance of the algorithm by recording the CPU time needed to generate these samples from the Parisian time. The experiment is implemented on an Intel Core i5-5200U CPU@2.20GHz processor, 8.00GB RAM, Windows 10, 64-bit Operating System and performed in Matlab R2019b. No parallel computing is used. Table 2 reports the results.

Discussion
We can apply this model to study the settlement delay in CHAPS. For an individual bank A, we assume that there are n counterparties in the system, namely bank 1, bank 2, . . . , bank n. We also assume that bank A uses an internal queue to manage its outgoing payments, and once the current payment is settled, it has probability P i to make another payment to bank i, i = 1, . . . , n. Let a reflected Brownian motion with drift µ i and dispersion σ i be the net balance between bank A and bank i.
To avoid the excessive exposure to liquidity risk, a time limit d i has been set for the duration of settlement delay between bank A and bank i. Both the central bank and the participating banks are interested in the first time that the limit is exceeded.
We model the net balance between bank A and the counterparties by the planar process X(t), and view the first exceeding time as the Parisian time of X(t). Using the results in the current paper, we can sample from this first exceeding time and estimate its distribution function numerically. We remark that this approach can be adopted by both the policymaker in the central bank and the credit control departments of the participating banks to lay down decisive actions. For example, the central bank may use time-dependent transaction fees to provide incentives to earlier settlements. Alternatively, the participating banks may also learn to coordinate their payments over time, creating non-binding behavioural conventions or implicit contracts.
In particular, an empirical method has been developed in Denbee and Zimmerman (2012) to detect the apparent 'free-riding' in the RTGS system, referring to the behaviour that the banks wait for incoming payments to fund subsequent outgoing payments and not supply an amount of liquidity to the system commensurate with the share of payments they are responsible for. Suppose the banks are required to hold buffers of liquid assets in order that they can make payments in a stress scenario, and the buffers are continuously calculated based on past activity. Banks may have an incentive to delay their payments so that the regulatory buffer will be reduced at subsequent recalibrations. The method in Denbee and Zimmerman (2012) could help to detect this behaviour and calibrate buffers independent of strategic actions. The study in the current paper provides another point of view towards this method. We can estimate the distribution of the settlement delay and take this into consideration when calculating the buffers.
It is also possible to extend the model in the current paper to the settlement systems other than CHAPS. For example, the structure of settlement delay in Interbank Electronic Payment System (SPEI operated by Banco de México) has been specified in Alexandrova-Kabadjova and Solis (2012) with real transactions data from 7 April to 7 May 2010. We may assume that the Markov model is adequate for SPEI, and use these data to calibrate the parameters of the model. Moreover, the observations in Alexandrova-Kabadjova and Solis (2012) suggest that low value payments do not increase the settlement delay in the system. This is reasonable under the assumption that the net balance between two banks follows a reflected Brownian motion with drift, because the process will make an infinite number of small excursions at the origin. This paper has focused on the model with one central bank (or agent) and several domestic participants, which is classified as a 'within border payment system' (see Bech et al. 2020). For a cross-border payment system, however, we need to consider a model containing two or more central banks, each with their own domestic participants. Assume that the system offers payment versus payment (PvP, see Bech et al. 2020) services, then the settlement delay may originate in any local system, and the first exceeding time of settlement delay of the whole system can be viewed as the joint distribution of the Parisian times of the local systems. With the technique developed in this paper, we are able to simulate the marginal distributions of the local exceeding time, but not the joint distribution. This is a topic for future research, and the result would be beneficial on a global scale.
In addition, our Brownian-type model reflects the random fluctuations of payments and delays, but the external events that can influence these are not taken into account. For example, the operational risks related to computer and telecommunication system breakdown may increase the settlement delay, see Rochet and Tirole (1996) for the impact of computer problem of the Bank of New York in 1985 and the San Francisco earthquake in 1989. More recently, many reports have suggested the impact of global pandemic in 2020 on the settlement systems. These might be interesting for a further study. Theorem A2 (Algorithm 4.4 of Dassios et al. 2020). Exact simulation of the subordinator Z(t) when η > 0. The inputs are (t, η).