Parisian Time of Reflected Brownian Motion with Drift on Rays and Its Application in Banking

Abstract

In this paper, we study the Parisian time of a reflected Brownian motion with drift on a finite 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 first time that the duration of settlement delay exceeds a predefined 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.

1. 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_ui 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_di 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.

2. 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:={_S1,…,_Sn}, and fix a distributionP:=_{{Pi}i=1,…,n}, so that_{∑i=1n Pi}=1. We also define the functionsμ(_Si):=_μiandσ(_Si):=_σifori=1,…,n, where_μi∈Rand_σi∈^R+ are constants (see Figure 1).

Consider a planar processX(t)on the system of rays S. We represent the position ofX(t)by(|X(t)|,Θ(t)), where|X(t)|denotes the distance betweenX(t)and the origin, andΘ(t)∈{_S1,…,_Sn}indicates the current ray of the process. LetU(t):=μ(Θ(t))t+σ(Θ(t))_Wtbe the “driving process”, and|X(t)|be the Skorokhod reflection ofU(t), i.e.,

|X(t)|=U(t)+max0≤s≤t^(−U(s))+,t≥0.

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).

We expectΘ(t)to be constant during each excursion ofX(t)away from the origin and has the same distribution as P whenX(t)returns to the origin. To this end, we initializeΘ(t)withP(Θ(0)=_Si)=_Pi,i=1,…,n, and letΘ(t)remain constant whenever|X(t)|≠0. Once|X(t)|=0,Θ(t)is randomized according to P, i.e.,

PΘ(t)=_Si∣|X(^t−)|=0=_Pi,i=1,…,n,∀t>0.

This means the coefficients ofU(t)remain constant whenever|X(t)|≠0, and the Skorokhod reflection ofU(t)has the same distribution as a reflected Brownian motion with drift_μiand dispersion_σion each ray_Si.

Therefore, we summarize the behaviour ofX(t)as follows. The initial state ofX(t)is distributed asP(X(0)=(0,_Si))=_Pi,i=1,…,n. Then it behaves as a Brownian motion with drift_μiand dispersion_σion ray_Si, 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,_Si)∣|X(^t−)|=0)=_Pi,i=1,…,n.

There are some special cases ofX(t). When_μi=0and_σi=1fori=1,…,n,X(t)becomes a Walsh Brownian motion. Whenn=2,_P1=α=1−_P2,_μ1=_μ2=0and_σ1=_σ2=1,X(t)recovers the skew Brownian motion. We also obtain a Brownian motion with driftμand dispersionσby settingn=2,_P1=_P2=12,_μ1=μ,_μ2=−μand_σ1=_σ2=σ; and a reflected Brownian motion by settingn=1,_P1=1,μ=0andσ=1.

Next, we define the last zero time and excursion time length ofX(t)asg(t):=sup{s≤t∣|X(s)|=0}andU(t):=t−g(t). ThenU(t)represents the time lengthX(t)has spent in the current ray since last time at the origin. On each ray_Si, there is a threshold_di>0for the excursion time length, our target is to find the first time that the threshold is exceeded byU(t). Thus, we are interested in the Parisian timeτdefined as

_τi:=inf{t≥0∣U(t)≥_di,Θ(^t−)=_Si},fori=1,…,n,τ:=mini=1,…,n_τi.

Note thatX(t)may make an excursion with infinite time length on a ray_Siif the drift_μion this ray is positive. Since our target is to study the Parisian timeτ, we are only interested in the excursion time length up to_di, even if the total length is infinite.

We need to calculate the excursion time length ofX(t), but the problem is there is no first excursion from zero; before anyt>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_Si,|^Xϵ(t)|behaves as a reflected Brownian motion with drift_μi, dispersion_σiand 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

P^Xϵ(t)=(ϵ,_Si)∣|^Xϵ(^t−)|=0=_Pi,i∈{1,…,n}.

Hence, we describe the behaviour of^Xϵ(t)as follows. The initial state of^Xϵ(t)is distributed asP(^Xϵ(0)=(ϵ,_Si))=_Pi,i=1,…,n. Then it behaves as a Brownian motion with drift_μi, dispersion_σiand starting fromϵon ray_Si, 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

_τiϵ:=inf{t≥0∣^Uϵ(t)≥_di,^Θϵ(^t−)=_Si},fori=1,…,n,^τϵ:=mini=1,…,n_τiϵ,

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 limitlimϵ→0E(^e−βτϵ)to calculate the Laplace transform of the Parisian timeτ.

3. Laplace Transform ofτ

