On the First Crossing of Two Boundaries by an Order Statistics Risk Process

We derive a closed form expression for the probability that a non-decreasing, pure jump stochastic risk process with the order statistics (OS) property will not exit the strip between two non-decreasing, possibly discontinuous, time-dependent boundaries, within a finite time interval. The result yields new expressions for the ruin probability in the insurance and the dual risk models with dependence between the claim severities or capital gains respectively.


Introduction
The probability that a stochastic process stays between two boundaries is of significant importance in many areas among which: statistics, for constructing confidence intervals for distribution functions, (see Wald and Wolfovitz 1939;Steck 1971) or in sequential analysis; finance, for pricing double-barrier options, (see Borovkov and Novikov 2005); actuarial science, in modelling the surplus of an insurance company (see Teunen and Goovaerts 1994), and also in probability (see e.g., Potzelberger and Wang 2001;Lotov 1996;Buonocore et al. 1990) and other areas.
The double-boundary crossing problem we are concerned with in this paper is as follows.We consider a stochastic process S t , constructed by the two independent sequences of random variables 0 ≤ T 1 ≤ T 2 ≤ . . .and 0 < Y 1 < Y 2 < . . ., with the assumption that at each instant T i the process jumps at the level Y i , i = 1, 2, . ... It is assumed that the instants T 1 , T 2 , . .., model the arrival of certain events, e.g., insurance claims, or capital gains or some other events in a particular application.The jumps, 0 < Y 1 < Y 2 < . . ., could be interpreted as reflecting the current level of some process of interest, e.g., the aggregate claims to an insurance company, or the accumulation of reserves of a bank.It is further assumed that the instants, T 1 , T 2 , . . .form a point process ξ, on [0, ∞), with a cumulative intensity function Λ([0, z]) = Λ(z) = Eξ[0, z] < ∞, such that, Λ(z) > 0, ∀z ∈ (0, ∞), Λ(0) ≥ 0, where ξ[0, z] denotes the number of arrivals in [0, z].
In summary, we consider a pure-jump stochastic process S t = Y ξ[0,t] with right-continuous, non-decreasing trajectories and are interested in the probability, P (g(t) ≤ S t ≤ h(t), 0 ≤ t ≤ z) that within a fixed time interval, [0, z], the process S t stays between an upper and a lower deterministic, non-decreasing time-dependent boundaries, h(t) and g(t).The interpretation of the boundaries is also dependent on the application.For example, h(t) and g(t) could model the cumulative insurance premium income and expense outgo, respectively in the insurance and dual risk models considered in the context of ruin theory.Or these could be interpreted as correspondingly, the maximum and minimum recommended bounds for a bank's loan-to-deposit ratio so that the bank maximizes its returns while at the same time remains liquid.Formal specification of the boundaries, h(t) and g(t), is provided in Section 2.
We shall make the following assumptions regarding the (risk) process S t = Y ξ[0,t] .It will be assumed that the point process, ξ, is a process with the so called order statistics (OS) property, or simply OS point process.The latter is characterized by the property that, given j arrivals in a finite interval [0, z], z > 0, the successive arrival times coincide in distribution with the order statistics of j independent and identically distributed random variables with a cumulative distribution function, F z (t), 0 ≤ t ≤ z, F z (0) = 0 and F z (z) = 1.
As noted by Lefèvre and Picard (2011), OS processes are particularly suitable for modelling claim arrivals with dependence as they capture situations where the arrival of a claim increases the likelihood of more claim-arrivals leading to non-overlapping clustering, see e.g., the birth process considered in Bühlmann (1970).In contrast, in a death point process the arrival of a claim decreases the likelihood of arrival of more claims, as is the case in modelling the arrival of claims from a group life policy covering a closed group of individuals.It should also be mentioned that OS processes are appealing in modelling population growth where (OS) birth and death processes are a key modelling tool (see e.g., Kendall. 1949;Haccou et al. 2005).
Recently, Dimitrova et al. (2014) considered more general OS point processes by allowing Λ(t) to be discontinuous at some fixed instants, i.e, F z (t) to be discontinuous at these instants, since as shown by Crump (1975), for an OS process, F z (t) = Λ(t)/Λ(z).This extends further the flexibility of the OS risk process, S t , allowing risk events to arrive not only at random moments but also at fixed instants with non-zero probability, possibly forming clusters.This is an appealing feature, e.g., in life and non-life insurance applications, see Section 4 and Dimitrova et al. (2014) for particular examples in the special case of ruin by exceeding an upper boundary.
Our aim in the present paper is to derive, under this more general definition of an OS risk process, S t , closed-form expressions for the non-exit probability P (g(t) ≤ S t ≤ h(t), 0 ≤ t ≤ z), under reasonably general assumptions on the boundaries, g(t) and h(t), allowing them to have jump discontinuities (see Theorem 1).
Its proof is based on extending two results of Steck (1971), which give the probability that the order statistics from a sample of j uniform r.v.s all lie in a j-dimensional rectangle.We show also how Theorem 1 directly yields new ruin probability formulas in the insurance and the dual risk models under reasonably general assumptions, which have recently attracted considerable amount of research (see e.g., Ignatov andKaishev 2000, 2004;Lefèvre andPicard 2011, 2014;Dimitrova et al. (2015)).We further note that our paper (see also an earlier version Dimitrova et al. (2015)) covers and extends the non-crossing probability results for parallel boundaries of Xu (2012), Goffard and Lefèvre (2017), Goffard (2017).
The paper is organized as follows.In Section 2, we prove our main result given by Theorem 1 and its discrete version, (see Corollary 2).In Section 3 we give Corollaries 3, 4, which establish new ruin probability formulas for the two important special cases, the insurance and the dual risk models with dependence between the claim severities or capital gains respectively.In Section 4, we consider the numerical implementation of the results of Theorem 1 and Corollary 2 and illustrate these results for some special cases of OS claim arrival processes, (see 1 and 2).

