Finite-horizon general insolvency risk measures in a regime-switching Sparre Andersen model

Insolvency risk measures play important role in the theory and practice of risk management. In this paper, we provide a numerical procedure to compute vectors of their exact values and prove for them new upper and/or lower bounds which are shown to be attainable. More precisely, we investigate a general insolvency risk measure for a regime-switching Sparre Andersen model in which the distributions of claims and/or wait times are driven by a Markov chain. The measure is defined as an arbitrary increasing function of the conditional expected harm of the deficit at ruin, given the initial state of the Markov chain. A vector-valued operator L, generated by the regime-switching process, is introduced and investigated. We show a close connection between the iterations of L and the risk measure in a finite horizon. The approach assumed in the paper enables to treat in a unified way several discrete and continuous time risk models as well as a variety of important vector-valued insolvency risk measures.


Introduction
In practice, insurers face the claims which can hardly be modeled as i.i.d. random variables. Natural catastrophes may rapidly cause a significant increase in large insurance claims by affecting the distribution of claims in several succeeding time periods. From a theoretical point of view, large insurance claims could be explained by means of the so-called heavy tail assumption. However, this would simultaneously make more difficult modelling the claims in those periods when no extreme events occur. Moreover, amount of the claims may be bounded from above by the policy's general terms and conditions, removing completely the claim distribution tails from consideration. Following Reinhard (1984) and Asmussen (1989), one can solve the fat tail assumption dilemma using Markov-modulated processes.
A crucial step forward was to let both the distributions of claims and wait times be driven by a Markov chain. Applying this approach to the Sparre Andersen model leads to a sophisticated regime-switching framework which can find some important applications from the Solvency II perspective (see Section 4 of Gajek and Rudź (2018b) for a simulation study including some details on how to determine the Solvency Capital Requirement).
Numerous regime-switching models have been proposed and investigated recently, see, e.g., Cheng and Zhan (2020), Jin et al. (2020), Shen et al. (2019), Cao et al. (2018), Liu and Privault (2018), Ignatieva et al. (2018), Xu et al. (2017), Wang et al. (2016), Landriault et al. (2015), Chen and Delong (2015), Chen et al. (2014), Guillou et al. (2013), Asmussen and Albrecher (2010), Kim and Kim (2007), Lu (2006), Bäuerle (1996), Asmussen (1989) and Reinhard (1984). Some useful suggestions and justifications for the application of Markov switching models in insurance mathematics can also be found in Rabehasaina (2009), De ¸bicka (2013) and D'Amico (2014. The latter article substantially extends the results of Kim and Kim (2007). Other references can be found, e.g., in Silvestrov (2014Silvestrov ( , 2015. For a discussion of recent developments of Markov switching models and the operator approach in ruin theory, we also refer the reader to Gajek and Rudź (2018a). In particular, a research methodology based on Banach Contraction Principle is used there to approximate a vector of ultimate ruin probabilities in a regime-switching Sparre Andersen model. A generalized Gerber's upper bound for finite-horizon ruin probabilities in this model can be found in Gajek and Rudź (2017). Lower and upper estimates for the vector of finite-horizon as well as ultimate deficit distributions at ruin has been proven in Gajek and Rudź (2018c). The generality of the regime-switching Sparre Andersen model allows an immediate application of the above mentioned methods in several continuous and discrete time risk models listed in Gajek and Rudź (2018a).
In this paper, we investigate general insolvency risk measures in the following Sparre Andersen risk model with a switch. Assume first that each stochastic object considered in the paper is defined on a common probability space (Ω, F , P). Let a positive random variable X k denote the amount of the kth claim; a positive T 1 -the moment when the first claim appears; a positive T k with k ∈ N\{1} -the time between the (k − 1)th claim and the kth one; a positive C k -the insurance premium rate per time T k ; and a non-negative real uthe insurer's surplus at 0. Let C k = c(I k−1 ), where c is a known positive function defined on a state space S = {1, 2, ..., s}, s ∈ N, of a Markov chain {I k } k=0,1,2,... (see Section 2 for details). Let us define the insolvency risk process {U(n, ·)} n∈N as follows: The moment of ruin, considered as a function of the initial surplus u, is given by τ (u) = inf{n ∈ N : U(n, u) < 0} under the convention that inf ∅ means ∞. Let us recall that Ψ i n (u) = P(τ (u) n|I 0 = i) is the conditional probability of ruin at or before the nth claim, given the initial state i ∈ S of {I k } k=0,1,2,... . The vector (Ψ 1 n , ..., Ψ s n ) will be denoted by Ψ n .
One-dimensional version of Ψ n is a traditional measure of the insolvency risk in a finite horizon. It has a disadvantage that it does not take into account the severity of ruin. On the other hand, if we use the conditional expectation of the deficit at ruin to measure the severity of ruin, given the initial state i, the probability of ruin is not reflected properly. For these reasons, one can combine both of the above quantities in the following way. Let us define the finite-horizon unconditional deficit, say D n , considered as a function of the initial surplus u, as follows where n ∈ N. One can use the conditional expectation of D n , given i ∈ S and u 0, as a measure of the insolvency risk. Then when they are used separately. Therefore, we will consider the following insolvency risk measure: where χ : [0, ∞) → [0, ∞) is a harm function and Δ : [0, ∞) → [0, ∞) is a strictly increasing function. From a practical point of view, the greater the deficit is, the more it harms, so we will assume that the function χ is non-decreasing. We will use the convention that δ i 0 (u) = 0. Set δ n (u) = (δ 1 n (u), ..., δ s n (u)), where n = 0, 1, ... and u 0. Clearly, δ n provides a variety of important vector-valued insolvency risk measures. In particular, when Δ(x) = x, we get a counterpart of the so-called spectral risk measure which is well-known in mathematical finance. A relationship between spectral and distorted risk measures can be found, e.g., in Gzyl and Mayoral (2006). In an insurance context a similar approach has been introduced and popularized by Gerber and Shiu (see Gerber and Shiu (1998) for details). Our approach, with an arbitrary strictly increasing function Δ, enables us to consider a wide class of insolvency risk measures that are listed below. Let us denote E i X = E[X|I 0 = i], P i (B) = P(B|I 0 = i), X + = max{X, 0} and 1 I A (x) = 1 if x ∈ A and 0 otherwise. Special cases of the risk measure (2) are: , when y → 0 + in (vii). All the above risk measures are commonly used, but their exact values in the Sparre Andersen risk model with a switch seem to be unknown. For this reason, we will provide an iterative method (with a simulation study given in Section 6) to establish their exact values. Moreover, several methods to find their upper and lower approximations will also be given. The main results of the paper are proved in Sections 3 and 4 and summarized in Section 5.
Other methods to compute explicit expressions for various ruin-related quantities can be found, e.g., in Woo and Liu (2018) or Lorek (2017). Markov additive processes are discussed, e.g, in Asmussen (2003) or Feng and Shimizu (2014).