We present the main result of the paper in this section. For simplicity, we denote the symmetric function

Ψ(x):=2πxΦ(2x)−πx+^e−x2,x∈R,

whereΦ(.)is the cumulative distribution function of standard normal distribution, and the constant

_Ci:=_Pi22π_{σi2 di}Ψ_μidi2_σi2+_{μi σi2},

where_μi,_σi,_Piand_di are defined in Section 2. For_μi∈R,_σi∈^R+,_Pi∈(0,1]and_di>0, we deduce from the definition that_Ci>0.

Theorem 1.

LetX(t)be a reflected Brownian motion with drift on a system of rays S, where_μi∈R,_σi∈^R+,_Pi∈(0,1]and_di>0are the drift, dispersion, entering probability and excursion time threshold of ray_Si,i=1,…,n. Forβ≥0, the Laplace transform of the Parisian time τ is

Ee−βτ=∑i=1n e−βdi Ci∑i=1n Ci+∑i=1n Pi ∫0di (1−e−βv)e−μi22σi2v12πσi2 v3dv,

and the expectation of τ is

E(τ)=_{∑i=1n di Ci}+_{∑i=1n Pi ∫0di} ^{e−_μi22_σi2v}12π_σi2vdv_{∑i=1n Ci}.

Proof.

We prepare some preliminary formulas for the proof. From Section 2, we know^Xϵ(t)starts fromϵon ray_Si, and behaves as a Brownian motion with drift_μiand dispersion_σias long as it does not return to the origin. Let_gi(ϵ,t)be the density of the first hitting time at 0 of such a Brownian motion, then

_gi(ϵ,t)=ϵ2π_σi2 ^{t3e−(ϵ+_μit)22_σi2t},for_μi∈R,_σi∈^R+,ϵ>0,t>0,i=1,…,n.

We define the following functions inϵ,

Li(ϵ):=∫0di e−βs gi(ϵ,t)dtandUi(ϵ):=∫di∞ e−βdi gi(ϵ,t)dt,

and call them_Liand_Uifor convenience. In their Laplace transforms respectively,_Lirepresents the excursion time length of^Xϵ(t)on ray_Siif it is shorter than the threshold_di, and_Uirepresents the excursion time length if it is longer than_di. In the latter case, we set the excursion time length to be_dibecause we are only interested in the excursion up to the threshold.

These functions have the limits

limϵ→0_Ui(ϵ)=0andlimϵ→0_Li(ϵ)=1.

Moreover, we calculate the limits of their derivatives to be

limϵ→0ddϵ_Ui(ϵ)=^e−β_di22π_{σi2 di}Ψ_μidi2_σi2+_{μi σi2},

limϵ→0ddϵ_Li(ϵ)=−22π_{σi2 di}Ψ_{μi2 di}2_σi2+β_di−_{μi σi2}=−22π_{σi2 di}Ψ_μidi2_σi2+_{μi σi2}−_∫0di (1−^e−βy)^{e−_μi22_σi2y}12π_σi2 ^y3dy,

the last equation can be checked usingΨ(x)=1+_∫01(1−^e−x2v)12^v3/2dv, which is obtained by a direct calculation from the definition ofΨ(x).

Now we study the Parisian time^τϵ. Define the sequence of random times

_ζ0=0,_ζm+1=inf{t>_ζm∣|^Xϵ(t)|=0},form∈_N0

recursively, and the mutually exclusive events

_Cm:={_ζm≤^τϵ<_ζm+1},form∈_N0.

Then,_Cmdenotes 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)=(ϵ,_Si)}for an arbitrary but fixed i, and calculate the Laplace transformsE(^e−βτϵ _??{Cm}∣^Xϵ(0)=(ϵ,_Si))form∈_N0.

Form=0, we interpret_C0as follows. Starting fromϵon ray_Si,^Xϵ(t)spends more than_ditime before hitting the origin, hence the exceeding occurs during the first excursion. This is equivalent to the event that a Brownian motion with drift_μiand dispersion_σispends more than_ditime to travel fromϵto 0, which has probability_{∫di∞ gi}(ϵ,t)dt. Thus(^τϵ _??{C0}∣^Xϵ(0)=(ϵ,_Si))=_di, and

Ee−βτϵ ??{C0}∣Xϵ(0)=(ϵ,Si)=∫di∞ e−βdi gi(ϵ,t)dt=Ui.