A Formula for the Non-Exit Probability
We assume that the two functions h(t) and g(t), t ∈ [0, +∞), have the following properties: h(t) and g(t) are non-decreasing functions, such that h(0) ≥ 0 and g(0) ≤ 0; h(t) ≥ g(t), ∀t ∈ [0, +∞); h(t) and g(t) may have jump-discontinuities, and we assume that h(t) is right-continuous and g(t) is left-continuous.We will define the inverse functions, h −1 (y) = inf{t : h(t) ≥ y}, g −1 (y) = sup{t : g(t) ≤ y}, y ∈ [0, +∞), so that h −1 (y) is left-continuous and g −1 (y) is right-continuous.We will further consider restrictions of h(t) and g(t) on [0, z], z > 0, denoted by h z (t) and g z (t), and define the corresponding restrictions of the inverse functions as h −1 z (y) = min(z, h −1 (y)) and g −1 z (y) = min(z, g −1 (y)).In view of these definitions, we will assume that the process S t does not exit if its trajectory touches either one of the two boundaries.
As noted, here we will adopt the more general definition of an OS arrival process, recently considered by Dimitrova et al. (2014).We assume that, ξ has an arbitrary cumulative intensity function Λ(z), such that Λ(z) > 0, z > 0, Λ(0) ≥ 0, and cumulative distribution function where z > 0. If X 1 , . . ., X j are independent and identically distributed random variables with distribution function F z (t) and if, X 1,j , . . ., X j,j denote their order statistics, then the process ξ is said to have the OS property if, given j arrivals in a finite interval [0, z], (i.e., given ξ[0, z] = j),j = 1, 2, . .., the successive points, 0 ≤ T 1 ≤ T 2 ≤ . . .≤ T j of the process ξ, have a conditional joint distribution which coincides with the joint distribution of the order statistics X 1,j , . . ., X j,j .
From the definition of F z (t) it can be seen that it is right-continuous and hence, with joint cdf, F y 1 , . . ., y j .We can now state our main result.
Theorem 1.The non-exit probability P(g z (t) ≤ S t ≤ h z (t), 0 ≤ t ≤ z), 0 < z < ∞, for an OS arrival process, ξ with an arbitrary distribution function, F z (t) defined as in (1), is given by: where and det δ is the determinant of an j × j matrix with elements: , and where for j = 0 we assume that the where for 1 ≤ m, l ≤ j and ρ , and where for j = 0 we assume that the multiple integral Remark 1. From ( 3) and ( 5), it follows that δ A matrix whose elements above or below the first subdiagonal are equal to zero (i.e., all elements, a ij = 0 if j − i > 1 or if i − j > 1) are called Hessenberg matrixes.The determinants of the latter are computed using a straightforward recurrence relation established by Lemma 2. For properties of Hessenberg matrixes and their determinants we refer to, e.g., Vein and Dale (1999).Since Appell polynomials and other types of Appell structures are also expressed as Hessenberg determinants (see e.g., Vein and Dale (1999), Ignatov and Kaishev 2000, 2016, and Dimitrova et al. (2014)), here we will interchangeably refer to det δ and det ρ as Appell-Hessenberg functions.
In order to prove Theorem 1, we will need an auxiliary result, stated as Lemma 1 bellow.Let X 1 , . . ., X j be independent and identically distributed random variables with distribution function F z (t), X 1,j , . . ., X j,j denote their order statistics and let, U 1 , . . ., U j be independent uniformly on [0, 1] distributed random variables with order statistics, U 1,j , . . ., U j,j .Let us also define Lemma 1.If F z (t) is a right continuous cdf then, for α 1 , . . ., α j and β 1 , . . ., β j such that 0 ≤ α i < β i ≤ z, i = 1, . . ., j and j, a natural number, we have For a proof of Lemma 1, see Appendix A.
Let us introduce the shorter notation The following Lemma gives a recurrence expression for the probability d j , obtained by Epanechnikov (1968), which directly follows, expanding the determinant in (11), with respect to the last column, with Lemma 2. For α 1 , . . ., α j and β 1 , . . ., where d 0 ≡ 1.
We will also need the following two results, due to Steck (1971).
Theorem 2. (G.P. Steck).Let u 1 , . . ., u j and v 1 , . . ., v j be such that 0 where det δ(j) is the determinant of an j × j matrix with element on the m-th row and the l-th column defined as: Corollary 1. (G.P. Steck).Under the conditions of Theorem 2, we have where We are now ready to prove Theorem 1.
Proof of Theorem 1. Recall that, since g z (t) is restricted to the interval [0, z] and g z (0) ≤ 0, two cases need to be considered, g −1 z (0) < z and g −1 z (0) = z.In the case of g −1 z (0) < z, applying the formula of total probability, we have where the summation starts from, j = 1, since for j = 0 the conditional probability is 0 when g −1 z (0) < z.For the conditional probability on the right-hand side of (13), we have where the last equality follows from Lemma 1, with ), noting that, T 1 , . . ., T j coincide in distribution with the order statistics, X 1,j , . . ., X j,j .The asserted equalities (2) and (4) now follow substituting, P g z (t) expressed by the last equality in ( 14), back in (13) and applying Theorem 2 and Corollary 1 to express In the case of g −1 z (0) = z, for j = 0, the conditional probability is 1 and the term P(ξ[0, z] = 0) should be added to the sum on the right-hand side of ( 13).This completes the proof of the assertion of Theorem 1.
It can directly be seen that Theorem 1 is also valid for discrete Y 1 , Y 2 , . . . in which case multiple integration is replaced by appropriate multiple summation, and (2) and (4) become exact, closed form expressions more appealing for numerical implementation.This is illustrated by Corollary 2 given next, where formula ( 15) is a discrete version of (4).
Corollary 2. The non-exit probability P(g z (t) ≤ S t ≤ h z (t), 0 ≤ t ≤ z), 0 < z < ∞, for an OS arrival process, ξ, with an arbitrary distribution function, F z (t) defined as in (1), and arbitrarily distributed (possibly dependent) discrete integer-valued Y 1 , Y 2 , . . .with joint probability mass function P(Y 1 = y 1 , . . ., Y j = y j ) = p(y 1 , . . ., y j ), is given by: is the determinant of an j × j matrix with elements, ρ (j) m,l defined as in ( 5) and where for j = 0 we assume that the multiple summation in ( 15) is equal to 0 Proof of Corollary 2. The proof directly follows the lines of the proof of Theorem 1, replacing multiple integration with the corresponding multiple summation.