A vector-valued risk operator
Let us denote by N the set of all positive integers. Set N 0 = N ∪ {0}, N 1 = N\{1}, R + = (0, ∞) and R 0 + = [0, ∞). Recall that X k denotes the amount of the kth claim; T 1 , the moment when the first claim appears; T k , the time between the (k − 1)th claim and the kth one. Let us denote by A n the moment when the nth claim appears. It is easy to see that A n = T 1 + ... + T n , n ∈ N 0 , with A 0 = 0. Let a random variable C k denote the insurance premium rate during the time interval [A k−1 , A k ). Throughout the paper, we will assume that the random variables C k , T k and X k are positive (a.s.) for every k ∈ N and their distributions have no singular parts.
Let {I k } k∈N 0 be a time-homogeneous Markov chain having the following properties: P2. An initial distribution (p i ) i∈S of {I k } k∈N 0 is determined by the probabilities p i = P(I 0 = i) that are positive for every i ∈ S. P3. A transition matrix P = p ij i,j ∈S of {I k } k∈N 0 is determined by the probabilities p ij = P (I k+1 = j |I k = i). The jump from I k−1 to I k can change the distribution of T k and/or X k at the moment A k only.
We will interpret {I k } k∈N 0 as 'switches' (see Gajek and Rudź (2017), p. 237 for a discussion on Property P3; a more detailed explanation of the considered model can be found in Gajek and Rudź (2018c)). Obviously, the amounts of claims and wait times are no longer assumed to be i.i.d. although some classical risk models can be derived. Let A, B ∈ F , and suppose that P(A) > 0. We will denote by P(B|A) the conditional probability of B given A; by EX, the expectation of a random variable X and by E[X|Z], the conditional expectation of X given the random variable Z. To simplify the notation, we will write: if p ij > 0 and 0 otherwise, where i, j ∈ S. With this notation, Δ −1 (δ n )(u) = (E 1 χ(D n (u)), ..., E s χ(D n (u))) for u 0 and F ij (x) = P ij (X 1 x), G ij (t) = P ij (T 1 t) for all i, j ∈ S and t, x ∈ R + .
Let R denote the set of all measurable functions defined on R 0 + and taking values in R 0 + . We use the symbol R s to denote the set {(ρ 1 , ..., ρ s ) : ρ i ∈ R for every i ∈ S}. We write the elements of R s in bold. Let ρ = (ρ 1 , ..., ρ s ) ∈ R s . Let us define operators i : R s → R by Let us denote = ( 1 , ..., s ). Given a non-decreasing harm function χ : R 0 + → R 0 + , let us define the risk operator L = (L 1 , ..., L s ) : R s → R s in the following way: for all i ∈ S and ρ ∈ R s . Throughout the paper, ∞ a means (a, ∞) for every a 0 and F (x) = 1 − F (x) for any distribution function F and every x 0.