Next, we consider the event_C1. In this case, the duration of the first excursion of^Xϵ(t)is shorter than_di, and the Laplace transform of the duration is_∫0di ^e−βs _gi(ϵ,t)dt. After the first excursion,^Xϵ(t)returns to the origin and jumps to(ϵ,_Sk)with probability_Pk, then exceeds the excursion time threshold_dkbefore returning to the origin. The behaviour of^Xϵ(t)during the second excursion is similar to what we described for_C0, with the index i replaced by k. Thus, we have

E^e−βτϵ _??{C1}∣^Xϵ(0)=(ϵ,_Si)=_∫0di ^e−βs _gi(ϵ,t)dt∑k=1n_PkE^e−βτϵ _??{C0}∣^Xϵ(0)=(ϵ,_Sk)=_Li∑k=1n_{Pk Uk}.

In the same way, we consider the event_C2. In this case, the duration of the first excursion of^Xϵ(t)is shorter than_di, with the Laplace transform_∫0di ^e−βs _gi(ϵ,t)dt. After the first excursion,^Xϵ(t)returns to the origin and jumps to(ϵ,_Sk)with probability_Pk. Restarting from(ϵ,_Sk),^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_C1, with the index i replaced by k. Hence,

E^e−βτϵ _??{C2}∣^Xϵ(0)=(ϵ,_Si)=_∫0di ^e−βs _gi(ϵ,t)dt∑k=1n_PkE^e−βτϵ _??{C1}∣^Xϵ(0)=(ϵ,_Sk)=_Li∑k=1n_{Pk Lk}∑j=1n_{Pj Uj}=_Li∑k=1n_{Pk Lk}∑j=1n_{Pj Uj}.

The same explanation applies to_Cmfor any positive integer m, i.e., the duration of the first excursion of^Xϵ(t)is shorter than_di, after that^Xϵ(t)restarts from(ϵ,_Sk)and exceeds the threshold exactly during the m-th excursion. Hence, we deduce that

E^e−βτϵ _??{Cm}∣^Xϵ(0)=(ϵ,_Si)=_Li∑k=1n_PkE^e−βτϵ _??{Cm−1}∣^Xϵ(0)=(ϵ,_Sk).

This implies a recursive structure between the Laplace transforms of^τϵconditioned on_Cmand_Cm−1, we solve for

E^e−βτϵ _??{Cm}∣^Xϵ(0)=(ϵ,_Si)=_Li ^{∑k=1n_{Pk Lk}m−1}∑j=1n_{Pj Uj},m=1,2,….

Since the exceeding of threshold may occur during any excursion of^Xϵ(t), we need to sum the result overm∈_N0, this gives

E^e−βτϵ∣^Xϵ(0)=(ϵ,_Si)=∑m=0∞E^e−βτϵ _??{Cm}∣^Xϵ(0)=(ϵ,_Si)=_Ui+∑m=1∞_Li ^{(∑k=1n_{Pk Lk})m−1}(∑j=1n_{Pj Uj})=_Ui+_Li(_{∑j=1n Pj Uj})1−_{∑k=1n Pk Lk},

the last equation holds because for each k,_gk(ϵ,t)is a probability density function on(0,∞), so for anyβ≥0,

0<∑k=1n_{Pk Lk}=∑k=1n_{Pk ∫0dk} ^e−βs _gk(ϵ,s)ds≤∑k=1n_{Pk ∫0dk gk}(ϵ,s)ds<∑k=1n_Pk=1.

Equation (6) boils down the Laplace transform of^τϵto the initial state of^Xϵ(t), which is distributed asP(^Xϵ(0)=(ϵ,_Si))=_Pi. Then we calculate the Laplace transform of^τϵto be

E^e−βτϵ=∑i=1n_PiE^e−βτϵ∣^Xϵ(0)=(ϵ,_Si)=∑i=1n_PiUi+_Li(_{∑j=1n Pj Uj})1−_{∑k=1n Pk Lk}=_{∑i=1n Pi Ui}(ϵ)1−_{∑k=1n Pk Lk}(ϵ).

