Eric C K Cheung
email: eric.cheung@unsw.edu.au

Landy Rabehasaina

Jae-Kyung Woo

Ran Xu
email: ran.xu@concordia.ca

Asymptotic Correlation Structure of Discounted Incurred But Not Reported Claims under Fractional Poisson Arrival Process

Keywords: Applied probability, Fractional Poisson process, Incurred But Not Reported (IBNR) claims, Infinite server queues, Correlation

This paper studies the joint moments of a compound discounted renewal process observed at different times with each arrival removed from the system after a random delay. This process can be used to describe the aggregate (discounted) Incurred But Not Reported claims in insurance and also the total number of customers in an infinite server queue. It is shown that the joint moments can be obtained recursively in terms of the renewal density, from which the covariance and correlation structures are derived. In particular, the fractional Poisson process defined via the renewal approach is also considered. Furthermore, the asymptotic behaviour of covariance and correlation coefficient of the aforementioned quantities is analyzed as the time horizon goes to infinity. Special attention is paid to the cases of exponential and Pareto delays. Some numerical examples in relation to our theoretical results are also presented.

Introduction

To model Incurred But Not Reported (IBNR) claims in insurance, in this paper it is first assumed that the claim arrivals follow a renewal process and the reporting delay for each claim is arbitrary. It is known (e.g. [START_REF] Willmot | A queueing theoretic approach to the analysis of the claims payment process[END_REF], [START_REF] Mikosch | Non-Life Insurance Mathematics: An Introduction with the Poisson Process[END_REF] and [START_REF] Ross | Introduction to Probability Models[END_REF]) that unreported claims in actuarial science can be connected to quantities considered in other fields such as the number of customers or particles in a system with an infinite server queue structure (i.e. customers are served immediately upon their arrivals by one of those many servers). There is a vast literature on systems involving infinite server queues. See e.g. [START_REF] Brown | Some results for infinite server queues[END_REF], [START_REF] Keilson | M/G/∞ with batch arrivals[END_REF], [START_REF] Liu | On the GI X /G/∞ system[END_REF], and [START_REF] Keilson | Networks of non-homogeneous M (t)/G/∞ systems[END_REF] for some classical results; and [START_REF] Ridder | Importance sampling algorithms for first passage time probabilities in the infinite server queue[END_REF], [START_REF] Pang | Infinite-server queues with batch arrivals and dependent service times[END_REF], [START_REF] Blom | A large-deviations analysis of Markov-modulated infinite-server queues[END_REF], [START_REF] Jansen | A large deviations principle for infinite-server queues in a random environment[END_REF], [START_REF] Moiseev | Queueing network MAP -(GI/∞) K with high-rate arrivals[END_REF], and [START_REF] Blom | Refined large deviations asymptotics for Markov-modulated infinite-server systems[END_REF] for some recent development of the subject. As mentioned in Mikosch (2009, Section 8.2.4), the counting process of the IBNR claims can be viewed as a generating model for the activities of processing packets in a large data network. In this case, the number of active sources at time t is the subject of interest of that model. Aggregate IBNR claim process is also related to the cumulative shock model with delayed termination (e.g. Finkelstein and Cha (2013, Section 4.4)), where the reporting time of a claim (i.e. the claim arrival time plus its reporting lag) corresponds to the time of the effective event which terminates fault. Interested readers are also referred to [START_REF] Blanchet | A heavy traffic approach to modeling large life insurance portfolios[END_REF] for the application of infinite serve queues in large insurance portfolios. See also Badescu et al. (2016a,b) for the theoretical properties and parameter estimations of a marked Cox model for IBNR claims.

In the analysis of such stochastic models, the knowledge of correlation structure of the underlying arrival process is often useful to understand the joint behaviour of the process of our interest observed at different times. Indeed, the correlation between two instants s and t when t → ∞ can be used to assess if the process has a long memory, i.e. if the state of the process at a given time has significant impact on its state at a much later time, and can serve to quantify such impact. This kind of information can be of paramount importance for decision making. In particular, the notion of long-range dependence appears when the correlation function exhibits the form of power-law decay. This definition arises in different fields such as finance, data network, and earthquake modelling. From a statistical point of view, the theoretical results allow the possibility of fitting a model to physical phenomena where longrange dependence or short-range dependence is practically observed, as explained in e.g. in Mikosch (2009, Section 8.2.4). In this paper, we shall focus on studying the covariance and correlation structures of compound renewal sums at different times (in the presence of a random delay in realization) as well as their asymptotic behaviours. To this end, the related joint moments are first derived, which may be of interest in its own right as they can in principle be used in moment-based approximations. It should be also noted that the joint moments and covariance structure of the aggregate discounted claims without any reporting delay were studied by Léveillé andAdékambi (2010, 2012) and [START_REF] Léveillé | A compound renewal model for medical malpractice insurance[END_REF].