Operator's formulae for the measures of risk
For every y 0, define W i n (u, y) = P i (D n (u) y), where i ∈ S, n ∈ N and u 0. Note that W i n (u, y), considered as a function of y, is the conditional distribution of D n (u) given the initial state i. Set: x ∈ R + and u, y ∈ R 0 + . In Theorem 1 below, we will assume that for all i, j ∈ S, k ∈ N 1 and t, x ∈ R + : A1. The conditional distribution of the random variables Z 2 , ..., Z k , given A simulation scheme for a regime-switching Sparre Andersen model, in which the above conditions hold, can be found in Gajek and Rudź (2018c).
, respectively, where δ i 0 (u) = 0 for all i ∈ S and u 0. To simplify the notation, we will denote χ(0) The next theorem provides a universal iterative method to establish the exact values of important insolvency risk measures. Exemplary 10 of them are listed in Section 1. Due to the flexibility of the considered probabilistic model, the methodology can be applied in several discrete and continuous time risk models as well. Numerical examples including some technical details on iterating the risk operator L are presented in Section 6.
Theorem 1 Assume that A1 and A2 hold. Let Δ : R 0 + → R 0 + be a strictly increasing function such that Δ(0) = 0. Assume that χ : R 0 + → R 0 + is a non-decreasing function such that χ(0) + = 0 and Eχ(D n ) < ∞ for every n ∈ N. If χ and W i n (·, y) have no common points of discontinuity for all i ∈ S and n ∈ N, then Proof We first show that, for all i ∈ S, n ∈ N 0 and u 0, it holds for all i ∈ S and u 0. Hence, by the law of total probability, Thus, Assumption A2 and the Fubini theorem imply Now, changing the integration variable and integrating by parts gives for all i ∈ S and u 0. Next, we will show that holds for all i ∈ S, n ∈ N and u, y ∈ R 0 + . Indeed, by the definitions of τ (u) and D n (u), In much the same way as above, We now apply Eq. 9 in order to prove (6). Indeed, integration by parts, Def. 3 and Eq. 9 imply altogether that Hence, the Fubini theorem and Def. 4 give for all i ∈ S, n ∈ N and u 0. By Eqs. 8 and 10, Eq. 6 holds for every n ∈ N 0 . Equivalently, Δ −1 (δ n+1 )(u) = LΔ −1 (δ n )(u). Repeating the same argument, we get Eq. 5, which completes the proof.
The next corollary will be used in the proof of Proposition 2 in Section 4. Let us denote P n,y (u) = (P 1 (D n (u) > y), ..., P s (D n (u) > y)), where P i (D n (u) > y) = P i (τ (u) n, −U (τ (u), u) > y), for n ∈ N, u 0 and y > 0. Clearly, P n,y (u) is the vector of the conditional probabilities that the deficit D n (u) is greater than y, provided that the initial state of the Markov chain is i = 1, ..., s.
In particular, when
Let us define and where i ∈ S and r ∈ R. Clearly, (M 1 (r), ..., M s (r)) is a vector-type counterpart of a one-dimensional modified moment generating function.
Given an arbitrary r > 0, let us denote now R i 0 (u, r) = e −ru and R 0 (u, r) = (R 1 0 (u, r), ..., R s 0 (u, r)) , where i ∈ S and u 0. For each i ∈ S, define iteratively a sequence {R i n (·, ·)} n∈N by (u, r), ..., R s n (u, r)). The issue of finding estimates for ruin-related quantities is important for many probabilists and practitioners. The regime-switching methodology results in several new bounds for multidimensional insolvency risk measures. The following proposition gives a twosided estimate for Δ −1 (δ i n )(u) in terms of m * (r) = min{m i (r) : i ∈ S} and M * (r) = max{M i (r) : i ∈ S}, where m i and M i are defined above. Example 4 shows that the resulting bound can be attainable.
Proposition 1 Let the assumptions of Theorem 1 hold. Then for all i ∈ S, n ∈ N and u 0.
Proof By Lemma 1, We will show first that the following inequalities where i ∈ S, r ∈ R + , u 0, hold for every n ∈ N. Indeed, so Inequalities 18 hold for n = 1. Assume now that they hold for some n ∈ N. By Eq. 11, one can verify that so Inequalities 18 hold for n + 1. Thus, by the induction principle, Formula 18 is valid for every n ∈ N. By Eq. 17 and Inequalities 18, Since r ∈ R + is arbitrary, Inequalities 16 hold as well.
Recall that F (x) = 1 − F (x) for any distribution function F and every x 0. The distribution F will be called new better than used in expectation with respect to the harm function χ (in short, χ -NBUE) if F (0) = 0 and the following inequality holds: for every x ∈ R + , provided that the integral on the right side of Inequality 19 is finite. Similarly, the distribution F will be called new worse than used in expectation with respect to the harm function χ (in short, χ -NWUE) if F (0) = 0, the opposite inequality in Bound 19 holds and its left side is finite.
Clearly, if F is NBU, it is χ -NBUE as well. Exponential distributions are both NBU and NWU. An arbitrary mixture of exponential distributions is NWU, but not NBU.
Focusing on ruin probabilities, one can show (see, e.g., Gajek and Rudź (2018b)) that the regime-switching Sparre Andersen model can be applied in the context of the EU directive Solvency II. In the case of general risk measures, the methodology is more complex. However, the vector Ψ n can be used to bound δ n . We provide now an upper and lower estimate of the insolvency measure δ i n in terms of the ruin probability Ψ i n . To this end, we will use the lower and the upper envelope distribution F * (y) = min {F ij (y) : (i, j ) ∈ S 2 } and F * (y) = max {F ij (y) : (i, j ) ∈ S 2 }, respectively.
Proposition 2 Let the assumptions of Theorem 1 hold.
Before giving the proof, let us consider an example showing that the bounds in Inequalities 20 and 21 are attainable.
Example 1 Assume that F ij are exponential with the same scale parameter, but the wait times have different distributions G ij depending on i, j ∈ S. Then both Inequalities 20 and 21 become equalities.
Another example confirming the precision of Inequality 21 is given below (see Example 2).
Proof of Proposition 2.
(i) We first show that Inequality 20 holds for n = 1. Indeed, Eq. 7 and the χ -NBUE assumption imply that Assume that Inequality 20 holds for some n ∈ N. We will show that it holds for n+1 as well. Indeed, Theorem 1, Def. 4, Eq. 8, Inequality 24 and Corollary 1 imply altogether that (ii) In much the same way as above, from Eq. 7 and the χ -NWUE assumption, we deduce that Theorem 1, Eq. 25, Inequality 27 and Corollary 1 imply altogether that By the induction principle, Inequalities 20 and 21 hold for every n ∈ N.
We will show now that Inequality 21 can be precise also for the distributions other than the one investigated in Example 1.
Example 2 Let us consider the classical discrete time risk model, in which the aggregated amount of premiums received each period is denoted by γ . Let F stands for the distribution function of the sums of claims in successive time periods. We assume that F is a mixture of two exponential distributions with scale parameters (0.5; 1), weights (0.9; 0.1), respectively, and γ = 2.1.
Assume first that Δ(x) = x and χ(x) = x, x ∈ R 0 + . Figure 1 shows that the insolvency risk measure δ 1 can be effectively bounded from below by Inequality 21. Figure 2 confirms the precision of Inequality 21 when Δ(x) = x 2 , x ∈ R 0 + .
generalizes the notion of the insolvency coefficient (or adjustment coefficient) which is well-known in a one-dimensional case. Theorem A.1 of Gajek and Rudź (2018a, p. 52) gives a sufficient condition for the existence of positive solutions of System 28. Let m i (r) be defined by Eq. 15 and m * (r) = max{m i (r) : i ∈ S}. Define μ i (r) = m * (r) 1 I (0, r * ) (r) + M i (r) 1 I [r * , ∞) (r), i ∈ S, r > 0, where r * = min{r i 0 : i ∈ S}. Corollary 1 of Gajek and Rudź (2018b) concerns the following inequality: for i ∈ S, n ∈ N, u 0. H. U. Gerber proved an upper bound to a one-dimensional finitehorizon ruin probability using a martingale approach for the Cramér-Lundberg model (see Gerber (1979), p. 139 for details). Grandell (1992) obtained an extension of this bound in ordinary renewal model. Gerber's bound finds some applications. For example, Centeno (2002) used it to investigate optimal reinsurance contracts. In Gajek and Rudź (2018b), it was proven that Inequality 29 is a generalization and improvement of Gerber's inequality in a regime-switching framework.
Corollary 2 Let the assumptions of Theorem 1 hold. Assume that F ij is χ -NBUE for all i, j ∈ S and the equations in System 28 have positive solutions. Then for all i ∈ S, n ∈ N and u 0.
Example 3 (Expected loss). Assume that F ij is NBUE for all i, j ∈ S. Under the assumptions of Proposition 2 (i) with Clearly, the above bound is attainable. Under the assumptions of Corollary 2, it also holds for all i ∈ S, n ∈ N and u 0.
For every i ∈ S, let us define now iteratively a sequence { R i n (·, ·)} n∈N by R i n (u, r) = L i R n−1 (u, r), n ∈ N, r > 0, u 0 , where R 0 (u, r) = R 0 (u, r), R n (u, r) = ( R 1 n (u, r), ..., R s n (u, r)) and L i is defined in Corollary 1. Recall that m i (r) = M i (r) − m i (r) and m * (r) = min{m i (r) : i ∈ S}. Proposition 1 generalizes the following two-sided bound on Ψ i n : where i ∈ S, n ∈ N and u 0. It originates from Theorem 4 of Gajek and Rudź (2018b), from which the idea of the proof of Proposition 1 also comes. Proposition 2 implies the following estimates of the insolvency measure δ i n : Corollary 3 Let the assumptions of Theorem 1 hold.
for all i ∈ S, n ∈ N and u 0.
for all i ∈ S, n ∈ N and u 0.
Proof Apply Inequalities 32 and Proposition 2.
In general, the bounds in Corollary 3 require one more assumption than estimates resulting from Proposition 1. On the other hand, the form of the associated risk operator L in Inequality 33 is slightly simpler than L used in Inequalities 16. We will present how to compute both these bounds in the following example.
Example 4 We will show first that both the upper bounds (16) and (33) are attainable even for n = s = 1 and an arbitrary Δ satisfying the assumptions of Theorem 1. Consider the model discussed in Example 2. Let β > 0 and F (x) = (1 − e −βx ) 1 I (0, ∞) (x) stands for the distribution function of the sums of claims in successive time periods.
By Proposition 1, and, by Corollary 3, Assume that γ = 0.40, u = 0.15, χ(x) = x, x ∈ R 0 + and β = 6. In this case so both the upper bounds resulting from Proposition 1 and Corollary 3 are attainable to at least three decimal places no matter what Δ has been chosen.
Consider now n = 2 in the above model. In this case  (16) and (33) are attainable to at least three decimal places for every Δ that satisfies the assumptions of Theorem 1.