Asϵ→0 , both numerator and denominator of the right hand side of (7) tend to 0, then we can calculate the limitlimϵ→0E(^e−βτϵ) using (4) and (5), and this givesE(^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.

Remark 1.

Whenn=2,_μ1=μ≥0,_μ2=−μ,_σ1=_σ2=1,_P1=_P2=12and_d1>0,_d2>0 , Equation (2) becomes the Laplace transform of the two-sided Parisian time of a Brownian motion with drift

Ee−βτ=e−βd1d2Ψ(μd12)+μd1 d2π2+e−βd2d1Ψ(μd22)−μd1 d2π2d2Ψ(β+μ22)d1+d1Ψ(β+μ22)d2,

this is the main result of Dassios and Wu (2010). Moreover, forn=2,_P1=_P2=12,_μ1=_μ2=0,_σ1=_σ2=1, we set_d2>0and let_d1→∞ , then Equation (2) gives the Laplace transform of the one-sided Parisian time of a standard Brownian motion

E(^e−βτ)→11+2πβ_d2exp(β_d2)Φ(2β_d2),

this was derived in Section 8.4.1 of Chesney et al. (1997).

4. 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,_Pi,_diand calculate_Cifori=1,…,n. Setλ=_{∑i=1n Ci}.

2.

Generate a multinomial random variable I whose probability function is

P(I=i)=_{Ci∑j=1n Cj}fori=1,…,n,

via the following steps:

(a)

Generate an uniform random variable_U1∼U[0,1].

(b)

SetP(I=0)=0. Fori=1,…,n, find the unique i such that

∑j=0i−1P(I=j)<_U1≤∑j=0iP(I=j),

then returnI=i.

3.

Generate a random variable^τ*via the following steps:

(a)

Generate an exponential random variableT∼exp(λ)by setting_U2∼U[0,1], then returnT=−1λln(1−_U2).

(b)

For eachi=1,…,n, generate the following subordinator:

• If_μi=0, generate a subordinator_Xiby settingα=12and

  _Xi=AlgorithmIT_Pi2π_{σi2 di}Γ(1−α)α;

• If_μi≠0, generate a subordinator_Xiby settingα=12and

  _Xi=AlgorithmIIT_Pi2π_{σi2 di}Γ(1−α)α,_{μi2 di}2_σi2.

(c)

Set^τ*=_{∑i=1n di Xi}.

4.

Outputτ=^τ*+_dI.

Proof.

For simplicity, we denote byM:=∑i=1n e−βdi Ciandλ:=∑i=1n Ci , then the Laplace transform (2) can be written as

E^e−βτ=Mλ+_{∑i=1n Pi ∫0di} (1−^e−βv)^{e−_μi22_σi2v}12π_σi2 ^v3dv,forβ≥0.

Since_Ci>0fori=1,…,n, we knowλ>0, and the denominator ofE(^e−βτ)is positive. This enables us to rewrite the Laplace transform in an integration format using the exponential function

E^e−βτ=M_∫0∞exp−tλ+∑i=1n_{Pi ∫0di} (1−^e−βv)^{e−_μi22_σi2v}12π_σi2 ^v3dvdt=Mλ_∫0∞λ^e−λtexp−∑i=1nt_Pi2π_{σi2 di∫01}(1−^e−β_diz)^{e−_{μi2 di}2_σi2z}1^z3/2dzdt.

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)=_{Ci∑j=1n Cj}fori=1,…,n,

then we can generate I using the strip method, this becomes Step 2. Note that the random variable_dI={_d1,…,_dn}has the Laplace transform

Ee−βdI=∑i=1ne−βdiCi∑j=1n Cj=Mλ.

Next, we denote by^τ*the random variable whose Laplace transform is

E^e−βτ*=_∫0∞λ^e−λtexp−∑i=1nt_Pi2π_{σi2 di∫01}(1−^e−β_diz)^{e−_{μi2 di}2_σi2z}1^z3/2dzdt

For each i, we interpret the expression

exp−t_Pi2π_{σi2 di∫01}(1−^e−β_diz)^{e−_{μi2 di}2_σi2z}1^z3/2dz

as the Laplace transform of the random variable_{di Xi}(t_Pi2π_{σi2 di}), where_Xi(t_Pi2π_{σi2 di})is a subordinator with truncated Lévy measure

^{e−_{μi2 di}2_σi2z}1^z3/2_??{0<z<1}dz

at timet_Pi2π_{σi2 di} . Comparing (10) with (A1), we know_Xi(.) can be generated via Algorithms 4.3 and 4.4 in Appendix A.

Moreover, (9) implies that^τ*=law_{∑i=1n di Xi}(T_Pi2π_{σi2 di}), whereT∼exp(λ)is an exponential random variable. Hence, we generate T in Step 3(a), sample from_Xi(T_Pi2π_{σi2 di})in Step 3(b) and calculate^τ*via Step 3(c).

Finally, sinceE(^e−βτ)=E(^e−β_dI)E(^e−βτ*), we have the representationτ=law_dI+^τ*, where_dIand^τ*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 setn=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;_P1=0.1,_P2=0.2,_P3=0.1,_P4=0.2,_P5=0.2,_P6=0.1,_P7=0.1;_d1=1,_d2=3,_d3=2.5,_d4=1.5,_d5=1.5,_d6=0.5,_d7=2.5.

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 be3.0534. Then we consider the following two standard measures for the associated error of the algorithm,