Application to Ruin Probability
In this section we will consider two corollaries of Theorem 1, which yield new explicit formulas for the probability of non-ruin in the insurance and dual risk models.

Ruin in the Insurance Risk Model
Ruin under the insurance risk model has been thoroughly investigated, but very few closed form expressions have been established.See e.g., Picard and Lefèvre (1997), Ignatov andKaishev 2000, 2004;Ignatov et al. (2001); Lefèvre and Loisel (2009) and Dimitrova et al. (2016) for finite-time ruin formulas under Poisson claim arrivals, involving Appell polynomials.Ruin formulas under more general arrival processes have been recently established: by Ignatov and Kaishev (2016) for the case when the claim arrival process ξ is a process with independent increments; and by Lefèvre and Picard 2011, 2014and Dimitrova et al. (2014) for the case when ξ is a process with the OS property.
In the insurance risk model, the upper boundary, h(t) models the cumulative premium income of an insurance company up to time t, and S t = Y ξ[0,t] is interpreted as its aggregate claim amount process with partial sums of the individual claims, Y 1 , Y 2 , . .., arriving at the instants, T 1 , T 2 , . ... Denote by R t = h(t) − S t the insurance company's risk process and define the instant of ruin, T as , where the lower boundary, g z (t) = 0 for every t ∈ [0, z], i.e., is such that, g −1 z (y) ≡ z, 0 ≤ y < ∞.We can now give the following.