Theoretical studies -a summary
Let us briefly summarize our main results and methods. We based our research methodology on a vector-valued risk operator L = (L 1 , ..., L s ) defined by Eq. 4, see Section 2 for details. In Section 3, we proved the following recursive vector-valued operator formula: which reflects a close connection between the iterations of the vector-valued risk operator L = (L 1 , ..., L s ) and the vector δ n (see Theorem 1 for details and Section 6 for numerical examples). Recall that Δ(L) and Δ −1 (δ n ) denote (Δ(L 1 ), ..., Δ(L s )) and (Δ −1 (δ 1 n ), ..., Δ −1 (δ s n )), respectively. In Section 4, we used Eq. 35 to prove a few upper and/or lower estimates of δ n . Whenever p ij = P (I k+1 = j |I k = i) is positive, we defined F ij (x) = P (X 1 x|I 0 = i, I 1 = j ) and G ij (t) = P(T 1 t|I 0 = i, I 1 = j) for all i, j ∈ S and t, x ∈ (0, ∞), see Section 2 for details. Recall that m i (r), M i (r) and m i (r) are given by Eqs. 13, 14 and 15, respectively.
Let R i 0 (u, r) = e −ru and R i n (u, r) be the ith coordinate of the nth iteration (see Section 4 for details) of L on R 0 (u, r) = (R 1 0 (u, r), ..., R s 0 (u, r)), r > 0. The following two-sided inequality holds for all i ∈ S, n ∈ N and u 0: where M * (r) = max{M i (r) : i ∈ S} and m * (r) = min{m i (r) : i ∈ S}. Example 4 clearly highlights that Inequality 36 can be attainable.
We have also used a kind of the NBUE (respectively, NWUE, cf. Brown and Li (2018)) property (see Inequality 19) to prove a few upper (respectively, lower) bounds on δ n . To this end, we defined the lower and the upper envelope distribution F * (x) = min {F ij (x) : (i, j ) ∈ S 2 } and F * (x) = max F ij (x) : (i, j ) ∈ S 2 , respectively. Surprisingly, δ i n can be estimated from above (respectively, from below) in terms of the probability of ruin Ψ i n and a linear functional of the envelope distribution F * (respectively, F * ). In Proposition 2 (i), we have shown that the following sharp inequality holds: Proposition 2 (ii) provides the corresponding lower bound on δ n in the NWUE case.
Assume that there exists a vector (r 1 0 , ..., r s 0 ) formed by positive solutions of the following equations: M i (r i 0 ) = 1, i ∈ S. Recall that r * = min{r i 0 : i ∈ S} and m * (r) = max{m i (r) : i ∈ S}. Proposition 2 (i) implies the following inequality: where μ i (r) = m * (r) 1 I (0, r * ) (r) + M i (r) 1 I [r * , ∞) (r), i ∈ S, r > 0 (see Corollary 2 for details). Inequality 37 is based on a generalized Gerber's bound. In Corollary 3 (i), another upper bound on δ i n is given under a kind of NBUE assumption. Corollary 3 (ii) provides a corresponding lower bound under the NWUE assumption.