• difference = sample average − true expectation

• standard error=samplestandarddeviationnumberofsamples

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 transformE(^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 [email protected] 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.

5. 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_Pito make another payment to bank i,i=1,…,n. Let a reflected Brownian motion with drift_μiand dispersion_σibe the net balance between bank A and bank i. To avoid the excessive exposure to liquidity risk, a time limit_dihas 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 processX(t), and view the first exceeding time as the Parisian time ofX(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.

Sample Size	Sample Average	Difference	Standard Error
1000	3.0666	0.0132	0.0616
4000	3.0302	−0.0232	0.0304
16,000	3.0470	−0.0065	0.0155
64,000	3.0509	−0.0025	0.0077
256,000	3.0520	−0.0014	0.0039
1,024,000	3.0538	0.0004	0.0019

Sample Size	CPU Time (in seconds)
1000	0.201831
4000	0.725738
16,000	3.080721
64,000	12.214876
256,000	52.715700
1,024,000	201.460605

Author Contributions

Methodology, A.D. and J.Z.; software, J.Z.; formal analysis, A.D. and J.Z.; writing-original draft preparation, J.Z.; writing-review and editing, A.D.; supervision, A.D.; All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The authors are grateful to the editor for handling this manuscript and both reviewers for giving us very helpful and inspiring comments, especially during the difficult time of pandemic.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Exact Simulation of Truncated Subordinator

In this appendix, we attach Algorithms 4.3 and 4.4 of Dassios et al. (2020). These algorithms exactly generate samples from the truncated subordinatorZ(t)with Laplace transform

E^e-vZ(t)=exp-αtΓ(1-α)_∫01(1-^e-vz)^{e-ηz zα+1}dz.

We first present two ancillary algorithms, namely Algorithms 4.1 and 4.2 of Dassios et al. (2020).

Lemma A1

(Algorithm 4.1 of Dassios et al. 2020). Exact simulation of(T,W).

1.

Setξ=Γ^(1-α)-1;_A0=(1-α)^αα1-α.

2.

minimiseC(λ)=_A0 ^{eξ1α λ1-1αα(1-α)1α-1 (_A0-λ)α-2}.

3.

record critical value^λ*; setC=C(^λ*).

4.

repeat {

5.

sampleU∼U[0,π];_U1∼U[0,1],

6.

setY=1-_U111-α;_AU=^{[sinα(αU)sin1-α((1-α)U)/sin(U)]11-α},

7.

sampleR∼Γ(2-α,_Au-λ);V∼U[0,1].

8.

if(V≤_AU ^{eξR1-α Yα e-λ*R (_AU-λ*)α-2 Yα-1}(1-^(1-Y)α)/C), break.

9.

}

10.

sample_U2∼U[0,1],

11.

setT=^{R1-α Yα};W=Y-1+^{[(1-Y)-α-_U2((1-Y)-α-1)]-1α}.

12.

return(T,W).

Lemma A2

(Algorithm 4.2 of Dassios et al. 2020). Exact simulation of{Z(t)|T>t}.

1.

sample_U1∼U[0,π]; set_AU1 =^{[sinα(α_U1)sin1-α((1-α)_U1)/sin(_U1)]11-α}.

2.

repeat {

3.

sample_U2∼U[0,1]; setZ=^{-log(_U2)_AU1 t11-α-1-αα}.

4.

if(Z<1), break.

5.

}

6.

return Z.

Next we provide the Algorithm 4.3 and 4.4 of Dassios et al. (2020).

Theorem A1

(Algorithm 4.3 of Dassios et al. 2020). Exact simulation of the subordinatorZ(t)whenη=0. The input is t.

1.

setZ=0;S=0.

2.

repeat {

3.

sample(T,W)via Algorithm 4.1; setS=S+T,Z=Z+1+W.

4.

if(S>t), break.

5.

}

6.

set_ZS-T=Z-1-W; sample_Zt-(S-T)via Algorithm 4.2.

7.

return_ZS-T+_Zt-(S-T).

Theorem A2

(Algorithm 4.4 of Dassios et al. 2020). Exact simulation of the subordinatorZ(t)whenη>0. The inputs are(t,η).

1.

repeat {

2.

sample_Ztvia Algorithm 4.3;V∼U[0,1].

3.

if(V≤exp(-η_Zt)), break.

4.

}

5.

return_Zt.

Proof.

For the proof as well as the motivation of the algorithms above, see Dassios et al. (2020). □