To model the long-range dependence property for the renewal process governing the claim arrival process, a fractional generalization of the Poisson process, known as the fractional Poisson process, will be utilized. It is known that there are different ways to define the fractional Poisson process such as using time-changed processes (e.g. [START_REF] Leonenko | Correlation structure of timechanged Lévy processes[END_REF]). Here we adopt the definition obtained by a renewal-type treatment (e.g. [START_REF] Mainardi | A fractional generalization of the Poisson Processes[END_REF], and [START_REF] Meerschaert | The fractional Poisson process and the inverse stable subordinator[END_REF], where the fractional Poisson process is viewed as a renewal process with Mittag-Leffler interarrival times. Consequently, the interarrival times are heavy-tailed and have infinite mean (e.g. [START_REF] Repin | Fractional Poisson law[END_REF]) as opposed to the light-tailed exponential interarrival times with finite mean in the Poisson process. The fractional Poisson process is an interesting choice in the context of actuarial science, and in particular it can be used to model the arrival of claims caused by rare and extreme events such as storms and high-magnitude earthquakes in line with the discussions in e.g. [START_REF] Benson | Recurrence of extreme events with power-law interarrival times[END_REF] and [START_REF] Biard | Fractional Poisson process: Long-range dependence and applications in ruin theory[END_REF]. Indeed, [START_REF] Benson | Recurrence of extreme events with power-law interarrival times[END_REF] commented that the use of this type of renewal process is a critical extension from the Poisson assumption because the application of a Poisson model for a geologic process with heavy-tailed waiting times between events will result in a significant misrepresentation of the associated risk. Our analysis will shed light on how this heavy-tail feature impacts the long-term correlation of the IBNR process.

In what follows, we shall define the model in the context of actuarial science. It is assumed that the claim counting process {N (t)} t≥0 is a renewal process with the sequence of arrival times {T n } ∞ n=1 where T n := n i=1 τ i for n ∈ N. Here N is the set of positive integers, and {τ i } ∞ i=1 represents the sequence of interarrival times that are independent and identically distributed (iid) with cumulative distribution function (cdf) F (t) and probability density function (pdf) f (t) = F (t). We adopt the usual convention that T 0 := 0. Therefore, the claim count is defined by N (t) := max{i ∈ N ∪ {0} : T i ≤ t} with N (0) := 0. The renewal function of the claim counting process shall be denoted by m(t) := E[N (t)] = ∞ i=1 F * (i) (t), where F * (i) (t) is the i-th fold convolution of the cdf F (t) with itself and f * (i) (t) := dF * (i) (t)/dt is its pdf. Let us introduce {X i } ∞ i=1 as the sequence of claim amounts which are also assumed to be iid. Denoting L i as the time-lag (or reporting delay) corresponding to the i-th claim arrival, the time-lags {L i } ∞ i=1 are assumed to form an iid sequence with common cdf W (t) := 1 -W (t) and pdf w(t) := W (t). It is assumed that the claim counting process {N (t)} t≥0 (or the interarrival times {τ i } ∞ i=1 ), the claim amounts {X i } ∞ i=1 and the delays {L i } ∞ i=1 are mutually independent. Such independence assumption is standard in the actuarial literature (see e.g. [START_REF] Karlsson | A stochastic model for time lag in reporting of claims[END_REF] and [START_REF] Willmot | A queueing theoretic approach to the analysis of the claims payment process[END_REF]), and we refer interested readers to [START_REF] Landriault | Analysis of IBNR claims in renewal insurance models[END_REF] for some discussions on possible dependence assumptions on the triplet (τ i , X i , L i ). For later use, the generic random variables of τ i , X i and L i are denoted by τ , X and L respectively. It will be seen that the moments of the claim amount are sufficient for our analysis, and therefore we define µ k := E[X k ] and µ 0 := 1 for convenience.

The total discounted IBNR claim process {Z δ (t)} t≥0 is now defined by

Z δ (t) := N (t) i=1 e -δ(T i +L i ) I {T i +L i >t} X i = ∞ i=1 e -δ(T i +L i ) I {T i ≤t<T i +L i } X i , (1.1)
where δ ≥ 0 is a constant force of interest. We shall adopt the convention that i j := 0 whenever i < j. For notational simplicity, we let Z(t) = Z 0 (t) when δ = 0. Note that, although the model is explained in an actuarial setting, it can be applied to a queueing context via a switch of terminology by changing 'claim arrival' and 'reporting delay' to 'customer arrival' and 'service time' respectively. In particular, if the distribution of X is assumed to be a point mass at one, then Z(t) can be interpreted as the total number of customers at time t in a G/G/∞ queue. (In such case, the assumption of independence between {N (t)} t≥0 and {L i } ∞ i=1 means that the service offered to customers does not depend on their arrival times, which is a standard assumption in G/G/∞ queue.) Consequently, our formulation under an actuarial context is more general than a G/G/∞ queue, as (i) insurance companies possibly take into account the effect of interest (allowing for δ > 0) when discounting claims; and (ii) actuaries associate a loss amount X to each claim event when calculating or estimating outstanding claims. In this paper, we are interested in the covariance and correlation structures of Z δ (s) and Z δ (t) for 0 < s ≤ t, which are defined as

Cov[Z δ (s), Z δ (t)] = E[Z δ (s)Z δ (t)]-E[Z δ (s)]E[Z δ (t)] and Corr[Z δ (s), Z δ (t)] = Cov[Z δ (s), Z δ (t)]/ Var[Z δ (s)]Var[Z δ (t)] respectively.
Concerning the renewal process for claim arrivals, in most of our analysis it is essential to know explicitly the renewal density m (t) for t > 0. The class of fractional Poisson processes turns out to be a good candidate thanks to its nice renewal density (see (4.1)). Certainly, some results will be simplified for the Poisson process as m (t) = λ > 0 does not depend on t.

Outline of results and long-range dependence of IBNR process

In order to have a global overview of the paper, our main contributions are presented here. In Section 3, the marginal moments of Z δ (t) and the joint moments of Z δ (s) and Z δ (t) for 0 < s ≤ t are derived in Theorems 1 and 2 using renewal arguments, and consequently the covariance is studied in Corollary 1 and Theorem 3. Section 4 is devoted to the asymptotic behaviours (as t → ∞) of the mean, variance and covariance in the case of fractional Poisson arrival process. To obtain simpler formulas, distributional assumptions on the random delays are made, where exact and asymptotic results are provided in Propositions 2-4 for the exponential case, and asymptotics in Propositions 5-7 for the Pareto case. We remark that the focus of the paper is to analyze the properties of the IBNR process and to provide probabilistic interpretations. Immediate practical applications of the results, including statistical estimation and fitting with a data set, are outside the scope of this paper and will be topics for future research. For legibility purpose, the asymptotic results obtained in the paper are summarized as follows when claims occur according to a fractional Poisson process with index α ∈ (0, 1). The notation g 1 (t) ∝ g 2 (t) means that g 1 (t) ∼ Cg 2 (t) as t → ∞ for some constant C > 0 that possibly depends on s but not on t. Moreover, L d = E(β) means that the delay L is exponentially distributed with mean 1/β whereas L d = Pareto(θ, η) means that the survival function of L is given by (4.27). We remark that the multiplicative constants are omitted in the summary below not only for the sake of brevity but also for the fact that the asymptotic correlation up to a multiplicative factor will be sufficient for determining whether a process possesses short-range or long-range dependence (see Definition 2). Nonetheless, the exact expressions for these multiplicative constants are given by (4.14), (4.19), (4.26), (4.32), (4.35) and (4.45) in the aforementioned propositions and can be computed or approximated (see Section 5).

L d = E(β) : E[Z δ (t)] ∝ e -δt t α-1 , Var[Z δ (t)] ∝ e -2δt t α-1 , Cov[Z δ (s), Z δ (t)] ∝ e -δt t α-2 (2.1) L d = Pareto(θ, η) : 0 < η < 1 η = 1 η > 1 E[Z(t)] ∝ t α-η ∝ t α-1 ln t ∝ t α-1 0 < η < α α ≤ η < 1 η = 1 η > 1 Var[Z(t)] ∝ t 2(α-η) ∝ t α-η ∝ t α-1 ln t ∝ t α-1 0 < η < 2 -α η ≥ 2 -α Cov[Z(s), Z(t)] ∝ t -η ∝ t α-2 (2.2)
We emphasize that our contributions are not only to actuarial science but also to queueing theory. In the queueing literature, only few papers have dealt with fractional Poisson arrivals because the techniques are usually very specific to fractional calculus. See e.g. [START_REF] Orsingher | On a fractional linear birth-death process[END_REF] for a particular queue with a birth and death structure modelled by fractional Poisson arrivals, which is one of the rare papers on the subject. However, there are no results in the literature concerning infinite server queues where the arrivals are modelled by a fractional Poisson process, to the best of our knowledge. The result (2.1) and the tables in (2.2) enable us to give some qualitative insight on whether the IBNR process (or the G/G/∞ queue as a special case) exhibits long-range dependence under the fractional Poisson setting. The behaviour will turn out to be different depending on the distribution of the delays, as displayed in (2.3)-(2.5) thereafter. The choice of the delays, namely exponential and Pareto, is motivated by the fact that these are representatives of light-tailed and heavy-tailed distributions respectively. These distributions were also chosen for technical purpose, and we however believe that this will shed light on whether behaviours resembling those in (2.3)-(2.5) are valid more generally for wider classes of light-tailed and heavy-tailed distributions in future research. Interested readers are also referred to [START_REF] Resnick | Self-similar communication models and very heavy tails[END_REF] for an infinite server queue with heavy-tailed service time and Poisson arrivals, where heavy-tailed service time was argued to be a good candidate for modelling the transfer of huge files across the World Wide Web. In the following, we first state the notion of long-range dependence defined by Maheshwari and Vellaisamy (2016, Section 2.2).

Definition 1 (Long-range dependence in Maheshwari and Vellaisamy (2016)) Let {Y (t)} t≥0 be a stochastic process. Suppose that the correlation function Corr[Y (s), Y (t)] satisfies, for all s > 0,

lim t→∞ Corr[Y (s), Y (t)] t -d = c(s)
for some d > 0 and c(s) > 0. If d ∈ (0, 1], then {Y (t)} t≥0 is said to possess long-range dependence property. On the other hand, if d ∈ (1, 2), then {Y (t)} t≥0 has short-range dependence.

From Leonenko et al. (2014, p.10), it is known that the fractional Poisson process with index α ∈ (0, 1) has long-range dependence property. However, the above Definition 1 only includes processes where the correlation function exhibits power decay (and it excludes the cases where the correlation decreases exponentially or follows other decay law). According to Mikosch (2009, p.283), one also has the following more general definition of long-range dependence, which is adapted to stochastic processes that are not necessarily stationary.

Definition 2 (General notion of long-range dependence) Let {Y (t)} t≥0 be a stochastic process. Then {Y (t)} t≥0 is said to possess long-range dependence property if one has for all s > 0 that

∞ 0 |Corr[Y (s), Y (t)]|dt = ∞.
On the other hand, {Y (t)} t≥0 has short-range dependence if, for all s > 0,

∞ 0 |Corr[Y (s), Y (t)]|dt < ∞.
For convenience, from now on long-range dependence and short-range dependence will be abbreviated as LRD and SRD respectively. It is instructive to note that if a process has the LRD (resp. SRD) property under Definition 1, then it also has the LRD (resp. SRD) property under Definition 2. In the upcoming discussion, Definition 2 will be adopted. In the case where the delays follow exponential distribution, one sees from (2.1) that

L d = E(β) : Corr[Z δ (s), Z δ (t)] ∝ t -3-α 2 , (2.3)
and therefore {Z δ (t)} t≥0 has SRD. When the delays are Pareto(θ, η) distributed, the asymptotic behaviour of the correlation can be readily obtained using (2.2), resulting in the following table.

L d = Pareto(θ, η) : 0 < η < α α ≤ η < 1 η = 1 1 < η ≤ 3-α 2 3-α 2 < η < 2 -α η ≥ 2 -α Corr[Z(s), Z(t)] ∝ t -α ∝ t -(η+α)/2 ∝ t -(1+α)/2 √ ln t ∝ t -η-(α-1)/2 ∝ t -η-(α-1)/2 ∝ t -(3-α)/2 Dependence LRD LRD LRD LRD SRD SRD (2.4)
The above results can be compared to the situation where arrivals occur according to a Poisson process (i.e. α = 1). We only look at the discount-free case (i.e. δ = 0) as follows. When L d = E(β), we use the covariance in (3.8) along with the variance in (3.11) to easily check that Corr[Z(s), Z(t)] ∝ e -βt in the Poisson case, and therefore the process {Z(t)} t≥0 has SRD. This is consistent with the fractional Poisson case (i.e. α ∈ (0, 1)) described in (2.3). In contrast, when L d = Pareto(θ, η), we omit the details and state that (3.8) implies Cov[Z(s), Z(t)] ∝ t -η in the Poisson case. Moreover, one can use (3.11) with the help of L'Hôpital's rule to verify that Var

[Z(t)] ∝ t -η+1 if 0 < η < 1; Var[Z(t)] ∝ ln t if η = 1;
and Var[Z(t)] ∝ 1 if η > 1, and consequently the following table is constructed.

L d = Pareto(θ, η) ; α = 1 : 0 < η < 1 η = 1 η > 1 Corr[Z(s), Z(t)] ∝ t -(η+1)/2 ∝ t -1 √ ln t ∝ t -η Dependence LRD LRD SRD (2.5)
Comparing (2.4) and (2.5), it is interesting to note that {Z(t)} t≥0 is SRD in the Poisson case but LRD in the fractional Poisson case when the Pareto delay parameter η is such that 1 < η ≤ 3-α 2 . For other values of η the dependence property in both cases are in agreement.

These correlation behaviours can be further interpreted as follows. When the delay distribution is exponential which is light-tailed, a claim that has occurred but has not been reported at time s is likely to be reported by time t > s when t grows large, regardless of whether the claim arrives according to a Poisson or a fractional Poisson process. This explains the SRD property. On the other hand, switching from exponential to Pareto delay which is heavy-tailed, an unreported claim at time s is more likely to remain unreported at time t > s (for large t) because E(β) ≤ st,∞ Pareto(θ, η). Here Y 1 ≤ st,∞ Y 2 means that P(Y 1 > x) ≤ P(Y 2 > x) for x large enough, which is a generalization of the classical stochastic order for random variables (see e.g. Shaked and Shanthikumar (2007, Chapter 1)). Consequently, the process {Z(t)} t≥0 may have a tendency to be LRD. Note that increasing the Pareto shape parameter η leads to stochastically smaller delay (in the sense that Pareto(θ 1 , η 1 ) ≤ st,∞ Pareto(θ 2 , η 2 ) for η 1 > η 2 ), which explains the observation from (2.4) and (2.5) that the asymptotic correlation between Z(s) and Z(t) decreases faster (as a function of t) as η increases, and {Z(t)} t≥0 eventually switches from LRD to SRD once η exceeds a certain threshold. In particular, for Poisson claim arrivals (see (2.5)), it is noticed that {Z(t)} t≥0 is LRD if and only if 0 < η ≤ 1, i.e. when the Pareto delay has infinite mean. For fractional Poisson claim arrivals with index α ∈ (0, 1), the interarrival times are heavy-tailed (see (4.2)) and have infinite mean, and hence it is natural that there is a competition between the index α for the interarrival time and the Pareto parameter η, and the threshold in this case is given by 3-α 2 as in (2.4).

Remark 1 Instead of analyzing the behaviour of Corr[Z(s), Z(t)] as t → ∞ (for fixed s > 0), it will also be interesting to study Corr[Z(s), Z(t)] as both s and t tend to infinity. For example, one may replace s and t by t and t + h respectively and look at Corr[Z(t), Z(t + h)], and our exact results (e.g. Corollary 1 and Equation (4.24)) are still valid, from which asymptotics as t → ∞ (for fixed h > 0) may be obtained. However, we leave this as future research because it is the asymptotics of Corr[Z(s), Z(t)] as t → ∞ that help us determine whether the IBNR process is SRD or LRD.

General results

In this section, general recursive formulas for the joint moments of the total discounted IBNR claims Z δ (s) and Z δ (t) are derived without any specific distributional assumptions on the interarrival times {τ i } ∞ i=1 , the reporting delays {L i } ∞ i=1 or the claim amounts {X i } ∞ i=1 (recall that the moments of the claim amounts suffice in all the analysis). Our results can be conveniently expressed in terms of the Dickson-Hipp (D-H) operator T (see [START_REF] Dickson | On the time to ruin for Erlang(2) risk processes[END_REF]). For any integrable real function g (with non-negative domain), the D-H operator is defined by T u g(v) := ∞ v e -u(t-v) g(t)dt with both u and v non-negative. We begin by looking at the expectation in the following proposition.

Proposition 1 For a renewal claim arrival process, the mean of the total discounted IBNR claims until time t > 0 is given by

E[Z δ (t)] = µ 1 e -δt t 0 T δ w(t -x)dm(x). (3.1)
In particular, when δ = 0 this reduces to

E[Z(t)] = µ 1 t 0 W (t -x)dm(x). (3.2)
Proof: By taking expectation on (1.1) with the use of independence assumptions, one finds

E[Z δ (t)] = E ∞ i=1 e -δ(T i +L i ) I {T i ≤t<T i +L i } X i = ∞ i=1 E e -δ(T i +L i ) I {T i ≤t<T i +L i } X i = ∞ i=1 µ 1 t 0 e -δx E e -δL i I {L i >t-x} dF * (i) (x) = µ 1 t 0 e -δx ∞ t-x e -δy dW (y) d ∞ i=1 F * (i) (x) = µ 1 t 0 e -δx ∞ t-x
e -δ(y-(t-x)) e -δ(t-x) w(y)dy dm(x), which simplifies to give (3.1). Then (3.2) follows immediately as T 0 w(x) = W (x).

This direct technique applied in Proposition 1, however, is not very helpful to derive the expression for higher-order moments. For higher moments, it is useful to employ the conditioning technique based on the first claim arrival time T 1 . It is thus convenient to define the new notation

Z 2,δ (t) := N (t) i=2 e -δ(T i -T 1 +L i ) I {T i +L i >t} X i .
Conditioning on T 1 = x ≤ t, we note that Z 2,δ (t) is also a discounted IBNR claim that has the same distribution as Z δ (t -x). Then Z δ (t) may be expressed in terms of Z 2,δ (t) as

Z δ (t) = e -δT 1 I {T 1 ≤t} e -δL 1 I {T 1 +L 1 >t} X 1 + Z 2,δ (t) . (3.3)
Subsequently, recursive formula for the marginal moments of Z δ (t), namely E[Z n δ (t)], can be derived for n ∈ N via a binomial expansion in Theorem 1. We remark that the same result can also be obtained by considering the moment generating function of Z δ (t) followed by differentiation: such a technique was used by [START_REF] Woo | On multivariate discounted compound renewal sums with time-dependent claims in the presence of reporting/payment delays[END_REF] who considered the joint moments of discounted compound renewal sums for k dependent business lines (or k different types of dependent claims) evaluated at the same time point t. However, our technique of using a binomial expansion is easily applicable in the present context for the joint moments of the discounted IBNR process (for a single line of business) evaluated at different time points s and t, namely Z δ (s) and Z δ (t), as in the proof of Theorem 2.

Theorem 1 For n ∈ N and t > 0, a recursive formula for evaluating E[Z n δ (t)] is given by

E[Z n δ (t)] = n-1 i=0 µ n-i n i e -(n-i)δt t 0 e -iδx T (n-i)δ w(t -x)E[Z i δ (t -x)]dm(x), (3.4)
where the starting value is

E[Z 0 δ (t)] = 1. Proof: See Appendix A.1.
Theorem 2 For 0 < s ≤ t and n, m ∈ N, the joint moments of Z δ (s) and Z δ (t) can be calculated recursively as

E[Z n δ (s)Z m δ (t)] = n-1 i=0 µ n-i n i e -(n-i)δs s 0 e -(m+i)δx E[Z i δ (s -x)Z m δ (t -x)]T (n-i)δ w(s -x)dm(x) + m-1 j=0 n i=0 µ n+m-i-j n i m j e -(n+m-i-j)δt s 0 e -(i+j)δx E[Z i δ (s-x)Z j δ (t-x)]T (n+m-i-j)δ w(t-x)dm(x), (3.5)
and this requires the application of Theorem 1 as a starting point.

Proof: See Appendix A.2.
Furthermore, from Theorems 1 and 2, the covariance formula for Z δ (s) and Z δ (t) is readily available as given below.

Corollary 1 For a renewal claim arrival process, the covariance of the total discounted IBNR claims at the time points s and t (where 0 < s ≤ t) is given by

Cov[Z δ (s), Z δ (t)] = µ 1 e -δs s 0 e -δx E[Z δ (t -x)]T δ w(s -x)dm(x) + µ 2 e -2δt s 0 T 2δ w(t -x)dm(x) + µ 1 e -δt s 0 e -δx E[Z δ (s -x)]T δ w(t -x)dm(x) -E[Z δ (s)]E[Z δ (t)]. (3.6)
If there is no discounting (i.e. δ = 0), this reduces to

Cov[Z(s), Z(t)] = µ 1 s 0 W (s -x)E[Z(t -x)] + E[Z(s -x)]W (t -x) dm(x) + µ 2 s 0 W (t -x)dm(x) -E[Z(s)]E[Z(t)]. (3.7)
Furthermore, for a Poisson process with rate λ > 0 (i.e. dm(x) = λdx), it is simplified to

Cov[Z(s), Z(t)] = λµ 2 s 0 W (t -x)dx, (3.8)
which is consistent with the covariance expression in Blom et al. (2014, Section 4).

Proof: As Cov[Z δ (s), Z δ (t)] = E[Z δ (s)Z δ (t)] -E[Z δ (s)]E[Z δ (t)],
setting n = m = 1 in Theorem 2 with a substitution of the result in Theorem 1 for s and t yields (3.6). When δ = 0, we have T 0 w(x) = W (x) and hence (3.7) is found. For a Poisson arrival process, by substituting E[Z(t)] = λµ 1 t 0 W (x)dx followed by interchanging the order of integral in the second term of (3.7), the first two terms can be rewritten as

µ 1 s 0 W (s -x)E[Z(t -x)] + E[Z(s -x)]W (t -x) dm(x) = λ 2 µ 2 1 s 0 W (x) t-s+x 0 W (y)dy + t t-s+x W (y)dy dx = λ 2 µ 2 1 s 0 W (x)dx t 0 W (y)dy ,
which is identical to the last term in (3.7). Therefore, (3.8) follows and the proof is complete.

The following corollary is immediately obtainable by letting s = t in Corollary 1.

Corollary 2 For a renewal claim arrival process, the variance of the total discounted IBNR claims at time t > 0 is given by

Var[Z δ (t)] = 2µ 1 e -δt t 0 e -δx E[Z δ (t -x)]T δ w(t -x)dm(x) + µ 2 e -2δt t 0 T 2δ w(t -x)dm(x) -{E[Z δ (t)]} 2 .
(3.9) When δ = 0, this reduces to

Var[Z(t)] = 2µ 1 t 0 W (t -x)E[Z(t -x)]dm(x) + µ 2 t 0 W (t -x)dm(x) -{E[Z(t)]} 2 .
(3.10)

Furthermore, for a Poisson arrival process with rate λ > 0, it is simplified to

Var[Z(t)] = λµ 2 t 0 W (x)dx. (3.11) Note that Var[Z(t)] → λµ 2 E[L] as t → ∞, where E[L] = ∞ 0 W (x)dx. Theorem 3 When δ = 0 and E[L] < ∞, the limit of Cov[Z(s), Z(t)] is given by Cov[Z(s), Z(t)] -→ 0 as t → ∞. (3.12) Proof: First, we use the result E[Z(∞)] := lim t→∞ E[Z(t)] = µ 1 E[L]/E[τ ]
which is obtainable by applying the Smith's Renewal Theorem to (3.2) as

E[Z(t)] -→ µ 1 E[τ ] ∞ 0 W (y)dy as t → ∞.
By taking the limit t → ∞ with the help of dominated convergence, the first term on the right-hand side of (3.7) in Corollary 1 becomes

lim t→∞ µ 1 s 0 W (s -x)E[Z(t -x)] + E[Z(s -x)]W (t -x) dm(x) = µ 2 1 E[L] E[τ ] s 0 W (s -x)dm(x),
where the second term in the integration converges to 0. This can be cancelled out with the third term in (3.7). Then, the desired result (3.12) follows.

The result in Theorem 3 shows that, provided that the delay time L is integrable, the IBNR claims are asymptotically uncorrelated, i.e. on a first approximation, one has

E[Z(t)Z(s)] ≈ E[Z(t)]E[Z(s)] as t → ∞.
More information on the speed of convergence of the covariance towards 0 will be given in the forthcoming Propositions 4 and 7 when additional distributional assumptions are made. In fact, one of the obstacles in the proofs in those propositions is to be able to provide a second order approximation for the joint moment E[Z(t)Z(s)] as t → ∞, in order to get some information on the rate of decrease of Cov[Z(s), Z(t)].

Fractional Poisson Process

In this section, we consider the fractional Poisson process for the claim counting process {N (t)} t≥0 , so that the renewal function and renewal density are given by, for t > 0,

m(t) = E[N (t)] = λt α Γ(1 + α) and m (t) = λαt α-1 Γ(1 + α) = λt α-1 Γ(α) , (4.1)
where Γ(x) is the (ordinary) Gamma function defined by Γ(x) = ∞ 0 y x-1 e -y dy for Re(x) > 0, and α and λ are parameters such that 0 < α ≤ 1 and λ > 0. From Equation (28) of [START_REF] Laskin | Fractional Poisson process[END_REF], its variance is given by

Var[N (t)] = m(t) + [m(t)] 2 αB(α, 1/2) 2 2α-1 -1 , where B(a, b) = Γ(a)Γ(b)/Γ(a + b) = 1 0 y a-1 (1 -y) b-1 dy is the Beta function.
Therefore, the variance is increasing in time according to a power law t 2α when 0 < α < 1. If α = 1, this process becomes the ordinary Poisson process with rate λ. From [START_REF] Repin | Fractional Poisson law[END_REF], the asymptotic behaviour for the survival function of the interarrival times satisfies (when 0 < α < 1)

1 -F (t) ∼ Ct -α as t → ∞, (4.2)
where C is some constant. For the estimation of the parameters α and λ of a fractional Poisson process, interested readers are referred to [START_REF] Cahoy | Parameter estimation for fractional Poisson processes[END_REF], who first derived the estimators by matching the first two moments of the log of the interarrival time τ and then proved their asymptotic normality.

In the sequel, two different distributional assumptions for the random delay variable L, namely exponential and Pareto, will be considered in detail. To proceed with our analysis, we first recall the definitions and properties of some special functions.

(i) The Gaussian (or ordinary) hypergeometric function 2 F 1 (a, b, c, z) is defined as the power series (for |z| ≤ 1)

2 F 1 (a, b, c, z) = ∞ n=0 (a) n (b) n (c) n z n n! , (4.3)
where (a) n = a(a + 1)

• • • (a + n -1)
for n ∈ N is the Pochhammer symbol with (a) 0 = 1. Its integral representation is given by

2 F 1 (a, b, c, z) = 1 B(b, c -b) 1 0 x b-1 (1 -x) c-b-1 (1 -zx) -a dx, Re(c) > Re(b) > 0, (4.4)
provided that the right-hand side converges. We also have the relationship between hypergeometric function and incomplete Beta function given by

B(a, b, x) = x 0 t a-1 (1 -t) b-1 dt = x a a 2 F 1 (a, 1 -b; a + 1; x). (4.5)
Another important identity for the Gaussian hypergeometric function is

2 F 1 (a, b; c; 1) = Γ(c)Γ(c -a -b) Γ(c -a)Γ(c -b) , Re(c) > Re(a + b). (4.6)
(ii) The confluent hypergeometric function of the first kind or the Kummer's function 1 F 1 (a, b, z) introduced by Kummer (1837) (see e.g. Abramowitz and Stegun (1972, Chapter 13)) is defined by

1 F 1 (a, b, z) = ∞ n=0 (a) n (b) n z n n! , (4.7)
and its integral representation is

1 F 1 (a, b, z) = 1 B(a, b -a) 1 0 e zx x a-1 (1 -x) b-a-1 dx, Re(b) > Re(a) > 0. (4.8)
We also have the transform relation between two Kummer's functions as

1 F 1 (a, b, z) = e z 1 F 1 (b -a, b, -z). (4.9)
The Kummer's function satisfies

1 F 1 (a, b, z) = Γ(b) Γ(a) e z z a-b {1 + O(|z| -1 )}, Re(z) > 0. ( 4

Exponential delay

In this section, it is assumed that W (x) = e -βx for x ≥ 0, where β > 0 is the exponential parameter for the delay L so that E[L] = 1/β and E[e -sL ] = β/(β + s).

Expectation of total discounted IBNR process. We first establish an expression for E[Z δ (t)] and provide its asymptotic behaviour as t → ∞.

Proposition 2 With exponential delay in a fractional Poisson claim arrival process, the mean of the total discounted IBNR claims Z δ (t) when t > 0 is given by, for δ = 0,

E[Z(t)] = µ 1 λe -βt t α Γ(1 + α) 1 F 1 (α, 1 + α, βt) = µ 1 e -βt 1 F 1 (α, 1 + α, βt)m(t), (4.11)
and for δ > 0,

E[Z δ (t)] = µ 1 λE[e -δL ]e -(β+δ)t t α Γ(1 + α) 1 F 1 (α, 1 + α, βt) = µ 1 E[e -δL ]e -(β+δ)t 1 F 1 (α, 1 + α, βt)m(t), (4.12)
where m(t) is the renewal function given in (4.1) for the fractional Poisson process. In addition, the corresponding asymptotic behaviours are

E[Z(t)] ∼ µ 1 λ βΓ(α) t α-1 = µ 1 E[L]m (t) as t → ∞, (4.13)
and

E[Z δ (t)] ∼ µ 1 λE[e -δL ] βΓ(α) e -δt t α-1 = µ 1 E[L]E[e -δL ]e -δt m (t) as t → ∞, (4.14)
respectively.

Proof: First, we note that for exponential delay, one has

T δ w(x) = β β+δ e -βx = E[e -δL ]e -βx = E[e -δL ]W (x). In particular, T 0 w(x) = W (x) when δ = 0. Hence, T δ w(x) = E[e -δL ]T 0 w(x)
, and as a result of Proposition 1 we have the relation

E[Z δ (t)] = e -δt E[e -δL ]E[Z(t)]. (4.15)
In other words, it is immediate to obtain the mean in the presence of discounting from the mean without discounting. When δ = 0, we have from Proposition 1 that E[Z(t)] = µ 1 e -βt t 0 e βx dm(x). Then, with the fractional Poisson process for claim arrivals, substitution of its renewal density (4.1) followed by a change of variable y = x/t yields

E[Z(t)] = µ 1 λαe -βt Γ(1 + α) t 0 e βx x α-1 dx = µ 1 λαe -βt Γ(1 + α) t α 1 0 e (βt)y y α-1 dy . (4.16)
Utilizing the Kummer's function in (4.8), the above expectation may be neatly expressed as (4.11). Consequently, using the relation (4.15), we easily obtain (4.12). It is indeed useful to have expressed E[Z(t)] and E[Z δ (t)] in terms of 1 F 1 (a, b, z) as its asymptotic behaviour is well established in the literature. More precisely, with the help of (4.10), one finds (4.13) and (4.14).

Example 1 (Poisson process) For notational convenience, we denote the present value of a t-year continuous annuity payable at rate $1 per year as āt δ := t 0 e -δx dx = (1 -e -δt )/δ. When the claims arrive according to a Poisson process with rate λ, from Proposition 1 we have E[Z(t)] = µ 1 λā t β and thus E[Z δ (t)] = µ 1 λe -δt E[e -δL ]ā t β from (4.15). Of course, the same results are obtainable by setting α = 1 in Proposition 2, since when α = 1 the fractional Poisson process is just a Poisson process. In particular, noting that e -βt t 1 F 1 (1, 2, βt) = āt β from (4.7), the result follows from (4.11) for α = 1. Also, asymptotic behaviours are given as

E[Z(t)] ∼ µ 1 λE[L] and E[Z δ (t)] ∼ µ 1 λE[L]E[e -δL ] e -δt as t → ∞.
Variance of total discounted IBNR process. Next, from (3.9), we find the variance of Z δ (t) and its asymptotic result in the following proposition.

Proposition 3 With exponential delay in a fractional Poisson claim arrival process, the variance of the total discounted IBNR claims Z δ (t) when t > 0 is given by

Var[Z δ (t)] = 2µ 2 1 λ 2 α{E[e -δL ]} 2 e -2(β+δ)t [Γ(1 + α)] 2 ∞ n=0 A n t 2α+n 1 F 1 (α, 2α + n + 1, 2βt) + µ 2 λE[e -2δL ]e -(β+2δ)t t α Γ(1 + α) 1 F 1 (α, 1 + α, βt) - µ 2 1 λ 2 {E[e -δL ]} 2 e -2(β+δ)t t 2α [Γ(1 + α)] 2 [ 1 F 1 (α, 1 + α, βt)] 2 , (4.17)
where

A n = αβ n n!(α + n) B(α, α + n + 1). (4.18)
In addition, the asymptotic behaviour when 0 < α < 1 is

Var[Z δ (t)] ∼ µ 2 1 λ 2 {E[e -δL ]} 2 β α+1 Γ(α) + µ 2 λE[e -2δL ] βΓ(α) e -2δt t α-1 as t → ∞. (4.

19)

Proof: The three terms on the right-hand side of (3.9) are expressed as follows. The first term can be calculated by using (4.12), T δ w(x) = β β+δ e -βx = E[e -δL ]e -βx and (4.1) as

2µ 1 e -δt t 0 e -δx E[Z δ (t -x)]T δ w(t -x)dm(x) = 2µ 2 1 λ 2 α{E[e -δL ]} 2 [Γ(1 + α)] 2 e -2δt t 0 e -2βx x α (t -x) α-1 1 F 1 (α, 1 + α, βx)dx. (4.20)
Using the infinite series representation of the Kummer's function in (4.7) with a change of variable from x/t to y, the integral term above can be expressed as

t 0 e -2βx x α (t -x) α-1 1 F 1 (α, 1 + α, βx)dx = ∞ n=0 αβ n n!(α + n) t 2α+n 1 0 e -2βty y α+n (1 -y) α-1 dy = ∞ n=0 A n t 2α+n 1 F 1 (α + n + 1, 2α + n + 1, -2βt) = ∞ n=0 A n t 2α+n e -2βt 1 F 1 (α, 2α + n + 1, 2βt), (4.21)
where A n is defined in (4.18) and the last equality is due to the Kummer's transformation (4.9). From (3.1) together with (4.1) and (4.12), the sum of the remaining two terms in (3.9) is given by

µ 2 E[Z 2δ (t)] µ 1 -{E[Z δ (t)]} 2 = µ 2 λE[e -2δL ]e -(β+2δ)t t α Γ(1 + α) 1 F 1 (α, 1 + α, βt) - µ 2 1 λ 2 (E[e -δL ]) 2 e -2(β+δ)t t 2α [Γ(1 + α)] 2 [ 1 F 1 (α, 1 + α, βt)] 2 . ( 4
t→∞ ∞ n=0 A n t 2α+n e -2βt 1 F 1 (α, 2α + n + 1, 2βt) t α-1 = C (4.23)
for some constant C. Note that (4.10) implies

lim z→∞ 1 F 1 (a, b, z) Γ(b) Γ(a) e z z a-b = 1.
In order to apply such an asymptotic result, we need to check the validity of interchanging limit and infinite summation on the left-hand side of (A.57). It can be proved that the related sequence is uniformly bounded (which gives a sufficient condition). This part of the proof is not trivial and is provided in Part 1 of the 'Supplementary materials'. We can now evaluate the limit on the left-hand side of (A.57) by interchanging the order of limit and infinite summation followed by substitution of (4.18), and this gives rise to

lim t→∞ ∞ n=0 A n t 2α+n e -2βt 1 F 1 (α, 2α + n + 1, 2βt) t α-1 = ∞ n=0 A n lim t→∞ t α+n+1 e -2βt 1 F 1 (α, 2α + n + 1, 2βt) = ∞ n=0 α α + n β n n! Γ(α)Γ(α + n + 1) Γ(2α + n + 1)   1 (2β) α+n+1 Γ(2α + n + 1) Γ(α) lim t→∞ 1 F 1 (α, 2α + n + 1, 2βt) Γ(2α+n+1) Γ(α) e 2βt (2βt) -α-n-1   = α (2β) α+1 ∞ n=0 Γ(α + n) n! 2 n = α (2β) α+1 ∞ n=0 ∞ 0 y α+n-1 e -y dy n! 2 n = α (2β) α+1 ∞ 0 y α-1 ∞ n=0 (y/2) n n! e -y dy = α (2β) α+1 ∞ 0 y α-1 e -y 2 dy = α (2β) α+1 2 α Γ(α) = Γ(α + 1) 2β α+1 .
Hence, we identify C := Γ(α+1) 2β α+1 in (A.57) so that (4.21) asymptotically behaves like Ct α-1 . Combining (4.20) and (4.22) with asymptotic results of (4.21) and (4.14) yields (4.19).

Covariance of total discounted IBNR processes. We finish this section by studying the covariance of Z δ (s) and Z δ (t) and its asymptotic behaviour as t → ∞.

Proposition 4 With exponential delay in a fractional Poisson claim arrival process, the covariance of the total discounted IBNR claims Z δ (s) and Z δ (t) for 0 < s ≤ t is given by

Cov[Z δ (s), Z δ (t)] = µ 2 1 λ 2 α{E[e -δL ]} 2 [Γ(1 + α)] 2 e -(β+δ)(s+t) ∞ n=0 A n s 2α+n 1 F 1 (α, 2α+n+1, 2βs) + s α W(2βs, n, t) + µ 2 λE[e -2δL ] Γ(1 + α) e -(β+2δ)t s α 1 F 1 (α, 1+α, βs) - µ 2 1 λ 2 {E[e -δL ]} 2 [Γ(1 + α)] 2 e -(β+δ)(s+t) (st) α 1 F 1 (α, 1+α, βs) 1 F 1 (α, 1+α, βt), (4.24)
where W(2βs, n, t) := 1 B(α, α + n + 1) 1 0 e 2βsy y α-1 (t -sy) α+n dy.

(4.25)

In addition, the asymptotic behaviour when 0 < α < 1 is

Cov[Z δ (s), Z δ (t)] ∼ µ 2 1 λ(1 -α){E[e -δL ]} 2 βΓ(α) e -(β+δ

Pareto delay

This section considers a heavy-tailed distribution for random delay. In particular, it is assumed that L follows a Pareto distribution with tail

W (x) = 1 -W (x) = θ θ + x η , x ≥ 0, (4.27)
where η, θ > 0 are the parameters. In this case, E[L] = θ η-1 when η > 1. However, if 0 < η ≤ 1 then L has infinite mean, and consequently W (x) is not a directly Riemann integrable function and the Smith's Renewal Theorem is no longer applicable. In the rest of this section, it is assumed that δ = 0. The proofs of all three propositions require the next lemma.

Lemma 1 Define, for y ∈ [0, 1),

G * (y) := y 0 (1 -x) -η x α-1 dx, (4.28)
where η ≥ 1 and α ∈ (0, 1). The above integral satisfies, as y → 1 -,

G * (y) = 1 η-1 (1 -y) -η+1 + o (1 -y) -η+1 , if η > 1, -ln(1 -y) + C + o(1), if η = 1, (4.29)
for some constant C. In particular, when 1 < η ≤ 2, finer asymptotic results can be given by, as

y → 1 -, G * (y) = 1 η-1 (1 -y) -η+1 + C * + o(1), if 1 < η < 2, (1 -y) -1 -(1 -α) ln(1 -y) + C * * + o(1), if η = 2, (4.30)
for some constants C * and C * * . Proof: See Part 2 of the 'Supplementary materials'.

Expectation of total IBNR process. We have the following proposition for the mean of Z(t).

Proposition 5 With Pareto delay in a fractional Poisson claim arrival process, the mean of the total IBNR claims Z(t) when t > 0 is given by

E[Z(t)] = µ 1 λθ η Γ(α)(θ + t) η-α t θ+t 0 (1 -y) -η y α-1 dy. (4.31)
Its asymptotic behaviour is, as t → ∞,

E[Z(t)] ∼                  µ 1 λθ η Γ(1 -η) Γ(α + 1 -η) t α-η , if 0 < η < 1, µ 1 λθ Γ(α) t α-1 ln t, if η = 1, µ 1 λθ (η -1)Γ(α) t α-1 , if η > 1.
(4.32)

Proof: Substitution of (4.1) and (4.27) into (3.2) yields

E[Z(t)] = µ 1 t 0 θ θ + t -x η λx α-1 Γ(α) dx = µ 1 λ Γ(α) θ θ + t η t 0 1 - x θ + t -η
x α-1 dx, from which (4.31) follows by a change of variable from x/(θ + t) to y.

Regarding the asymptotic behaviour, we separate the analysis into three cases as follows. Case 1. 0 < η < 1: The integral in (4.31) can be expressed in terms of the incomplete Beta function (4.5), i.e.

E[Z(t)] = µ 1 λθ η Γ(α)(θ + t) η-α B α, 1 -η, t θ + t . As lim t→∞ B(α, 1 -η, t θ+t ) = B(α, 1 -η) = Γ(α)Γ(1-η) Γ(α+1-η)
, one finds the first asymptotic result in (4.32). Case 2. η = 1: Replacing y by t/(θ + t) in Lemma 1 and using the second result of (4.29) gives

t θ+t 0 (1 -y) -1 y α-1 dy = -ln θ θ + t + C + o(1) ∼ ln t as t → ∞.
Moreover, it is obvious that 1/(θ + t) 1-α ∼ t α-1 as t → ∞. By substituting the these asymptotics into (4.31), we arrive at the second result in (4.32). Case 3. η > 1: Similar to Case 2, applying the first result of (4.29) yields

t θ+t 0 (1 -y) -η y α-1 dy = 1 η -1 θ θ + t -η+1 + o θ θ + t -η+1 ∼ θ -η+1 η -1 t η-1 as t → ∞.
Combining with the fact that 1/(θ + t) η-α ∼ t α-η as t → ∞ gives the third asymptotic formula of (4.32) thanks to (4.31).

Remark 2 Consider the fractional Poisson claim arrival process with 0 < α < 1. One particular consequence of the asymptotic results (4.32) in Proposition 5 is that

lim t→∞ E[Z(t)] =      0, if 0 < α < η, µ 1 λθ η Γ(1 -η), if 0 < η = α < 1, ∞, if 0 < η < α < 1.
The above result can be intuitively understood in the following way. Since the interarrival time τ and the delay L exhibit power decay according to (4.2) and (4.27), one checks easily that τ ≤ st,∞ L (resp. L ≤ st,∞ τ ) when α > η (resp. when η > α),

Hence, if α > η then delays are too long and cannot make up for the faster arriving claims, which justifies why the expected value of IBNR claims tends to infinity as time goes by.

Variance of total IBNR process. The upcoming Proposition 6 is concerned with the asymptotic behaviour of Var[Z(t)]. We first introduce the following lemma in preparation for its proof.

Lemma 2 Let G : [0, 1] → R be a continuous function, and define the integral, for s > 0,

J(s) := s 0 (v + s -x) -γ x -ξ G s -x v + s -x dx, s > 0 (4.33)
where v > 0 is fixed. It is assumed that ξ < 1 and γ ≥ 1 so that the integral must converge. We have the asymptotic result, as s → ∞,

J(s) ∼ s -ξ v 1-γ 1 0 (1 -z) γ-2 G(z)dz, if γ > 1, s -ξ (ln s)G(1), if γ = 1. ( 4

.34)

Proof: See Part 3 of the 'Supplementary materials'.

Proposition 6 With Pareto delay in a fractional Poisson claim arrival process, the asymptotic behaviour of the variance of the total IBNR claims Z(t) when 0 < α < 1 is, as t → ∞,

Var[Z(t)] ∼                                  2µ 2 1 λ 2 θ 2η Γ(α + 1 -2η)Γ(1 -η) Γ(2α + 1 -2η)Γ(α + 1 -η) - µ 2 1 λ 2 θ 2η [Γ(1 -η)] 2 [Γ(α + 1 -η)] 2 t 2(α-η) , if 0 < η < α, µ 2 1 λ 2 θ 2α [Γ(1 -α)] 2 + µ 2 λθ α Γ(1 -α), if η = α, µ 2 λθ η Γ(1 -η) Γ(α + 1 -η) t α-η , if α < η < 1, µ 2 λθ Γ(α) t α-1 ln t, if η = 1, 2µ 2 1 λ 2 θ α+1 [Γ(α)] 2 1 0 (1 -y) 2η-α-2 G * (y)dy + µ 2 λθ (η -1)Γ(α) t α-1 , if η > 1, (4.35)
where G * (y) is given by (4.28).

Proof: Recall that the variance of Z(t) is given by (3.10). Application of (4.1), (4.27) and (4.31) yields that its first term equals

2µ 1 t 0 W (t -x)E[Z(t -x)]dm(x) = 2µ 1 t 0 θ θ + t -x η µ 1 λθ η Γ(α)(θ + t -x) η-α t-x θ+t-x 0 (1 -y) -η y α-1 dy λx α-1 Γ(α) dx = 2µ 2 1 λ 2 θ 2η [Γ(α)] 2 t 0 x α-1 (θ + t -x) 2η-α t-x θ+t-x 0 (1 -y) -η y α-1 dy dx. ( 4.36)
To analyze the asymptotic behaviour of Var[Z(t)], we mainly focus on the term (4.36). The analysis is separated into three cases as follows. Case 1. 0 < η < 1: With the help of (4.5) and (4.3), (4.36) can be rewritten as

2µ 2 1 λ 2 θ 2η [Γ(α)] 2 t 0 x α-1 (θ + t -x) 2η-α B α, 1 -η, t -x θ + t -x dx (4.37) = 2µ 2 1 λ 2 θ 2η [Γ(α)] 2 t 0 x α-1 (θ + t -x) 2η-α ( t-x θ+t-x ) α α 2 F 1 α, η, α + 1, t -x θ + t -x dx = 2µ 2 1 λ 2 θ 2η [Γ(α)] 2 t 0 x α-1 (t -x) α (θ + t -x) 2η ∞ n=0 (η) n n!(α + n) t -x θ + t -x n dx = 2µ 2 1 λ 2 θ 2η [Γ(α)] 2 ∞ n=0 (η) n n!(α + n) 1 0 t 2α+n (θ + t) 2η+n y α-1 (1 -y) α+n 1 - t θ + t y -(2η+n) dy = 2µ 2 1 λ 2 [Γ(α)] 2 θ θ + t 2η t 2α ∞ n=0 t θ + t n (η) n n!(α + n) B(α, α + n + 1) 2 F 1 2η + n, α, 2α + n + 1, t θ + t , (4.38)
where the last two lines are due to a change of variable from x/t to y together with (4.4). Note that the second and the third terms in (3.10) are equivalent to µ 2 E[Z(t)]/µ 1 and {E[Z(t)]} 2 respectively, and their asymptotics are given by (4.32) as

µ 2 E[Z(t)] µ 1 ∼ µ 2 λθ η Γ(1 -η) Γ(α + 1 -η) t α-η and {E[Z(t)]} 2 ∼ µ 2 1 λ 2 θ 2η [Γ(1 -η)] 2 [Γ(α + 1 -η)] 2 t 2(α-η) as t → ∞. (4.39)
This case can be further subdivided into two cases. Case 1a. 0 < η < α+1 2 : Now, we shall proceed to prove that (4.38) is asymptotically proportional to t 2(α-η) as t → ∞. To see this, we let

χ * n (t) := t θ + t n (η) n n!(α + n) B(α, α + n + 1) 2 F 1 2η + n, α, 2α + n + 1, t θ + t (4.40)
and look at

lim t→∞ 2µ 2 1 λ 2 [Γ(α)] 2 θ θ+t 2η t 2α ∞ n=0 χ * n (t) t 2(α-η) = 2µ 2 1 λ 2 θ 2η [Γ(α)] 2 lim t→∞ t θ + t 2η ∞ n=0 χ * n (t) = 2µ 2 1 λ 2 θ 2η [Γ(α)] 2 lim t→∞ ∞ n=0 χ * n (t) (4.41) = 2µ 2 1 λ 2 θ 2η [Γ(α)] 2 B(α, α + 1 -2η)B(α, 1 -η). ( 4
v = θ, s = t, γ = 2η -α, ξ = 1 -α, and G(x) = B(α, 1 -η, x) yields, as t → ∞, t 0 x α-1 (θ + t -x) 2η-α B α, 1 -η, t -x θ + t -x dx ∼ Ct α-1 , if α+1 2 < η < 1 t α-1 (ln t)B(α, 1 -η), if η = α+1 2 ,
where C := θ α-2η+1 1 0 (1 -z) 2η-α-2 B(α, 1 -η, z)dz. Because the assumption α+1 2 ≤ η implies α < η as 0 < α < 1, combining with the results in (4.39) we observe that the second term in (3.10) dominates the asymptotic behaviour which is proportional to t α-η . Therefore, the part of the third result of (4.35) when α+1 2 ≤ η < 1 follows. Case 2. η ≥ 1: In this case, we can re-express the integral in (4.36) in terms of (4.28) as

A * (t) := t 0 x α-1 (θ + t -x) 2η-α t-x θ+t-x 0 (1 -y) -η y α-1 dy dx = t 0 x α-1 (θ + t -x) 2η-α G * t -x θ + t -x dx. (4.43) It is shown in Appendix A.4 that A * (t) ∼ t α-1 θ α-2η+1 1 0 (1 -y) 2η-α-2 G * (y)dy as t → ∞. (4.44)
Case 2a. η = 1: According to (4.32), the asymptotics of the second and the third terms in (3.10) are

µ 2 E[Z(t)] µ 1 ∼ µ 2 λθ Γ(α) t α-1 ln t and {E[Z(t)]} 2 ∼ µ 2 1 λ 2 θ 2 [Γ(α)] 2 t 2α-2 (ln t) 2 as t → ∞.
Comparing with (4.44), it is clear that the second term is the dominant one, and one obtains the fourth result of (4.35). Case 2b. η > 1: Again, (4.32) implies that the second and the third terms in (3.10) asymptotically behave like

µ 2 E[Z(t)] µ 1 ∼ µ 2 λθ (η -1)Γ(α) t α-1 and {E[Z(t)]} 2 ∼ µ 2 1 λ 2 θ 2 (η -1) 2 [Γ(α)] 2 t 2α-2 as t → ∞.
Comparison with (4.44) reveals that the dominant terms in (3.10) are the first two, giving rise to the final result in (4.35).

Covariance of total IBNR processes. While the proof of the previous proposition requires the asymptotic results (4.29), the finer asymptotics (4.30) are needed in the proof of the upcoming proposition. Proposition 7 Let s > 0. With Pareto delay in a fractional Poisson claim arrival process, the asymptotic behaviour of the covariance of the total IBNR claims Z(t) and Z(s) when 0 < α < 1 is, as t → ∞,

Cov[Z(s), Z(t)] ∼      D 1 t -η , if 0 < η < 2 -α, D 2 t α-2 , if η = 2 -α, D 3 t α-2 , if η > 2 -α, (4.45)
where

D 1 := µ 2 λθ η s α Γ(α + 1) + µ 1 λθ η Γ(α) s 0 E[Z(s -x)]x α-1 dx, D 2 := µ 2 λθ 2-α s α Γ(α + 1) + µ 1 λθ 2-α Γ(α) s 0 E[Z(s -x)]x α-1 dx + µ 2 1 λθ Γ(α) s 0 W (s -x)xdm(x), D 3 := µ 2 1 λθ(1 -α) (η -1)Γ(α) s 0 W (s -x)xdm(x).
The expectation E[Z(s -x)] appearing in the above constants is given by (4.31). Proof: See Appendix A.5.

Numerical Examples

This section is dedicated to numerical illustrations of the theoretical results regarding some exact and asymptotic expressions of the mean, variance, covariance and correlation of IBNR processes for fractional Poisson claim arrivals with exponential and Pareto delay distributions. In all examples, it is assumed that the first two moments of the claim amount X are µ 1 = 1 and µ 2 = 4 respectively (so that Var[X] = 3 and the coefficient of variation of √ 3 indicates claim with high variability). In addition, the parameter λ of the fractional Poisson process is assumed to be λ = 1.5, and except for Example 5 we assume a force of interest of δ = 0.

Example 2 (Impact of parameters of delay distribution) In this example, the underlying claim arrival process is assumed to have fractional Poisson parameter α = 0.6, and both E(β) and Pareto(θ, η) delays will be considered under various choices of parameters. We first calculate the asymptotic mean and variance of the IBNR claims Z(t) by plugging in a large t (where t = 100, 000) into our asymptotic expressions. The results are listed in Table 1 andTable 2. The corresponding exact values, whenever they can be calculated explicitly using our theoretical results, are quoted in parenthesis. The asymptotic and the exact values are very close in each case. Table 1 shows that E[Z(t)] decreases as β increases for exponential delays. Intuitively, this is due to the fact that E(β 1 ) ≤ st E(β 2 ) for β 1 > β 2 (where Y 1 ≤ st Y 2 means that usual stochastic ordering P(Y 1 > x) ≤ P(Y 2 > x) for all x ≥ 0), as stochastically smaller delays imply that incurred claims are reported quicker and there are less IBNR claims. Similarly, when one switches to Pareto delay, the monotonicity of E[Z(t)] in θ and η can be explained by the stochastic orderings Pareto(θ, η 1 ) ≤ st Pareto(θ, η 2 ) for η 1 > η 2 and Pareto(θ 1 , η) ≤ st Pareto(θ 2 , η) for θ 1 < θ 2 . Note that the values of Var[Z(t)] in Table 2 show the same pattern as those in Table 1. Moreover, the values of E[Z(t)] and Var[Z(t)] in Table 1 are of small magnitude for exponential delays, and this is consistent with the implication of the summary in (2.1) that E[Z(t)] and Var[Z(t)] tend to zero as t tends to infinity. However, for Pareto delays, the magnitude of E[Z(t)] and Var[Z(t)] is large when η = 0.2, 0.4 but is small when η = 1.0, 1.2, 1.4. Such an observation is consistent with (2.2) which indicates that both E[Z(t)] and Var[Z(t)] tend to infinity or zero depending on whether η < α or η > α (recall that α is fixed at 0.6). See also Remark 2.

Next, the asymptotic covariance and correlation between Z(s) and Z(t) are provided in Table 3 and Table 4 respectively. In each table, we consider two scenarios with different initial time s but the same time frame t -s. Specifically, in Case 1 we assume s = 10, 000 and t = 100, 000 while in Case 2 we let s = 20, 000 and t = 110, 000. Although the asymptotic covariances in Table 3 show similar pattern to those in Table 1 and Table 2 concerning means and variances, it is important to note that the covariance should always tend to zero as t tends to infinity according to (2.1) and (2.2). The high values of the covariance for Pareto delays when η = 0.2, 0.4 are attributed to the fact that the proportionality constant D 1 is large in these cases (see (4.45)). In calculating the asymptotic correlation in Table 4, we are dividing the results in Table 3 by the square root of the product of the asymptotic variances Var[Z(s)] and Var[Z(t)] because both s and t are large. It is observed that correlation decreases when β increases for E(β) delay or when η increases for Pareto(θ, η) delay. However, the behaviour of the correlation in response to a change in θ for Pareto(θ, η) delay depends on the fixed value of η. In particular, the correlation decreases in θ when η = 0.2, 0.4 but increases in θ when η = 1.0, 1.2, 1.4. For both exponential and Pareto delays, the correlation increases as the initial time s increases (while keeping the time frame t -s fixed).

Exponential

Pareto Table 3: Asymptotic covariance between Z(s) and Z(t) Example 3 (Impact of fractional Poisson parameter α) In this example, the effect of α (where 0.1 < α < 0.9) on the asymptotics of the IBNR claims is examined. We focus on the case of Pareto delays with t = 100, 000 (and additionally s = 10, 000 for covariance and correlation). Moreover, we let θ = 1.0 and consider three scenarios where η = 0.4, 1.0, 1.4. The results are given in Figures 1234. For each fixed η, the asymptotic mean in Figure 1 is increasing in α. Intuitively, as α increases, the interarrival times of the claims becomes less heavy-tailed (see (4.2)) so that claims arrive more frequently, thereby increasing the IBNR claims. It is also noted that the asymptotic mean is continuous in α, as evident from each of the three pieces in (4.32). The asymptotic variance in Figure 2 and the asymptotic covariance in Figure 3 show similar trend at a first glance. However, we remark that for the case η = 0.4 (which corresponds to the first subfigure of Figure 2) there is indeed a discontinuity at α = η, which can hardly be seen due to the scale of the plot. Such a curve is plotted using the first three pieces of (4.35). In contrast, the second and third subfigures of Figure 2 for the cases η = 1.0 and η = 1.4 are continuous in α and are plotted using the fourth and fifth pieces of (4.35). Similarly, the first two subfigures of Figure 3 concerning the covariance are continuous in α and are both plotted using the first piece of (4.45), but the third subfigure is discontinuous at α = 0.6 and is plotted using all three pieces of (4.45). The aforementioned discontinuities explain those observed in the first and third subfigures of Figure 4. Figure 4: Asymptotic correlation between Z(s) and Z(t) against α 100, 000 to 200, 000. The results for α = 0.3, 0.5, 0.8 are plotted in Figure 5 for both E(β = 1) and Pareto(θ = 1, η = 1) delays. The correlation always decreases in t because an IBNR claim at time s is more likely to be reported by time t as t increases. For exponential delays (where the IBNR process is always SRD), the correlation increases as the fractional Poisson parameter α increases. This can be explained by (2.3) where the correlation is always asymptotically proportional to t -3-α 2 which is increasing in α. Note that the proportionality constant also depends on α, but one expects that the effect of a change in α on t -3-α 2 dominates for large t. On the other hand, the correlation under Pareto delays decreases in α when α increases from 0.3 through 0.5 to 0.8 in this example. However, one should note from Example 3 that in general the asymptotic correlation is not necessarily monotone in α (see Figure 4).

β E[Z(t)] η E[Z(t)] θ = 0.5 θ = 1.0 θ = 2.0 0.
Example 5 (Impact of force of interest δ) This final example briefly illustrates the effect of δ on the exact values of mean and variance of the IBNR claims for exponential delays with parameter β = 1.0. We let α = 0.6 and consider four time values t = 20, 50, 100, 200 in Figure 6. It can be seen that both E[Z δ (t)] and Var[Z δ (t)] quickly approaches zero as δ increases due to heavier discounting. Such an effect is more pronounced for a larger value of t, which is expected because of the presence of e -δt and e -2δt in the asymptotic formulas for the mean and variance respectively (see (2.1)).

A.2 Proof of Theorem 2

First, from (3.3) with the application of binomial expansion it follows that, for 0 < s ≤ t and n, m ∈ N,

Z n δ (s)Z m δ (t) = e -(n+m)δT 1 I {T 1 ≤s} I {T 1 ≤t} e -δL 1 I {T 1 +L 1 >s} X 1 + Z 2,δ (s) n e -δL 1 I {T 1 +L 1 >t} X 1 + Z 2,δ (t) m = e -(n+m)δT 1 I {T 1 ≤s} Z n 2,δ (s) + n-1 i=0 n i e -(n-i)δL 1 I {T 1 +L 1 >s} X n-i 1 Z i 2,δ (s) × Z m 2,δ (t) + m-1 j=0 m j e -(m-j)δL 1 I {T 1 +L 1 >t} X m-j 1 Z j 2,δ (t) = e -(n+m)δT 1 I {T 1 ≤s} Z n 2,δ (s)Z m 2,δ (t) + I {T 1 +L 1 >s} n-1 i=0 n i e -(n-i)δL 1 X n-i 1 Z i 2,δ (s)Z m 2,δ (t)
+ I {T 1 +L 1 >t} m-1 j=0 n i=0 n i m j e -(n+m-i-j)δL 1 X n+m-i-j 1 Z i 2,δ (s)Z j 2,δ (t) .
Taking expectations on the both sides of above equation via conditioning on

T 1 = x results in E[Z n δ (s)Z m δ (t)] = s 0 e -(n+m)δx E[Z n δ (s-x)Z m δ (t-x)] + n-1 i=0 µ n-i n i E[Z i δ (s-x)Z m δ (t-x)] ∞ s-x e -(n-i)δy dW (y) + m-1 j=0 n i=0 µ n+m-i-j n i m j E[Z i δ (s-x)Z j δ (t-x)] ∞ t-x e -(n+m-i-j)δy dW (y) dF (x) = s 0 e -(n+m)δx E[Z n δ (s-x)Z m δ (t-x)] + n-1 i=0 µ n-i n i E[Z i δ (s-x)Z m δ (t-x)]e -(n-i)δ(s-x) T (n-i)δ w(s-x) + m-1 j=0 n i=0 µ n+m-i-j n i m j E[Z i δ (s-x)Z j δ (t-x)]e -(n+m-i-j)δ(t-x) T (n+m-i-j)δ w(t-x) dF (x), (A.3)
where the notation of D-H operator has been used. Similar to the arguments used in the proof of Theorem 1, we proceed by defining g n,m,h (s

) := E[Z n δ (s)Z m δ (s + h)] and g n,m,h (s)
:= e (n+m)δs g n,m,h (s) for s > 0 and fixed h ≥ 0. Replacing t by s + h in (A.3) and multiplying both sides by e (n+m)δs , we arrive at the renewal equation

g n,m,h (s) = g n,m,h * F (s) + a n,m,h * F (s) (A.4)
for s > 0, where a n,m,h (s) := e (n+m)δs a n,m,h (s) with

a n,m,h (s) := n-1 i=0 µ n-i n i E[Z i δ (s)Z m δ (s+h)]e -(n-i)δs T (n-i)δ w(s) + m-1 j=0 n i=0 µ n+m-i-j n i m j E[Z i δ (s)Z j δ (s+h)]e -(n+m-i-j)δ(s+h) T (n+m-i-j)δ w(s+h).
Analogous to (A.2), the solution of (A.4) is

g n,m,h (s) = s 0 a n,m,h (s -x)dm(x),
which is identical to (3.5) upon dividing by e (n+m)δs and then replacing h with t -s.

A.3 Proof of Proposition 4

1. Proving the exact formula (4.24). We start by proving (4.24). First, plugging (4.12), T δ w(x) = β β+δ e -βx = E[e -δL ]e -βx and (4.1) into the first term on the right-hand side of (3.6) yields

µ 1 e -δs s 0 e -δx E[Z δ (t -x)]T δ w(s -x)dm(x) = µ 2 1 λ 2 α{E[e -δL ]} 2 [Γ(1 + α)] 2 e -(β+δ)(s+t) s 0 e 2βx x α-1 (t -x) α 1 F 1 (α, 1 + α, β(t -x))dx. (A.5)
Again, with the application of (4.7) followed by changing a variable x/s to y, we find that the integral in (A.5) can be expressed as

∞ n=0 αβ n n!(α + n) s α 1 0 e 2βsy y α-1 (t -sy) α+n dy = ∞ n=0 A n s α W(2βs, n, t), (A.6)
where A n and W(2βs, n, t) are given by (4.18) and (4.25) respectively. Similarly, the third term of (3.6) can be represented as

µ 1 e -δt s 0 e -δx E[Z δ (s -x)]T δ w(t -x)dm(x) = µ 2 1 λ 2 α{E[e -δL ]} 2 [Γ(1 + α)] 2 e -(β+δ)(s+t) s 0 e 2βx x α-1 (s -x) α 1 F 1 (α, 1 + α, β(s -x))dx = µ 2 1 λ 2 α{E[e -δL ]} 2 [Γ(1 + α)] 2 e -(β+δ)(s+t) ∞ n=0 A n s 2α+n 1 F 1 (α, 2α + n + 1, 2βs), (A.7)
where the last equality is due to (4.21). Lastly, using (3.1) and the fact that T δ w(t-x) = e -β(t-s) T δ w(sx), we observe that the second term in (3.6) is simply

µ 2 e -2δt s 0 T 2δ w(t -x)dm(x) = µ 2 µ 1 e -(β+2δ)(t-s) E[Z 2δ (s)], (A.8)
which corresponds to the second term on the right-hand side of (4.24) thanks to (4.12). It is also clear that the last terms in (3.6) and (4.24) are identical again because of (4.12). Therefore, combining this with (A.5)-(A.7) yields (4.24).

2. Proving the asymptotic formula (4.26). We now proceed to prove the asymptotic result (4.26). From (A.7) and (A.8), it is noted that the third and the second terms of (3.6) are O(e -(β+δ)t ) and O(e -(β+2δ)t ) respectively. In what follows, it will be shown that these two terms are asymptotically dominated by the first and the last terms of (3.6). Using (3.1), (4.15) and T δ w(s-x) = E[e -δL ]e -β(s-x) , the sum of the first and the last terms of (3.6) can be expressed as

µ 1 e -δs s 0 e -δx E[Z δ (t -x)]T δ w(s -x)dm(x) -E[Z δ (s)]E[Z δ (t)] = µ 1 e -δs s 0 e -δx E[Z δ (t -x)] -E[Z δ (t)] T δ w(s -x)dm(x) = µ 1 {E[e -δL ]} 2 e -(β+δ)s e -δt s 0 e βx E[Z(t -x)] -E[Z(t)] dm(x). (A.9)
In view of the integrand above, one needs to first study the asymptotics of E[Z(t -x)] -E[Z(t)] as t → ∞ for all x ∈ [0, s] (where s ∈ (0, t] is fixed). By applying (4.16), we rewrite

E[Z(t -x)] -E[Z(t)] as, for x ∈ [0, s], E[Z(t -x)] -E[Z(t)] = µ 1 λ Γ(α) [t α-2 Ξ(t, x) -Π(t, x)], (A.10)
where

Ξ(t, x) := 1 0 [e βx(1-y) -1]t 2 e -βt(1-y) y α-1 dy, (A.11) Π(t, x) := [t α -(t -x) α ] 1 0
e -β(t-x)(1-y) y α-1 dy, (A.12) will be analyzed separately.

Term Ξ(t, x). Performing the change of variable v = t(1 -y) to the term Ξ(t, x) in (A.11) yields

Ξ(t, x) = t 0 (e βxv/t -1)te -βv 1 - v t α-1 dv = t/2 0 + t t/2 := Ξ 1 (t, x) + Ξ 2 (t, x) (A.13)
with the obvious definitions of Ξ 1 (t, x) and Ξ 2 (t, x). We first look at Ξ 1 (t, x). Note that (e βxv/t -1)t = ∞ j=1 (βxv) j /(j!t j-1 ) is decreasing in t > 0. Therefore, for t large enough such that t ≥ 2s, one has that, for x ∈ [0, s],

I {0≤v≤t/2} (e βxv/t -1)te -βv 1 - v t α-1 ≤ I {0≤v≤t/2} (e βxv/2s -1)2se -βv 1 - v t α-1 ≤ I {0≤v≤t/2} (e βv/2 -1)2se -βv 1 2 α-1 ≤ s 2 α-2 (e -βv/2 -e -βv ), (A.14)
which is integrable with respect to v on [0, ∞). For all v ≥ 0, (e βxv/t -1)t and (1 -v t ) α-1 tend to βxv and 1 respectively as t → ∞. Then, by dominated convergence we arrive at

lim t→∞ Ξ 1 (t, x) = ∞ 0 lim t→∞ I {0≤v≤t/2} (e βxv/t -1)te -βv 1 - v t α-1 dv = ∞ 0 βxve -βv dv = x β for all x ∈ [0, s]. (A.15) As for Ξ 2 (t, x), changing variable z = v/t leads to Ξ 2 (t, x) = 1 1/2 (e βxz -1)t 2 e -βtz (1 -z) α-1 dz,
which can be upper bounded as

|Ξ 2 (t, x)| ≤ (e βs -1)t 2 e -βt/2 1 1/2
(1 -z) α-1 dz for all x ∈ [0, s] and t ≥ s, (A. [START_REF] Keilson | M/G/∞ with batch arrivals[END_REF] which tends to 0 as t → ∞. Combining with (A.15), we conclude that the limit of (A.

|Π(t, x)| ≤ C * x t (C Π * 1 + C Π * 2 )t α-1 ≤ C Π t α-2
for all x ∈ [0, s] and t ≥ max{s, t * }, (A. [START_REF] Léveillé | Joint moments of discounted compound renewal sums[END_REF] where

C Π := C * s(C Π * 1 + C Π * 2 ).
The result. Plugging (A.17) and (A.22) into (A.10) yields the asymptotic equivalent, for all x ∈ [0, s],

E[Z(t -x)] -E[Z(t)] ∼ µ 1 λ Γ(α) 1 -α β xt α-2 as t → ∞. (A.24)
We shall now establish an equivalent for the integral term in (A.9) and consider

s 0 e βx E[Z(t -x)] -E[Z(t)] dm(x) t α-2 = s 0 e βx E[Z(t -x)] -E[Z(t)] t α-2 dm(x). (A.25)
Utilizing the definition (A.10) together with the bounds (A.18) and (A.23), it is clear that

E[Z(t -x)] -E[Z(t)] t α-2 ≤ µ 1 λ Γ(α) (C Ξ + C Π )
for all x ∈ [0, s] and t ≥ max{2s, t * }.

Taking limit as t → ∞ on both sides of (A.25) along with the use of dominated convergence and (A. [START_REF] Léveillé | A compound renewal model for medical malpractice insurance[END_REF], it is easy to see that

lim t→∞ s 0 e βx E[Z(t -x)] -E[Z(t)] dm(x) t α-2 = µ 1 λ Γ(α) 1 -α β s 0 e βx xdm(x). (A.26)
It now remains to conclude how the asymptotic formula (4.26) is obtained. Recall that the second and the third terms of (3.6) are respectively O(e -(β+2δ)t ) and O(e -(β+δ)t ). These are negligible compared to the sum of the first and the last terms, namely (A.9), which is O(e -δt t α-2 ) because of (A.26). Applying (A.9) and (A.26) again yields the result (4.26).

A.4 Proof of auxiliary results for Proposition 6

C * n = (η) n n!(α + n) Γ(α)Γ(α + 1 -2η) Γ(2α + 1 -2η) = (η) n n!(α + n) B(α, α + 1 -2η).
Then, via application of (4.3) and (4.6) (requiring α + 1 > α + η, i.e. η < 1 which is satisfied), we can sum C * n 's as

∞ n=0 C * n = B(α, α + 1 -2η) ∞ n=0 (η) n n!(α + n) = B(α, α + 1 -2η) 1 α 2 F 1 (α, η, α + 1, 1) = B(α, α + 1 -2η) Γ(α)Γ(1 -η) Γ(α + 1 -η) = B(α, α + 1 -2η)B(α, 1 -η) < ∞. (A.28)
Hence, uniform convergence is proved. Then, further using the definitions of χ * n (t) and C * n in (4.40) and (A.27), we can now evaluate the limit on the right-hand side of (4.41) as

lim t→∞ ∞ n=0 χ * n (t) = ∞ n=0 lim t→∞ χ * n (t) = ∞ n=0 (η) n n!(α + n) B(α, α + n + 1) lim t→∞ t θ + t n 2 F 1 2η + n, α, 2α + n + 1, t θ + t = ∞ n=0 (η) n n!(α + n) B(α, α + n + 1) 2 F 1 (2η + n, α, 2α + n + 1, 1) = ∞ n=0 C * n = B(α, α + 1 -2η)B(α, 1 -η), (A.
A * (t) = θ α-2η+1 t α-1 1 0 (1 -y) α-1 1 - t θ + t y 2η-2α-1 G * t θ + t y dy t θ + t . (A.30)
Based on the fact that t θ+t → 1 -as t → ∞, we shall establish the limit of A * (t) t α-1 as t → ∞. It is first noted that for all y ∈ [0, 1), G * ( t θ+t y) increases towards G * (y) as t → ∞. In order to determine the limit of A * (t) t α-1 , we shall use the asymptotic result of G * (y) as y → 1 -in Lemma 1. In what follows, we focus on the case η > 1 as the argument will be similar when η = 1. We shall need the following lemma.

Lemma 3 Given θ > 0. For every z ∈ (0, 1) one has the limiting result

lim t→∞ 1 z (1 -y) α-1 1 - t θ + t y η-2α dy = 1 z (1 -y) η-α-1 dy,
which is a convergent integral for η > 1 and α ∈ (0, 1). Proof: We distinguish two cases according to the sign of η -2α. Case 1. η ≤ 2α: Because 1 -t θ+t y η-2α increases towards (1 -y) η-2α as t → ∞, the result follows by monotone convergence. Case 2. η > 2α: For all t ≥ 0 and y ∈ (z, 1), it is clear that

(1 -y) α-1 1 - t θ + t y η-2α ≤ (1 -y) α-1 .
Since 1 z (1 -y) α-1 dy is finite, application of the dominated convergence theorem yields the result. From (4.29), one has that there exists x 0 ∈ (0, 1) such that

|G * (x)| ≤ 3 2 1 η -1 |1 -x| -η+1 for all x ∈ (x 0 , 1).
This implies that we can first fix a z ∈ (x 0 , 1), and then there exists t 0 large enough such that

G * t θ + t y ≤ 3 2 1 η -1 1 - t θ + t y -η+1
for all y ∈ (z, 1) and t ≥ t 0 . (A.31)

For example, one may pick t 0 such that t 0 θ+t 0 z = x 0 , i.e. t 0 = θx 0 z-x 0 . Next, due to (A.30), we turn our attention the asymptotic behaviour of

A * (t) t α-1 1 θ α-2η+1 θ + t t = 1 0 (1-y) α-1 1- t θ + t y 2η-2α-1 G * t θ + t y dy = z 0 + 1 z := J 1 (z, t)+J 2 (z, t), (A.32)
where J 1 (z, t) and J 2 (z, t) have obvious definitions, and z ∈ (x 0 , 1) is fixed. Starting with the term J 1 (z, t), it is noted that |[1 -t θ+t y] 2η-2α-1 G * ( t θ+t y)| is upper bounded by max{1, (1 -z) 2η-2α-1 }G * (z) for all t ≥ 0 and y ∈ [0, z] and the integral z 0 (1 -y) α-1 dy is finite. Application of the dominated convergence theorem gives rise to

lim t→∞ J 1 (z, t) = z 0 (1 -y) α-1 (1 -y) 2η-2α-1 G * (y)dy = z 0 (1 -y) 2η-α-2 G * (y)dy. (A.33)
As for J 2 (z, t), using (A.31) one has the upper bound

|J 2 (z, t)| ≤ 3 2 1 η -1 1 z (1 -y) α-1 1 - t θ + t y 2η-2α-1 1 - t θ + t y -η+1 dy, = 3 2 1 η -1 1 z (1 -y) α-1 1 - t θ + t y η-2α
dy for all t ≥ t 0 .

Thus, it follows from Lemma 3 that

lim sup t→∞ |J 2 (z, t)| ≤ 3 2 1 η -1 1 z (1 -y) η-α-1 dy, which is equivalent to lim sup t→∞ A * (t) t α-1 1 θ α-2η+1 θ + t t -J 1 (z, t) ≤ 3 2 1 η -1 1 z (1 -y) η-α-1 dy (A.34)
thanks to (A.32). Let us now observe that, by the triangular inequality,

A * (t) t α-1 1 θ α-2η+1 θ + t t - 1 0 (1 -y) 2η-α-2 G * (y)dy ≤ A * (t) t α-1 1 θ α-2η+1 θ + t t -J 1 (z, t) + J 1 (z, t) - 1 0 (1 -y) 2η-α-2 G * (y)dy ,
and therefore, with the help of (A.34) and (A.33) we arrive at

lim sup t→∞ A * (t) t α-1 1 θ α-2η+1 θ + t t - 1 0 (1 -y) 2η-α-2 G * (y)dy ≤ 3 2 1 η -1 1 z (1 -y) η-α-1 dy + 1 z (1 -y) 2η-α-2 G * (y)dy. (A.35)
While the first integral on the right-hand side is clearly finite, the second integral is also convergent thanks to the estimate (4.29). Since z ∈ (x 0 , 1) is arbitrary, we can let z → 1 -in the above inequality, ending with the right-hand side converging to 0. We thus obtain (4.44) for η > 1.

Finally, we remark that the procedure to derive (4.44) in the case η = 1 is similar. The splitting of (A.32) into J 1 (z, t) and J 2 (z, t) can still be adopted but the major difference is that (A.31) needs to be modified using the second case of (4.29). The limit (A.33) is still valid, and one can obtain a similar upper bound to (A.35) which converges to 0. The details are omitted.

A.5 Proof of Proposition 7

Recall that s is fixed such that 0 < s ≤ t. As in the proof of Proposition 4 for the case of exponential delays, to analyze the asymptotic behaviour of Cov[Z(s), Z(t)] as t → ∞ we shall look at the terms appearing in the general expression (3.7). We can write Cov[Z(s), Z(t)] = J * 1 (t) + J * 2 (t) + J * 3 (t) for 0 < s ≤ t, where

J * 1 (t) := µ 2 s 0 W (t -x)dm(x), (A.36) J * 2 (t) := µ 1 s 0 E[Z(s -x)]W (t -x)dm(x), (A.37)
where G * (•) is defined in (4.28) (with the definition extended to allow for 0 < η < 1 as well). Hence, one has E[Z(t)] = µ 1 λ Γ(α) t α-η ϕ θ t thanks to (4.31). Defining

ψ(t, x) := ϕ θ t -x -ϕ θ t (A.43) for x ∈ [0, s], we can write E[Z(t -x)] -E[Z(t)] = µ 1 λ Γ(α) (t -x) α-η ψ(t, x) + [(t -x) α-η -t α-η ]ϕ θ t ,
and then J * 3 (t) in (A.41) can be decomposed as J * 3 (t) = J * 31 (t) + J * 32 (t) where Note that the limiting behaviour of t η-2 as t → ∞ depends on the value of η. It will be seen that L'Hôpital's rule can be applied when η = 2 where the limit on the left-hand side of (A.46) is of the indeterminate form 0/0 or ∞/∞. This involves differentiating ψ(x, t) in (A.43) with respect to t, yielding

J * 31 (t) := µ 2 1 λ Γ(α) s 0 W (s -x)(t -x) α-η ψ(t, x)dm(x), (A.44) J * 32 (t) := µ 2 1 λ Γ(α) ϕ θ t s 0 W (s -x)[(t -x) α-η -t α-η ]dm(x
∂ ∂t ψ(t, x) = - θ (t -x) 2 ϕ θ t -x + θ t 2 ϕ θ t = - θ (t -x) 2 + θ t 2 ϕ θ t -x + θ t 2 ϕ θ t -ϕ θ t -x ,
(A.47) where the derivative of ϕ(•) in (A.42) is given by

ϕ (v) = -θ η v -η (v + 1) -1 + θ η (α -η)(v + 1) α-η-1 G * 1 v + 1 .
Four cases are considered as follows.

Case 1. η ∈ (1, 2) ∪ (2, ∞): Utilizing (4.29) in Lemma 1, it is not difficult to see that G * 1 v+1 = 1 η-1 v -η+1 + o(v -η+1
) as v → 0 + , and consequently one has the expansion

ϕ (v) = -θ η v -η + θ η v -η [1 -(v + 1) -1 ] + θ η (α -η)(v + 1) α-η-1 1 η -1 v -η+1 + o(v -η+1 ) = -θ η v -η + θ η v -η+1 + θ η α -η η -1 v -η+1 + o(v -η+1 ) = -θ η v -η -θ η 1 -α η -1 v -η+1 + o(v -η+1 ) as v → 0 + . (A.48)
The first term on the right-hand side of (A.47) thus satisfies the asymptotic formula

- θ (t -x) 2 + θ t 2 ϕ θ t -x ∼ -- θ (t -x) 2 + θ t 2 θ η θ t -η = θt η-2 (1 -x/t) -2 -1 ∼ 2θxt η-3 as t → ∞. Since 1 t 2 [t η -(t -x) η ] ∼ ηxt η-3 and 1 t 2 t η-1 -(t -x) η-1 ∼ (η -1)xt η-4 = o(t η-3
) as t → ∞, the second term on the right-hand side of (A.47) can again be estimated using (A.48) via

θ t 2 ϕ θ t -ϕ θ t -x = - θ t 2 [t η -(t -x) η ] - θ 2 t 2 1 -α η -1 t η-1 -(t -x) η-1 + o(t η-3 ) = -ηθxt η-3 + o(t η-3 ) as t → ∞.
Plugging the above two estimates into (A.47) yields, for all x ∈ [0, s],

∂ ∂t ψ(t, x) ∼ -(η -2)θxt η-3 as t → ∞. (A.
ϕ θ t = θ η θ t + 1 α-η 1 η -1 θ/t θ/t + 1 -η+1 + C * + o(1) as t → ∞
for some constant C * . Further expanding θ t + 1 a = 1 + a θ t + o 1 t at a = α -η and a = η -1, we obtain In view of the definition (A.43), taking difference of the above equation with (A.51) results in ψ(t, x) = o(1) as t → ∞, i.e. the left-hand side of (A.46) is in the form 0/0. Consequently, the L'Hôpital's rule can be applied as in (A.50), and further use of (A.49) leads to (A.46). Case 2. η ∈ (0, 1): In this case, it is clear that G * 1 v+1 = B α, 1 -η, 1 v+1 = B(α, 1 -η) + o(1) as v → 0 + . Therefore, it is clear from (A.42) and (A.43) that ψ(t, x) → 0 as t → ∞, implying that the left-hand side of (A.46) is 0/0. Note that (A.48) is now replaced by ϕ 1) as v → 0 + where C * is a constant. Note that the first term is the dominant one, and the procedure leading to (A.49) is still valid. It remains to show that ψ(t, x) → 0 as t → ∞ so that L'Hôpital's rule in (A.50) can be applied and (A.46) holds true. To see this, we substitute the above estimation of G * 1 v+1 into (A.42) (under v = θ/t) so that, similar to (A.51),

(v) = -θ η v -η + θ η (α -η)B(α, 1 -η) + o(
(v) = -θv -1 + θ(1 -α) ln v + C * + o(
ϕ θ t = θ θ t + 1 α-1 -ln θ t + C + o(1) = θ[ln t + C * + o(1)] as t → ∞
for some constant C * . One also has the same equation but with t replaced by t -x, and therefore taking difference results in

ψ(t, x) = θ ln t -x t + o(1) = o(1) as t → ∞.
Case 4. η = 2: This final case is different from the previous ones in the sense that t η-2 in the denominator on the left-hand side of (A.46) equals to one. We shall proceed directly to calculate the limit lim t→∞ ψ(t, x) as follows. First, application of the finer asymptotic result (4.30) to (A.42) with v = θ/t yields, for some constant C * * , Having proved (A.46) for η > 0, we are ready to conclude on J * 31 (t). For all x ∈ [0, s] and t ≥ 2s, two observations are made as follows. First, it is clear that |(1 -x/t) α-η | ≤ max(1, 1/2 α-η ). Second, the fact that ϕ(v) is decreasing in v implies |ψ(t, x)/t η-2 | ≤ |ψ(t, s)/t η-2 |, which is upper bounded as it is continuous and convergent to a finite limit as t → ∞ thanks to (A.46). Using the definition (A.44), we can thus apply the dominated convergence theorem to obtain

ϕ θ t = θ 2 θ t + 1 α-2 θ/t θ/t + 1 -1 -(1 -α) ln θ/t θ/t + 1 + C * * + o(1) = θ 2 θ t -1 θ t + 1 α-1 -(1 -α) θ t + 1 α-2 ln θ t + (1 -α) θ t + 1 α-2 ln θ t + 1 + C * * θ t + 1 α-2 + o(1)
J * 31 (t) t α-2 = µ 2 1 λ Γ(α) s 0 W (s -x) 1 - x t
α-η ψ(t, x) t η-2 dm(x) -→ -

µ 2 1 λθ Γ(α) s 0 W (s -x)xdm(x) as t → ∞.
(A.53) Term J * 32 (t) of J * 3 (t). Note that J * 32 (t) defined in (A.45) is identical to zero if η = α, and therefore it is sufficient to focus on the case η = α. For all x ∈ [0, s] and t ≥ s, one obtains easily that 1 t α-η-1 [(t -x) α-η -t α-η ] → (η -α)x as t → ∞. Moreover, one has the inequality An asymptotic formula of J * 32 (t) will follow once we find one for ϕ θ t . Recall that (A.42) implies ϕ θ t ∼ θ η G * 1 θ/t+1 as t → ∞, where G * 1 θ/t+1 is asymptotically equivalent to 1 η-1 θ t -η+1 (resp. ln t) when η > 1 (resp. η = 1) according to (4.29) in Lemma 1, and tends to B(α, 1-η) when 0 < η < 1. Consolidating these results with (A.54), we arrive at the asymptotic behaviour for J * 32 (t) given by, as t → ∞, (κ-1) 2 are non-negative constants. A crucial point concerning these two constants is that they only depend on κ (that is to be fixed later on) and are always finite under κ > 1. Incorporating the inequalities (A.75) and (A.76) into (A.74) leads to Then, with the help of (A.78), χ n,21 (t) in (A.71) can be upper bounded as, for κ > 1,

| 1 t α-η-1 [(t -x) α-η - t α-η ]| ≤ | 1 t α-η-1 [(t -s) α-η -t α-η ]|
J * 32 (t) ∼                    µ 2 1 λθ η (η -α) Γ(α) B(α, 1 -η) s 0 W (s -x)xdm(x) t α-η-1 , if 0 < η < 1, µ 2

. 10 )

 10 See Equation (13.1.4) on p.504 of Abramowitz and Stegun (1972).

)s s 0 e

 0 βx xdm(x) e -δt t α-2 as t → ∞. (4.26) Proof: See Appendix A.3.

4 Figure 1 : 4 Figure 2 : 4 Figure 3 :

 414243 Figure 1: Asymptotic mean of Z(t) against α

Figure 5 :Figure 6 :

 56 Figure 5: Asymptotic correlation over time t for exponential (left) and Pareto (right) delays t = 20 t = 50 t = 100 t = 200

  [START_REF] Mikosch | Non-Life Insurance Mathematics: An Introduction with the Poisson Process[END_REF] from which (4.42) follows. Note that 2 F 1 (a, b, c, z) is continuous at z = 1 when Re(c -a -b) > 0 (i.e. η < α+1 2 in our case), and (A.28) has been applied. 2. Proving (4.44) when η ≥ 1. Via a change of variable with y = t+θ t t-x θ+t-x in (4.43), i.e. x = (θ + t) 1 -θ θ+t(1-y) , so that x ∈ [0, t] implies y ∈ [0, 1], and dx = -(θ + t) θt [θ+t(1-y)] 2 dy, we obtain

  1), and one can check that the arguments leading to (A.49) and (A.46) still hold. Case 3. η = 1: According to (4.29) in Lemma 1, one has G * 1 v+1 = -ln v + C + o(1) as v → 0 + for some constant C. With some simple algebra it can be shown that (A.48) is replaced by ϕ

0 W

 0 which is an upper bounded quantity. The dominated convergence theorem thus yields from (A.45) the equivalent (s -x)xdm(x) t α-η-1 as t → ∞.(A.54)

2 κ- 1 [B 1 (n + 1 ) + C 1 ] n + 1 ≤ M * 1 ( 1 .β 1 2β 1 2

 2111111111 2βt) α for 2βt > κ(n + 1) with κ > 1,(A.77) whereM * 1 := 2κα κ-1 exp 2(B 1 +C 1 )κ-With (A.73) and (A.77), we obtain the upper bound|U(α, 2α + n + 1, 2βt)| ≤ 1 (2βt) α + |ε n,21 (2βt)| ≤ M * * 1 (2βt) α for 2βt > κ(n + 1) with κ > 1, (A.78) where M * * 1 := M * 1 + 1. Next, noting that e -x x n+1 (as a function of x) is decreasing on [n + 1, ∞), we have e -2βt (2βt) n+1 ≤ e -κ(n+1) [κ(n + 1)] n+1 for 2βt > κ(n + 1) with κ > 1.Consequently,t n+1 e -2βt β n = 1 n+1 e -2βt (2βt) n+1 ≤ 1 n+1 e -κ(n+1) [κ(n + 1)] n+1 for 2βt > κ(n + 1) with κ > 1.

  |χ n,21 (t)| ≤ I {2βt>κ(n+1)} t α 1 β 1 2 n+1 e -κ(n+1) [κ(n + 1)] n+1 α α + n as n → ∞.As a result, application of ratio test gives, for κ > 1last inequality follows from the fact that xe -x+1 /2 (as a function of x) is decreasing on [1, ∞). Hence,∞ n=0 V n,21 < ∞. (A.80) Term χ n,22 (t) of χ n,2(t). Lastly, we consider the term χ n,22 (t) in (A.72). Applying again Olver et al. (2006, Chapter 13.7(ii)) with n = 1, a = α + n + 1 and b = 2α + n + 1 implies U(α + n + 1, 2α + n + 1, e ±iπ (2βt)) = 1 [e ±iπ (2βt)] α+n+1 + ε n,22 (e ±iπ (2βt)).(A.81)The term e ±iπ (2βt) in the above Tricomi's function this time tends to -∞ regardless of whether one chooses to multiply 2βt by e iπ or e -iπ . Even though the value of U(α + n + 1, 2α + n + 1, e ±iπ (2βt)) is different depending on whether we choose the argument e iπ (2βt) or e -iπ (2βt), the upcoming bounds (A.83) and (A.84) are the same as they only involve the modulus of e ±iπ (2βt

  To obtain the asymptotic formula (4.19) as t → ∞, we focus on the first term of Var[Z δ (t)] in (4.17) (or (4.21)) and would like to show that lim

.22) Hence, combining (4.20)-(4.22) results in (4.17).

Table 1 :

 1 Asymptotic mean of Z(t) for t = 100, 000

	Exponential	Pareto
	β Var[Z(t)] η	Var[Z(t)] θ = 0.5 θ = 1.0 θ = 2.0
	0.1	1.004398 0.2 (0.994294)	14,755.3 19,469.8 25,690.5
	0.2	0.399871 0.4 (0.397343)	210.101 365.806 636.906
	0.5	0.126382 1.0 (0.125977)	0.231930 0.463859 0.927718
	1.0	0.055399 1.2 (0.055298)	0.120935 0.262715 0.588618
	2.0	0.025129 1.4 (0.025104)	0.061263 0.133768 0.301616

Table 2 :

 2 Asymptotic variance of Z(t) for t = 100, 000

		Exponential					Pareto
	β	Case 1: s = 10, 000 Case 1: s = 20, 000 t = 100, 000 t = 110, 000	η	Case 1: s = 10, 000 t = 100, 000 Case 2: s = 20, 000 t = 110, 000 θ = 0.5 θ = 1.0 θ = 2.0 θ = 0.5 θ = 1.0 θ = 2.0
	0.1	1.02e-3	1.35e-3	0.2 2,131.32 2,786.52	3,646.48	4,112.99 5,389.03	7,066.23
	0.2	2.55e-4	3.38e-4	0.4	45.4345	73.6313	120.865	73.4282	119.888	198.112
	0.5	4.08e-5	5.41e-5	1.0	8.81e-3	1.82e-2	3.88e-2	1.20e-2	2.48e-2	5.22e-2
	1.0	1.02e-5	1.35e-5	1.2	7.49e-4	1.75e-3	4.16e-3	1.01e-3	2.35e-3	5.54e-3
	2.0	2.55e-6	3.38e-5	1.4	8.01e-5	2.34e-4	7.07e-4	1.06e-4	3.10e-4	9.36e-4

Table 4 :

 4 Asymptotic correlation between Z(s) and Z(t)

  Note that the inequality (A.14) implies the uniform upper bound |Ξ 1 (t, x)| ≤ C Ξ 1 where C Ξ 1 := s 2 α-2 β for all x ∈ [0, s] and t ≥ 2s. Moreover, because t 2 e -βt/2 is bounded on t ≥ 0, (A.[START_REF] Keilson | M/G/∞ with batch arrivals[END_REF]) also means that |Ξ 2 (t, x)| ≤ C Ξ 2 for all x ∈ [0, s] and t ≥ s where C Ξ 2 is some constant. As a result, we arrive at the upper bound |Ξ(t, x)| ≤ C Ξ for all x ∈ [0, s] and t ≥ 2s, (A.[START_REF] Kummer | De integralibus quibusdam definitis et seriebus infinitis[END_REF] whereC Ξ := C Ξ 1 + C Ξ 2 .Term Π(t, x). Concerning Π(t, x) in (A.12), we first consider, by a change of variable v = t(1 -y), -x/t) α | ≤ C * x/t for x ∈ [0, s] and t large enough (say t ≥ t * ) where C * is some constant, we can upper bound (A.20) as

																13) is given by
									lim t→∞	Ξ(t, x) =	x β	for all x ∈ [0, s].	(A.17)
	t	0	1	e -β(t-x)(1-y) y α-1 dy =	0	t	e -βv e βxv/t 1 -	v t	α-1	dv =	0	t/2	+	t/2 t	:= Π * 1 (t, x) + Π * 2 (t, x) (A.19)
	where Π * 1 (t, x) and Π * 2 (t, x) have obvious definitions. Note that Π(t, x) can be rewritten as
						Π(t, x) =	1 -1 -x t 1	α	[Π * 1 (t, x) + Π * 2 (t, x)]t α-2 .	(A.20)
												t
	Applying similar arguments used for analyzing Ξ 1 (t, x) and Ξ 2 (t, x) gives rise to
						lim t→∞	Π * 1 (t, x) =	0	∞	e -βv dv =	1 β	for all x ∈ [0, s]
	by dominated convergence, and
				1												1
	|Π * 2 (t, x)| = t	e -β(t-x)z (1 -z) α-1 dz ≤ e βs te -βt	(1 -z) α-1 dz for all x ∈ [0, s] and t ≥ s,
				1/2												1/2
																(A.21)
	which tends to 0 as t → ∞. Consequently, (A.19) converges to 1/β as t → ∞, and then one uses
	(A.20) to evaluate
					lim t→∞	Π(t, x) t α-2 = lim t→∞	1 -1 -x t 1 t	α	1 β	=	αx β	for all x ∈ [0, s].	(A.22)
		In addition, for all x ∈ [0, s] and t ≥ s, one may check that |Π * 1 (t, x)| ≤ C Π * 1 where C Π * 1 :=
	e βs/2 2 α-1 β while (A.21) implies that |Π * 2 (t, x)| ≤ C Π * 2 for some constant C Π * 2 . Coupled with the fact that
	|1 -(1

  49)This is further split into two cases as follows. Case 1a. η ∈ (2, ∞): In view of the expression on the right-hand side of (A.49), it is noted that (η -2)θxt η-3 dt = θxt η-2 → ∞ as t → ∞, asserting that ψ(t, x) → ∞ as t → ∞. As a result, the left-hand side of (A.46) is in the form ∞/∞, and application of L'Hôpital's rule yields

	lim t→∞	ψ(t, x) t η-2 = lim t→∞	∂ ∂t ψ(t, x) (η -2)t η-3 for all x ∈ [0, s],	(A.50)

which proves (A.46) thanks to (A.49). Case 1b. η ∈ (1, 2): Using the finer asymptotic result (4.30), one finds that (A.42) with v = θ/t can be expressed as

  where C * * * and C * * * * are constants. The equation is also valid when t is replaced by t -x, and therefore taking difference results inψ(t, x) = θ 2 -x θ + (1 -α) ln t -x t + o(1) = -θx + o(1) as t → ∞,which is (A.46).

													as t → ∞.					(A.52)
	We further expand each term in the square bracket as
				θ t	-1 θ t	+ 1	α-1	=	θ t	-1	+ (α -1) + o(1),
				θ t	+ 1		α-2	ln	θ t	= ln	θ t	+ (α -2)	θ t	ln	θ t	+ o	θ t	ln	θ t	= ln	θ t	+ o(1),
			θ t	+ 1	α-2	ln		θ t	+ 1 = o(1),
					C * * θ t	+ 1	α-2	= C * * + o(1).
	Therefore, (A.52) is reduced to
	ϕ	θ t	= θ 2	t θ	-1

-(1 -α) ln θ t + C * * * + o(1) = θ 2 t θ + (1 -α) ln t + C * * * * + o(1) as t → ∞,

  Asymptotics for J * 3 (t). We can now gather the asymptotic results (A.53) and (A.55) to give the asymptotics for J * 3 (t) as t → ∞. First, if η = α, then only J * 31 (t) matters as J * 32 (t) equals zero. Second, if η ∈ (0, α) ∪ (α, 1], then J * 32 (t) asymptotically dominates J * 31 (t). Third, if η > 1, then both J * 31 (t) and J * 32 (t) asymptotically behave like t α-2 and therefore J * 3 (t) ∼ J * 31 (t) + J * 32 (t). These are all summed up and C 1 := α(1 -α) +

				κ+ 1 4
					(A.55)
	µ 2 1 λθ(η -α) (η -1)Γ(α)	0	s	W (s -x)xdm(x) t α-2 ,	if η > 1.

1 λθ(1 -α) Γ(α) s 0 W (s -x)xdm(x) t α-2 ln t, if η = 1,

  ). Hence, we choose to keep the notation e ±iπ for presentation purpose. Thus, the error term ε n,22 (e ±iπ (2βt)) in (A.81) satisfies the inequality|ε n,22 (e ±iπ (2βt))| ≤ 2ω n,22 (t)ζ n,22 (t) (α + n + 1)(1 -α) (2βt) α+n+2 exp 2ω n,22(t)ρ n,22 (t)ζ n,22 (t) 2βt (A.82) for 2βt > 2(n + 1). Note that, contrary to (A.74), the above inequality corresponds this time to e ±iπ (2βt) which belongs to region R 3 ∪ R 3 in Olver et al. (2006, Chapter 13.7(ii)), and hence it is valid for 2βt > 2(n + 1). The functions σ n,22 (t), ω n,22 (t) and ρ n,22 (t) are given by their Equations (13.7.8) and (13.7.9) as σ n,22 (t) := n+1 + n + 1) 2 -2(α + n + 1)(2α + n + 1) + 2α + n + 1| + ν n,22 (t)σ n,22 (t) 1 + ν n,22 (t)σ n,22 (t) 4 [1 -ν n,22 (t)σ n,22 (t)] 2

			2βt ≤ 1 2 ,
			ν n,22 (t) :=	1 2	+	1 2	1 -4[σ n,22 (t)] 2	-1 2	≤	√	2,
		ζ n,22 (t) :=	π 2	+ σ n,22 (t)[ν n,22 (t)] 2 ν n,22 (t) ≤	π 2	+ 1	√	2,
		ω n,22 (t) :=	1 1 -ν n,22 (t)σ n,22 (t)	≤	1 1 -	2 √ 2	= 2 +	√	2,
	and
	ρ n,22 (t) := |2(α ≤ 1 2 1 2 |2α(1 -α) + (n + 1)(1 -2α)| +	19 + 14 4	√	2	≤ B 2 (n + 1) + C 2 ,
	where B 2 := |1-2α| 2	and C 2 := α(1 -α) + 19+14 √ 4	2	are non-negative constants. By substituting the
	above results into (A.82), we arrive at
	|ε n,22 (e ±iπ (2βt))| ≤ 2(1 +	√	2)(π + 2)	(α + n + 1)(1 -α) (2βt) α+n+2	exp	2(1 +	√	2)(π + 2)[B 2 (n + 1) + C 2 ] 2βt
	≤ 2(1 +	√	2)(π + 2)(1 -α)	1 (2βt) α+n+1	α + n + 1 2βt	exp	(1 +	√	2)(π + 2)[B

2 (n + 1) + C 2 ] n + 1

Acknowledgements

The authors would like to thank the anonymous reviewers for helpful comments and suggestions which greatly improved an earlier version of the paper as well as Dr. Olde Daalhuis for insightful discussions concerning Part 1 of the supplementary materials in relation to the proof of Proposition 3.

Appendix A A.1 Proof of Theorem 1

We begin by utilizing (3.3) to expand Z n δ (t) (for n ∈ N and t > 0) into a binomial sum as

Conditioning on T 1 = x, its expectation is given by

e -(n-i)δy dW (y) dF (x)

)]e -(n-i)δ(t-x) T (n-i)δ w(t -x) dF (x).

(A.1)

For convenience, for t > 0 we define g n (t) := e nδt E[Z n δ (t)] and a n (t) := e nδt a n (t), where

Further define the convolution operator * as b 1 * b 2 (t) := t 0 b 1 (t -x)db 2 (x) for t > 0 where b 1 (•) and b 2 (•) have positive domain. Via multiplication of (A.1) by e nδt , one obtains g n (t) = g n * F (t) + a n * F (t) for t > 0, which is a renewal equation satisfied by g n (•). Therefore, its solution is given by

where the last line follows from the representation m(x) = ∞ i=1 F * (i) (x) of the renewal function. Dividing both sides by e nδt yields the desired result (3.4). Γ(α + 1)

By squeezing principle one obtains the asymptotic behaviour

Term J * 2 (t). Again, putting (4.1) and (4.27) into (A.37) yields, in the same spirit as (A.39),

Term J * 3 (t). Utilizing (3.2), one sees that J * 3 (t) in (A.38) can be written as

In order to derive the asymptotics for J * 3 (t) as t → ∞, we proceed by studying the asymptotic behaviour of E[Z(t -x)] -E[Z(t)] for all x ∈ [0, s] and then apply dominated convergence. Let us first define, for v ≥ 0,

which is decreasing in v. By changing variable y = z/(v + 1), we get that

(A.56)

The result. From (A.39) and (A.40), it is clear that J * 1 (t) and J * 2 (t) are both asymptotically proportional to t -η as t → ∞. By comparing with (A.56), one sees that J * 3 (t) is dominated by

Supplementary materials

Part 1: Proving the validity of interchanging limit and infinite summation in Equation (4.23) in Proposition 3

In the proof of Proposition 3 , we have determined the constant C in Equation (4.23), namely

(where A n is defined in (4.18)) by assuming that the limit and the infinite summation on the left-hand side can be interchanged. Here we shall prove that this is valid. Setting

the left-hand side of (A.57) equals lim t→∞ ∞ n=0 χ n (t). To interchange the order of limit and infinite summation, a sufficient condition is that the series

Upper bounding |χ n (t)| by a well defined V n is, in view of (A.58), only possible if one is able to provide some fine asymptotics for the Kummer's function 1 F 1 (a, b, z) as z → ∞ (see (4.10)). We need some precise information on the term O(|z| -1 ), and in particular it is important to analyze how O(|z| -1 ) depends on the parameter b since b = 2α + n + 1 varies with n. Refined expansions are only available in the literature for the Tricomi's functions (to the best of our knowledge) but not for 1 F 1 (a, b, z). Hence, we will use those expansions coupled with the identity (see Olver et al. (2006, Equation (13.2.41)))

where U is the Tricomi's function (or the confluent hypergeometric function of the second kind). We recall that U(a, b, z) is a multivalued complex function, and hence the term e ±iπ , multiplied by z, changes its argument by ±iπ and then changes the value of the term U(b -a, b, e ±iπ z).

We first decompose (A.58) into two parts as

where

In the above definitions, κ > 0 is a constant that will be chosen 'large enough' later on. In other words, we break down χ n (t) into two terms χ n,1 (t) and χ n,2 (t) according to whether t grows slower or faster than O(n) towards ∞, and the terms χ n,1 (t) and χ n,2 (t) will be studied separately. Term χ n,1 (t). First, we look at the term χ n,1 (t). Using (4.18) and (4.8), (A.63) can be rewritten as

It can be readily checked that e -x x α+n (as a function of x for x ≥ 0) achieves maximum at x = α + n.

Therefore, one has that e -2βty (2βty) α+n ≤ e -(α+n) (α + n) α+n . Further noting that tI {2βt≤κ(n+1)} ≤ κ 2β (n + 1), (A.65) can be upper bounded as

and therefore

Then, we apply the ratio test to check that

and thus regarding bounds on U(a, b, z). In particular, we shall put n = 1, a = α or a = α + n + 1, and b = 2α + n + 1 into their Equation (13.7.4). It is noted that the bounds are of slightly different forms depending on the way |z| tends to infinity. In particular, we will use their result by putting z = ±2βt (i.e. z is real), which tends to ±∞ as t → ∞. (t) , and

Recall that κ is to be chosen 'large enough'. Here we shall choose 2βt > κ(n + 1) for some κ that we, at this point, impose strictly larger than 1 (so that the constraint in (A.74) is satisfied) to arrive at the inequalities σ n,21 (t)

Thus, for κ ≥ 2,

Consequently, the ratio test verifies that, for κ ≥ 2,

and hence,

86)

Verifying (A.59) and (A.60). Combining (A.62), (A.70), (A.66), (A.79) and (A.85), it is clear that by choosing κ ≥ 2 in (A.63) and (A.64), we have

i.e. the upper bound (A.59) is obtained by defining the sequence [START_REF] Léveillé | Covariance of discounted compound renewal sums with a stochastic interest rate[END_REF] . Then, the condition (A.60) is satisfied thanks to (A.69), (A.80) and (A.86).

Part 2: Proof of Lemma 1

First, we split G * (y) defined in (4.28) into two parts as

where the first integral is always finite. Therefore, we focus on the divergent second integral and distinguish between the cases η > 1 and η = 1 for ease of presentation. Case 1. η > 1: Performing integration by parts, one obtains

The analysis is further separated into three cases as follows.

Case 1a. 1 < η < 2: In this case, -η + 1 > -1 and therefore the integral

Incorporating such an observation into (A.88), it is clear that

for some constant C. Further expanding

, we obtain

where the last line is due to -η+2 > 0. Combining (A.87), (A.89) and (A.91) yields the first asymptotic result in (4.30). Case 1b. η = 2: Equation (A.88) is still valid, but the integral

x α-2 dx on the right-hand side can easily shown to be divergent at y = 1. Integration by parts gives

where the integral y 1/2 x α-3 ln(1 -x)dx converges as y → 1 -because 1 1/2 ln(1 -x)dx is finite. Incorporating (A.88) and (A.92) into (A.87) and utilizing the expansion y a = [1 + (y -1)] a-1 = 1 + a(y -1) + o(y -1) at a = α -1 and a = α -2, one finds (for some constant C)

which simplifies to the second asymptotic formula in (4.30) because (y -1) ln(1 -y) = o(1). Case 1c. η > 2: The integral on the right-hand side of (A.88) satisfies

By applying the expansion (A.90) (which is also valid for η > 2) and the above result to (A.88), one proves the first result in (4.29) when η > 2.

It is instructive to note that the proven result (4.30) also asserts that the first asymptotic formula in (4.29) is also valid when 1 < η ≤ 2. Case 2. η = 1: Similar to (A.92), integrating by parts on the second integral in (A.87) results in

To analyze the asymptotic behaviour of J(s) as s → ∞ (so that v v+s → 0 + ), we shall study the integral, for small h > 0,

where

1/2 y -γ (1 -y) -ξ dy as h → 0 + by dominated convergence (as G(1 -h y ) is bounded on h ∈ (0, 1 2 ] and y ∈ [ 1 2 , 1] and

We need to consider two cases as follows.

Case 1. γ -2 > -1 ⇔ γ > 1: We rewrite (A.96) as

where ω(z, h)

It can be seen that ω(z, h) is bounded by some constant for all z ∈ [0, 1] and a small enough h. Besides, one has that (1 -2 -t h 1-t ) -ξ G(1 -(2h) t )dt.

(A.99)

For t ∈ [0, 1] and h small enough (less than 1/2), it is observed that both (1-2 -t h 1-t ) -ξ and G(1-(2h) t ) are upper bounded by some constant. Therefore, one checks easily that

Hence, by the dominated convergence theorem together with the fact that -ln 2h = -ln 2 -ln h ∼ -ln h as h → 0 + , one finds that the asymptotic behaviour of (A.99) (or (A.96) with γ = 1) is given by U 1 (h) ∼ (-ln h)G(1) as h → 0 + .

(A.100) Therefore, combining the two cases from (A.98) and (A.100), we have that U 1 (h) is O(h 1-γ ) for γ > 1 but O(-ln h) for γ = 1. Since U 2 (h) in (A.95) is O(1), it can be concluded that (A.94) asymptotically behaves like, as s → ∞,

Finally, plugging h = v v+s in (A.93) results in (4.34).