Ruin in the Dual Risk Model
The dual risk model has been intensively studied in the ruin probability literature which is summarized in Dimitrova et al. (2015), where an instructive link between the insurance and the dual risk models is established and exploited to obtain new ruin formulas, for the dual model.
In the dual model, the lower boundary, g(t) models the cumulative expense outgo of a company up to time t, and S t = Y ξ[0,t] is interpreted as its aggregate capital gain process with partial sums of the individual capital gains, Y 1 , Y 2 , . .., arriving at the instants, T 1 , T 2 , . ... Such companies could typically, be mineral exploration companies or companies specializing in research and development, where capital gains occur at random instants of e.g., mineral finds or technological breakthroughs.
It is not difficult to see that the probability of non-ruin of the company, P (g(t) ≤ S t , 0 ≤ t ≤ z) coincides with the non-exit probability P (g z (t) ≤ S t ≤ h z (t), 0 ≤ t ≤ z), where the upper boundary, h z (t) = +∞ for every t ∈ [0, z], i.e., is such that, h −1 z (y) ≡ 0, 0 ≤ y ≤ ∞.We can now give the following.

Numerical Considerations and Examples
In this section, we give details on how our main result established by Theorem 1 can be implemented numerically and illustrate it based on two particular examples of OS point processes modelling insurance claim arrivals.
It should be noted that for the efficient evaluation of P (g z (t) ≤ S t ≤ h z (t), 0 ≤ t ≤ z), following (2), or (4), it is essential to be able to: 1. appropriately truncate the infinite summation; 2. compute the underlying (multiple) j-dimensional integrals, for possibly large values of j; 3. efficiently compute the integrand functions det δ , or equivalently det ρ (j) m,l 1≤m,l≤j .
We will address items 1-3 as follows.
1.It should be mentioned that explicit expressions for the probability P(ξ[0, z] = j) exist for the well-established OS point processes in the literature (see e.g., the four cases discussed in Goffard and Lefèvre ( 2017)), and it can directly be substituted in either (2), or (4).This allows for truncating the infinite sum with respect to j up to a maximum value say j * above which the probability ∑ ∞ j * +1 P(ξ[0, z] = j) < is negligibly small leading to negligibly small sum of the product terms in (13)), hence in (2), or (4).Similar approach has been exploited in Section 4 of Dimitrova et al. (2016) for the Poisson special case.If an expression for P(ξ[0, z] = j) is not available for an OS point process, one could view (2), or (4), as an expectation with respect to P(ξ[0, z] = j) and apply Monte Carlo simulation, generating ξ[0, z], calculating the multiple integral of dimension ξ[0, z] and then averaging over the number of simulations.2. In order to compute the j-dimensional integral in (2), or (4), one can view it as an appropriate expectation with respect to the order statistics Y 1,j , Y 2,j , . . ., Y j,j of j independent uniformly distributed on [0, h z (z)] random variables and apply a Monte Carlo approach (employing only order statistics of uniforms) to compute the integral.This approach has been applied and thoroughly investigated in Section 5 of Dimitrova et al. (2016) for computing similar multiple integrals.More precisely, following Dimitrova et al. (2016), one can rewrite (4) as × f (y 1 , . . ., y j ) × 1 V j dy j . . .dy 1 , where As is well known, 1/V j can be viewed as the joint probability density function of the order statistics, Y 1,j , Y 2,j , . . ., Y j,j , of j independent uniformly distributed on [0, h z (z)] random variables.Therefore, the multiple integrals in (2), or (4), can be interpreted as expectations with respect to the order statistics Y 1,j , Y 2,j , . . ., Y j,j and so, a Monte-Carlo simulation approach can be used to evaluate these.Thus, one may simulate, say N > 0, samples from the j-tuples Y 1,j , Y 2,j , . . ., Y j,j and utilizing the method for truncating the infinite sum so that a prescribed accuracy is achieved, as outlined in item 1 above (see also section 4 of Dimitrova et al. (2016) , we arrive at It should be noted that this Monte Carlo approach can be applied also when the Y i 's are discrete, i.e., with respect to (15), since although it is exact, its evaluation may still be prohibitively time consuming for very large values of h z (z) .See Dimitrova et al. (2016) for numerical details and comparison with the direct Monte Carlo approach.3. Finally, as noted in Remark 1, the determinants in (2) and ( 4) are of a Appell-Hessenberg type, they can be efficiently evaluated using recurrence formula (8), noting that with Next for illustrative purposes we specify formula (4), for the case of two particular OS arrival processes, the non-homogeneous linear birth process, known also as the Pólya-Lundberg process for which F z (t) is continuous, and an (OS) Poisson process with Poisson clusters at fixed instants, where F z (t) has jump discontinuities.
First we consider a double boundary non-crossing problem equivalent to a non-ruin probability model, as described in Section 3.1.We assume the arrival of insurance claims is modelled by a Pólya-Lundberg process with parameters λ > 0 and b > 0. The latter is a non-homogeneous linear birth process for which ξ[0, z] has a negative binomial distribution, i.e., As known, the Pólya-Lundberg process has cdf F z (t) = t/z for 0 ≤ t ≤ z with F z (t) = 1 for t ≥ z.It admits representation as a mixed poisson Process, i.e., ξ(t) = N(Wt), where N(t) is a homogeneous Poisson process with unit rate and W is Gamma distributed with shape parameter 1/b and scale parameter 1/λb, i.e., W ∼ Gamma(1/b, 1/λb).Assuming such an OS process for the claim arrivals, we specify the remaining parameters of the underlying OS risk model as follows.
Example 1.Consider a double boundary non-crossing problem equivalent to a non-ruin probability model as described in Section 3.1, where the arrival of claim amounts is modelled by a Pólya-Lundberg process with parameters λ = 2 and b = 1.Assume the upper boundary equivalent to the premium income function is h z (t) = t 2 + 1.5, the time horizon z = 2, the lower boundary function g z (t) ≡ 0, ∀t ≥ 0 and the partial claim sums, Y 1 , Y 2 , . . . to be the integers 1, 2, . .., respectively.Utilizing expression (4) of Theorem 1, and more precisely Corollary 2, for the non-ruin probability we have where we have used that h z (z) = 5, g −1 z (0) = z and where the number 0.568265 is computed using the recurrence formula (8) to evaluate the determinants, following (21) and using (17) substituting in it F z (h −1 z (y l )−) = min(z, y l − 1.5)/z for y l ≥ 1.5 or zero otherwise.
Let us note that this example has been considered by Goffard and Lefèvre (2017) and our result of P(T > z) = 0.568265 coincides with theirs.
Second, to demonstrate the generality and flexibility of our model, consider an OS point process, ξ, with a continuous component and a pure jump component in the underlying cdf, F z (x), i.e., ξ(0, x] = η(0, x] + ξ 1 1 (0,x] (t 1 ) + . . .+ ξ j 1 (0,x] (t j ), where η(0, x] is a Poisson process with unit rate, defined on (0, ∞] and independent of the Poisson random variables, ξ 1 , . . ., ξ j , ξ i ∈ P (µ i ), i = 1, . . .j, assumed also mutually independent.By construction, ξ is an OS process with independent increments.It could be particularly suitable for applications, especially when data comes from two (or more) independent insurance portfolios (lines of business), among which one with claim frequency data at fixed instants 0 < t 1 < t 2 < . . .< t j (e.g., annual observations) and a second one with data at policy level of the instants of claiming.In view of the Solvency II requirements, it would be instructive to be able to evaluate the probability of non-ruin in a finite time interval, (0, z], due to claims coming from all lines of business.Without loss of generality, we will illustrate this, based on the following example. Example 2. Let 0 < t 1 < t 2 ≤ z.Then, the cumulative intensity function, Λ(x) is x + µ 1 for t 1 ≤ x < t 2 , x + µ 1 + µ 2 for t 2 ≤ x ≤ z, From (A1), we have The third equality in (A2) follows since, for a right-continuous cdf F z (t) the inequality u ≤ F z (t) holds iff F −1 z (u) ≤ t (see e.g., David and Nagaraja (2003)), where one can replace u with U i,j .Consider now the probability on the left-hand side of (7).We have where in the second equality of (A3) we have used (A2) and the last equality follows since U i,j , i = 1, . . ., j are absolutely continuous.
m,l 1≤m,l≤j are lower Hessenberg matrix.