Numerical studies
Using a fast algorithm of numerical iterations of L, we will show now how to calculate exact values of δ n in a Markov-switching framework.
Theorem 1 gives an exact recursive formula for the vector δ n of insolvency risk measures in a regime-switching model. It involves the vector-valued risk operator L. However, iterating L entails some numerical difficulties. To compute any coordinate of an arbitrary L's iteration, one has to calculate all coordinates of each previous iteration. Moreover, all the iterations of L should be properly discretised in a specific way which is discussed below. These computational complications arise even when calculating ruin probabilities (see, e.g, Gajek and Rudź (2018b)). Surprisingly, the problem of iterating the risk operator is complex even for risk models without a switch, cf. Gajek and Rudź (2013). We will discuss it for an arbitrary general insolvency risk measure δ n within the following discrete time regime-switching model.
Let a positive real m be such that P(T k = m) = 1 for each k ∈ N. We interpret T 1 , T 2 , ... as a fixed time period m (e.g. a quarter) and non-negative (a.s.) X k -as sum of the claims during the kth period. Recall that c is a known positive function defined on the state space S. Then a positive random variable γ k = c(I k−1 )m can be interpreted as the total amount of premiums during the kth period. Let us denote γ (·) = c(·)m. Then, the total amount of premiums during the first period, γ 1 , given the initial state i, equals a positive real γ (i).
The risk operator associated with the above risk process can be written as follows: and for all i ∈ S, ρ= (ρ 1 , ..., ρ s ) ∈ R s and u 0.
To calculate general insolvency risk measures resulting from Eq. 5, it is necessary to iterate L first. We will present now a numerical procedure for it. First of all, we will show how to iterate L on an arbitrary function ρ 0 = (ρ 0 1 , ..., ρ 0 s ) ∈ R s . Then, we will illustrate the use of Eq. 5 to calculate δ n .
Assume that each coordinate ρ 0 1 (u), ..., ρ 0 s (u) of the initial approximation ρ 0 (u) is discretised 1 with a fixed step h > 0 on a sufficiently large interval [0, u ], u kγ * . To calculate the first iteration ρ 1 , one has to discretise each of its s coordinates: on the interval [0 , u − γ * ] with the same step h. In much the same way as above, one calculates ρ 2 , ..., ρ k−1 . Finally, to determine the kth iteration ρ k of L on ρ 0 , one has to discretise the following s coordinates: on the interval [0 , u − kγ * ] with the step h.
In general, to calculate ρ k , one has to compute and store in a file all s coordinates of each previous iteration: ρ k−1 , ρ k−2 , ..., ρ 0 . Other numerical complications may arise from the argument u + γ (i) − x of the above integrands. Namely, the estimates determined in each previous iteration should be appropriately discretised (at the same points) on the interval larger by γ (i) (or by γ * , in the worst case) than in the current one. It means that the initial interval [0, u ] should be sufficiently large.
We illustrate the above methodology within the following example.