- Email: [email protected]
Computational methods in risk theory: A matrix-algorithmic approach
Insurance: Mathematics North-Holland
and Economics
259
10 (1991) 2.59-274
Computational methods in risk theory: A matrix-algorithmic approach Soren Asmussen Institute of Electronic Systems, Aalborg Unicersity, DK-9220 Aalborg, Denmark
Tomasz Rolski * Mathematical Institute, Wroclaw University, 50-370 Wroctaw, Poland Received
April
1991
Abstract: This paper is concerned with the numerical computation of the probability t)(u) of ruin with initial reserve u. The basic assumption states that the claim size distribution is phase-type in the sense of Neuts. The models considered are: the classical compound Poisson risk process, the Sparre Anderse process and varying environments which are either governed by a Markov process or exhibit periodic fluctuations. The computational steps involve the iterative solution of a non-linear matrix equation numerical examples are presented. Q = cY(Q) as well as the evaluation of matrix-exponentials e QU.A number of worked-out Keywords: Ruin probability, ment. Periodic environment.
Phase-type
distribution,
Matriu-algorithmic
methods,
Non-linear
matrix
iteration,
Markovian
environ-
1. Introduction
This paper is concerned with the numerical computation
of infinite horizon ruin probabilities
(1-l) where (S,} is the claim surplus process, u is the initial reserve and T.(U)= inf{t L 0: S, > u} is the time of ruin (first passage time to level u). Thus N
s, =
c
l_J-pt,
i=l
where IV, is the number of claims up to time t, the Ui are claim sizes and p is the premium rate. For example, in the classical compound Poisson model, {IV,}is a homogenous Poisson process and the Ui are independent of {IV,) and i.i.d. with common distribution B; even in this simple case, the numerical computation is non-trivial. At the core of our discussion is the fact that risk theory and queueing theory are mathematically closely related areas. For the compound Poisson model with p = 1 the relation is t/l(u) =
qw>u>,
(1.2)
where W is the stationary waiting time in the M/G/l queue with the same arrival process and service times distributed as the claim sizes. The quantities asked for in queueing applications are, however, * The second
author
0167-6687/92/$05.00
was supported
in part by The Philip Lundberg
0 1992 - Elsevier
Science
Publishers
Foundation
B.V. All rights reserved
S. Asmussen,
260
T. Rolski
/ Computational
methods
in risk theory
typically more simple than the whole distribution of W as in (1.2). Most often only the mean (El+’or the first few moments are required, and in the M/G/l case these are explicit in terms of moments of U. The GI/M/l queue is one more example, where explicit solutions are available, but more general systems present difficulties. In queueing theory, the development of computational methods has, however, been rapidly progressing during the last one or two decades, due not least to the application potential in areas like computer networks. Of particular importance for our discussion is the matrix-algorithmical method initiated by Neuts [see Neuts (1981, 1989)], and the extensive bibliographies therein. Briefly speaking the purpose of this paper is to explain that this approach is relevant in risk theory as well, to carry out the necessary adaptation, and to present some worked-out numerical examples. Much of the literature on ruin probabilities deals with approximations, for example the celebratedCram&-Lundberg approximation
(1.3) Our concern is here, however, the computation of the exact values. The basis assumption, which makes this possible, is restricting the claims to be phase-type; the crucial fact is that phase-type assumptions on the claims again provide a phase-type structure of the ruin probabilities. The phase-type set-up covers for example finite mixtures of exponential distributions and Gamma distributions with integer parameter. Whereas many papers in the literature deal with such examples, the general class of phase-type distributions have hardly been exploited at all in risk theory. Essentially, the distributions in the literature which do not fit into the framework are heavy-tailed ones like the lognormal [see Embrechts and Veraverbeke (1982)]. We start the paper by a brief summary of phase-type distributions in Section 2. Section 3 then deals with the classical risk model with Poisson arrivals of claims. For this case, the relevant expression for e(u) is a translation of a queueing result given in Neuts (1981) and, despite the simplicity of the solution and that it is well-known in queueing theory, it covers many special cases treated separately in the literature on ruin probabilities. The rest of the paper then relies on some more recent advances in algorithmical queueing theory, see in particular Asmussen (1991, 1992) Sengupta (1989, 1990), Ramaswami (1990) and Lucantoni (1990). The essential feature is that whereas largely the original matrix-geometric method as surveyed in Neuts (1981) is best suited for studying discrete processes (like the numbers of customers present in the queue), the papers just listed are able to deal directly with continuous-state problems. In Section 3, we study the risk process with renewal arrival of claims (sometimes called the Sparre Andersen process); the material is a translation of Asmussen (1992) which in turn has some features similar to Sengupta (1989). Section 4 deals with risk processes in a Markovian environment as studied by Asmussen (1989); compared to that paper which deals largely with approximations, the class of models for which we can present an explicit (or at least computationally tractable) exact solution is considerably widened. Finally, in Section 5 we give some remarks on models with a periodically varying Poisson intensity and how they can be approximated by a Markovian environment model. We remark at the end, that from the point of risk theory it may not feel all that satisfying to derive the results and algorithms by translation of the queueing literature. However, at least within the framework of this paper, we see both areas as special cases of random walks and their extensions (e.g. Markov additive processes); the approach that we present just happened to be initiated in queueing theory rather than in risk theory. A more self-contained (and therefore also more lengthy) exposition will be given elsewhere [Asmussen (1992)l. In the rest of the paper we scale for simplicity the claim surplus process to have premium rate p equal 1 and suppose that the safety loading 77 is positive; this implies that the surplus process tends to --OOwith probability 1. The traffic intensity p of queueing theory is related to 77 by p = l/(1
+ 7~)
or
77 = (l/p)
- 1.
(1.4)
S. Asmussen, T. Rolski / Computational methods in risk theory
2. Phase-type
261
distributions
That the distribution B of the generic claim U is of phase-type with representation that the d x d matrix T is the restriction to 11,. . . , d} of the intensity matrix
(r, T, d) means
of a Markov jump process on IO,. . . , d) with 0 as absorbing state, and that U is distributed as the time to absorption in 0 given the initial distribution r (written as row vector in the following). In particular, the vector t, of exit rates is t, = - Te where e is the column vector of ones. We assume that T is a proper subintensity matrix. For example a hyperexponential density (i.e. a finite mixture of exponential densities) Cfrjai epalx is phase-type with T = (a,, . . ., adjdiae, and an Erlang density adxd-‘/Cd l)! e- uX is phase-type with representation I-,
(Y --(Y
0 p=(l
Fundamental 1 -B(x)
0
..a
o),
. .
T=
0 G? ‘.
.
... ...
0 0
0 0
.
(2.1)
(j
(j
0
..:
_La
0
0
0
*.*
0
b(x)
=B’(x)
CY
-a
analytical identities are =lmB(dy) X
=reTxe,
=reTXtO,
(2.2)
T)-‘to,
(2.3)
-1
e=r(
pL(g)= on”B(dn) /
= ( - l)“n!rT-“e.
-sZ-
(2.4)
See Neuts (1981) or Asmussen (1987, Ch. III 6) for more detail. A crucial fact of these formulae is their amenability to numerical computations: all that is needed is (a) matrix inverses [the standard algorithm is Gauss-Jordan elimination with full pivoting, cf. Press et al. (198811 and (b) matrix exponentials for which a variety of computational recipes are available as well, see the appendix. A further important property of phase-type distributions is denseness: any given distribution on (0, w) can be approximated arbitrarily close by a phase-type distribution. This point is illustrated in Figures 1 and 2, where we have plotted phase-type fits to two truncated lognormal distributions, more precisely the distributions of eUU-” I2 A 2 where U is standard normal. These examples are the ones to be used in the rest of the paper, and thus some further comments are in order. Example 2.1 (Figure 1). The choice u = 1.8 in Figure 1 is from Benckert and Jung (1974) who found this to provide an adequate fit to some Swedish fire insurance data, and the distribution appears again in Thorin and Wikstad (1973); the truncation at 2 could correspond, e.g., to reinsurance. The parameters of the fitted phase-type distribution B, (say) in Figure 1 are d = 2,
.H = ( 0.5614
0.4386))
T=
- 8.640 0.101
(2.5)
The mean is the same as for the truncated log-normal, viz. 0.6015 (note that the unrestricted log-normal distribution has mean 1 and thus the truncation matters more than a brief look at Figure 1 might
262
S. Asmussen, T Rolski / Computational methods in risk theory
i
3
-
- -
- -
lognormal;
(T =
1.8
Bl
2
1 :\\
-4
--__
I
-----mm-__
I 2
1 Fig. 1.
indicate). The squared coefficient of variation (s.c.v.) is 1.8988. The standard comparison is the exponential distribution, where the S.C.V.is 1, and since a rough rule of thumb asserts that the risk increases with any added stochastic variability, we thus expect larger ruin probabilities than for exponentially distributed claim sizes with same mean.
Example 2.2 (Figure 2). The choice u = 0.8 in Figure 2 was made to obtain a distribution with a different shape than in Figure 1, but is otherwise arbitrary. The parameters of the distribution B, on Figure 2 are d = 4, m= (0.9731 T=
i
0.0152
- 28.648 0.102 0.113 0.100
0.0106
28.532 - 8.255 0.107 0.102
0.00 1
0.089 8.063 - 5.807 0.111
(2.6)
3-
--
2-
-
- -
lognormal;
0 = 0.8
B2
1
Fig. 2.
2
263
S. Asmusen, T. Rolski / Computational methods in risk theory
The mean of l3, is 0.8849 and the S.C.V.is 0.5132, and thus we can expect smaller for the exponential case.
ruin probabilities
than
Some basic references for fitting phase-type distributions are Bux and Herzog (1977), Johnson and Taaffe (1989) and Bobbio and Cumani (1990); the distributions B,, B, were obtained by the method of Asmussen and Nerman (1990). It should be noted, however, that our aims here are not - to survey just a few _ to discuss - to neglect distribution.
statistical aspects of phase-type fitting, only to point out that phase-type distributions with phases may catch the shape of a given (reasonably smooth) claim size distribution very well; whether phase-type distributions are appropriate for fitting heavy-tailed distributions; the problem that ruin probabilities are sensitive to even small changes in the claim size
Nevertheless, we may argue that if one wants to distinguish between lognormal in the figures, one needs a substantial amount of statistical data, presumably most other parametric models would fail a goodness of fit test anyway. With to explain how to calculate ruin probabilities using phase-type assumptions, purposes in the rest of the paper through (2.5), (2.6). For our computational the following structure of the ruin probabilites 4(u) will be found for a variety claims:
or phase-type distribution so much that these and these remarks, we proceed exemplified for numerical approach it is crucial that of models with phase-type
(u, S, d) if for some d-dimensional Definition 2.1. A function I,%(U)is of phase-type with representation subprobability vector u and some d x d subintensity matrix S we have I++(U)= v esue.
3. The classical compound
Poisson model
For the classical compound Poisson model with claim size distribution B of phase-type and arrival intensity p, the corresponding queue is M/PH/l and was considered by Neuts (1981, pp. 52-57). In risk theoretical terms the result is the following: Proposition
3.1.
T, d),
the ruin probability is phase-type with representation CT+, Q, d), where and
Q=T+tp+,
T+= -@rT-‘.
Proof. This follows directly by combining follows. First note that by standard matrix -rT-r
=
mrer* /0
(3.1) (1.2) and Neuts calculus
(1981). A short self-contained
dx
proof goes as
(3.2)
[the convergence in (3.2) is ensured by the Perron-Frobenius theorem because T is a proper subintensity matrix] so that the components of rr+ are nonnegative and the total mass IIr+ II = r+e of r+ is I(rr+ 1)= IIpuma
eTx dx II = p&)
Next, by the Pollaczek-Khinchine P(W4u) where
pB, P
= (1 -p)
PB
formula
[see e.g. Asmussen
(1987)1, (3.3)
5 p”Bo*“(u), n=O
is the (defective)
1 -B(x) (1)
= p < 1.
distribution
= /?P eTxe = p( -rT-‘)
with density e”(
- Te) = r+ eTXtO.
S. Asmussen, T. Rolski / Computational methods in risk theory
264
_______.
B
1
. . . . . . . . . . . . . . . . . . . . exponential
25’EU
WEU
Fig. 3.
the That is, pB, is defective phase-type with representation (P,, T, d). In (3.3) we can thus represent distribution of W by a geometric-compound of a (P,, T, d) phase-type distribution which is phase-type q (a,, T + c,,r+, d) as is readily seen by a probabilistic argument. Example 3.1. We shall reinspect an example in Gerber arrival intensity is p = 3 and the claim size density is
(1979) in view of Proposition
3.1. Here
the
b(~)=+.3e-~~++.7e-~~. Thus
B is hyperexponential
r+=
The eigenvalues
Q(U)
=
in accordance
with r = (+ i), T = (-3
(3
eQ"e =
.
Q=T+t,,r+=
&),
of Q are - 1 and
r+
- 7)diag, so that
- 6, and standard
diagonalization
techniques
yield
$ edu + & ep6’
with Gerber
(1979).
For numerical illustration, we fix from now on the safety loading 77 to 10% (so that p = 0.9091) and consider $(u) in the range 0 I u _<50 p, where p is the average claim size (for this simple model, we have just p = [EU). The ruin probability @l(u) is plotted on a logarithmic scale in Figure 3 for Poisson 2.1, 2.2 and the arrivals and three different claim size distributions, namely B,, B, of Examples exponential distribution with mean one (note that the arrival rate /? is determined by either of n, p sot that for B, we get p = p/p = 1.496 and for B, we get p = 1.017). The basic parameters of (3.1) were calculated as d = 2, ?r+= ( 0.10757
0.80152))
Q = ( - ;:;;;;
_;‘;;;;)
S. Asmussen, T. Rolski / Computational methods in risk theory
265
for Example 2.1, and as d = 4, r+=
(0.038
0.143
( -28.648
Q= \
0.102 0.124 0.171
0.210
28 532 - 8.254 0.149 0.368
0.518), 0.089 8.064 - 5.746 0.503
0.027 0.088 5.446 - 1.210
for Example 2.2. It is seen from Figure 3 that our expectations of which claim size distribution has the higher risk were confirmed. The almost linear appearence of the curves on Figure 3 and in later figures of the paper reflect the fact that the ruin probabilities are close to being exponential (the scale is logaritmical). There are essentially two reasons for this. The Cramer-Lundberg approximation (1.3) implies that the curves have linear asymptotes; thus possible deviations from linearity have to be sought for small values of U. However, also here the curves look almost linear, which we would explain by the fact that n = 10% is rather small so that we are in the domain of the diffusion approximation, cf. Grandell (1977).
4. The Sparre Andersen process In this section we suppose that the number of claims up to time t is a renewal process with interarrival distribution A and that the claim sizes are i.i.d. and independent of IN,] with common distribution B of phase-type with representation (w, T, d). The renewal assumption was introduced by Sparre Andersen (1957) and extensively studied by Thorin (1970, 1971, 1974, 1975) and Thorin and Wikstad (1973, 1974). We consider here two cases. In the first one we suppose that the process starts at a claim arrival (the arrivals form a non-delayed renewal process) and in the second case [Thorin (1975)] that the arrival process is a stationary renewal process. The ruin probability for the non-delayed case is denoted by I/J(U) and for the stationary case by 4*(u). The corresponding queue is GI/PH/l and was studied by Asmussen (1992) [see also Sengupta (1989)]. Before we recall the result from Asmussen (1992), we introduce some notation. For a subintensity matrix K (this provides the convergence), define AIK] = 10”eKXA(dx), and let V, be the nonlinear matrix operator defined by V,(K)
= T-t&&K],
(4.1)
where t, = - Te. Proposition 4.1. For the renewal model with claim size distribution of phase-type, the rain probability $(u) is of phase-type CT+, Q, d), where the subintensity matrix Q solves the non-linear matrix equation
Q = qd’l(e)
(4.2)
and rr+ = e’(Q - T)/e’t,; the ruin probability I/I*(u> is of phase-type (p y, Q, d), where y is the distribution soluing y(T + t,lr) = 0, ye = 1. Proof. The result concerning $(u) follows from Asmussen (1992) and the common interpretation of $(u) and p(W > U) as P(M > U) where M is the maximum of the random walk with generic increments distributed as the independent difference between a claim size (service time) and an interarrival time. For ICI*(u), one notes similarly that e*(u) = P(V> u), where V is the virtual waiting time of the GI/PH/l queue. 0
It is important
to notice that Q may be computed
Q n+l = !P,(Q,>, (n = 0, 1, . . .I.
as the limit Q = lim, +_Q,,, where Q, = T and
Two interarrival distributions, A,, A, were considered, one phase-type and one deterministic, respectively. We shall explain next how to carry out the crucial steps of calculating AIK] in (4.1). Consider first the phase-type case. Then (see the appendix for the Kronecker notation):
266
S. Asmussen, T. Rolski / Computational methods in risk theory
Proposition
d[K]
4.2.
Zf A is phase-type with representation (Y, S, q), then
= (z@v)(-K6Ls~‘(zc3s”),
(4.3)
where sg = -Se. Proof.
We have eKXveS”s,=(Z~Y)(eKx~esx)(Z~Sg)
(4.4)
and hence using eKx @ eSXe(“*‘jX we get
a[K]
=
irn eKxv eSXsOdx
= (ZC3v)im(eKr
@esX)
dx(Z@ss,)
= (Z@pY)( -KfBS)-‘(z@ss,).
0
(4.5)
In the deterministic case, where A, is degenerate at l/p (say), we have AlK] = eKla; to compute the matrix exponential, we used uniformization, see the appendix. We performed computations for the claim size distribution B, and the following interarrival distributions: A, - the hyperexponential distribution with density 0.4 .5c~ e-5a + 0.6. (yepa; A, - the distribution degenerate at l/p; and, for the sake of comparison, also the exponential distribution (which is already included in Figure 3). The scaling factors a, p were determined by 77= 10%. The illustration is in Figure 4. Again, the comparison of the risks involved in the different cases correspond nicely to the fact that the hyperexponential distribution has a S.C.V. which is higher than for the exponential distribution which is in turn higher than for the deterministic case. The basic parameters of Proposition 4.1 are d = 2, r+=
( 0.159
for the hyperexponential 7r+= (0.013
Q= ( -;$“6
0.780))
inter-arrival
-
t I
times, whereas for the deterministic
Q= ( -;$;
0.814),
0.001
_;‘i;;)
-
I
-
I
_;‘;;;).
-
I
case we got
-
I
-
-
I
25’EU Fig. 4.
__-_____
B
....................
exponential
-
I
-
I
-
I
1
_
-
I
50’EU
267
S. Asmussen, T. Rolski / Computational methods in risk theory
5. Markov modulated
arrivals
In this section we suppose that the environment is governed by an irreducible Markov process with states i= l,..., p and intensity matrix A = {hii}, say. This means that if the governing process (environment) is in state i, then claims are Poisson with intensity pi and claim sizes U”’ have common distribution C”’ of phase-type with representation (rci), Sci), 9”‘) (i = 1,. . . , p). For such a process the average arrival rate is p = Zip_, yipi, where y is the probability solution of the stationarity equation r/i = 0, the average claim size is
(5.1) and the average amount of claims per unit time is p = pp. Thus, as before, (1 -p)/p, and letting pi = p,lEU”‘, ni = l/p, - 1, we can rewrite (5.1) as P =
5
YiPi
,)./ =
i=l
5
ESi.
i=l
P
the safety loading is
(5.2)
This formula expresses overall characteristics of the Markov-modulated process in terms of the p standard compound Poisson risk processes that we obtain by fixing the parameters of the ith to be pi, Cc’). Note, however, that whereas we need to assume that the overall safety loading n is > 0, it is possible that vi < 0 for one or more i; such states i represent environmental conditions which are unfavourable in the sense that when present, a loss is suffered. The model has been studied in risk theory by Janssen and Reinhard (1985), Janssen (1980), Reinhard (1984) and Asmussen (1989), whereas some main early queueing references are Burman and Smith (1986) and Regterschot and de Smit (1986). We follow here Asmussen (1991). Note that the connection between the queue and the risk process is slightly more complicated than in (1.2) by involving a time-reversion of A; however, the basic fact is that both the waiting distribution for the queue and the ruin probabilities can be expressed in terms of maxima of a certain continuous-time random walk with Markov-dependent increments. We denote by J/i(~) the ruin probability subject to the initial condition that the reserve is u and the environment is i(i= l,..., p). Proposition 5.1. For a risk process in a Markovian environment with phase-type claim size distributions, the ruin probability *i(U) is phase type with representation CT?), Q, E) where: p, j’=l,...,
E={(j,j’):j=l,...,
q(j)},
and p~j;=Bj(el~rr(j))(-K~S(j))-‘(ej~Z), where the K is the solution of the nonlinear matrix equation K+A-(P)diag+
W(K).
(5.3)
Here W(K) denotes the matrix with the i&h element
The result is proved in Asmussen (1991) with W(K) = lFeKXs(dX), where S(dx) denotes integration of the jth column with respect to /3jC”‘(dx); the further calculations are as in Proposition 4.2. q Proof.
268
S. Asmussen,
T. Rolski / Computational
methods in risk theory ecn = 24 . . . . . . . . . . . . . . . . . . . . ecn = 16 I
I
O.Ol-------__-___
0.001
-
-
-
I
-
-
-
I
I
-
I
-
I 25’EU
-
I
-
I
_
I
_
_
I 50’EU
Fig. 5.
We performed computations for the following (thus 4 (I) = 2 , q(*) = 4) and the intensity matrix
data: p = 2, PI = 3.325, & = 0.781, C(r) = B,, C’*’ = B,
The value of & is chosen to obtain 7, = - 50% (or equivalently pr = p, pr = 2) and p2 is then determined from n = 10% by equation (5.2). The role of the parameter A is as a measure of the degree of Markov-modulation. As A + co, the Markov-modulated process tends to the homogeneous Poissonian claim arrival stream [this was noted in Asmussen (1989) and proved formally in Rolski (199O)l whereas for A + 0, the initial state becomes predominant. The three curves of the ruin function in Figure 5 correspond to three different values of h, which in turn were determined by letting the mean number ecn = - C?= ryipiAz; ’ of claims till completion of an environmental cycle be ecn = 4, 16, 64, respectively. Thus, the case of eta = 64 should correspond to an almost homogeneous Poissonian input, whereas the case ecn = 4 should be almost equivalent to the classical compound Poisson model with parameters given by the initial condition. The basic parameters of Proposition 5.1 computed for case ecn = 4 were K= P:)=
P$?=
Q=
0.511 2.554
-0.246 1.229 1 ’
( TT$~, : j = 1, 2, j’ = 1,. . . , q(j)) = (0.173
0.575
0.002
0.018
0.042
0.142),
(0.013
0.238
0.029
0.105
0.151
0.365 ),
1 -7.489 0.273 0.000 0.000 0.004 \ 0.024
5.813 - 0.525 0.000 0.001 0.069 0.443
0.012 0.012 - 28.648 0.102 0.121 0.153
0.122 0.018 28 532 - 8.255 0.137 0.297
0.278 0.042 0.089 8.064 - 5.763 0.393
0.945 0.141 0.027 0.087 5.402 - 1.495
269
S. Asmussen, T. Rolski / Computational methods in risk theory
The ordering environmental
of the states in E is lexicographical, so that e.g. the fifth element state 2 and phase 3 of B,. For ecn = 16 we have
K= ( 0.318 1.598
0.709
0.001
0.006
0.015
0.056),
P?‘= (0.009
0.211
0.029
0.107
0.157
0.383),
’ -7.380
Q=
0.289
6.709 - 0.390
0.000 0.003 0.018
0.001 0.061 0.392
0.000
\
to
-0.070 0.347 1’
(0.190
4!)=
corresponds
0 .ooo
0.004 0.001 - 28.648 0.102 0.121 0.154
0.042 0.006 28.532 - 8.255 0.138 0.301
0.375 \ 0.056 0.027 , 0.088 5.407 - 1.464 /
0.100 0.015 0.089 8.064 - 5.761 0.403
and for ecn = 64, K=
- 1.317 0.263
rY)= (0.195
0.774
0.000
0.002
0.004
0.017))
@=
(0.008
0.193
0.029
0.108
0.159
0.391),
Q=
l-7.342 0.295 0 .ooo 0 .ooo 0.002 \ 0.015
7.138 - 0.326 0.000 0.001 0.057 0.368
0.012 0.002 28.532 - 8.255 0.138 0.303
0.001 0.000 - 28.648 0.102 0.121 0.154
0.029 0.004 0.089 8.064 - 5.761 0.407
0.112 ’ 0.017 0.027 0.088 . 5.410 -1.4491
In Figure 5, the initial condition is state 1 of the environment, whereas in Figure 6, we have plotted @r(u) as well as I/J&U) for the intermediate case ecn = 16. Note that Figure 6 gives an example where the non-exponentiality of the ruin probabilities for small u is apparent; the version of the Cram&-Lundberg
ecn =I6 .. . . . . . . . . . . . . . . . . . .
0.1
-
____-------
0.01
-
___--------
0.00,
_
_
_
I
-
I
-
I
-
-
I
I 25’EU Fig. 6.
-
-
I
-
I
-
I
y
-
-
I 50’EU
270
S. Asmussen, T. Rolski / Computational methods in risk theory
approximation given in Asmussen (1989) implies that the asymptotic linear slope of the two curves is the same. Figure 5 reconfirms the findings of Asmussen (1989) (and much of the literature on similar models): the risk is considerably increased as the environment becomes more variable.
6. Periodic environment The claim arrival intensity may depend on some seasonal effect, in which case we can model the stream of claims by a periodic Poisson process. We scale for simplicity the model to have period 1 and assume that the claim sizes (Ui} are i.i.d. with common distribution B. Then, if we ask for the ruin probability starting at a given time of season S, N, = N”(t) is a Poisson process with a periodic arrival rate p’(s + .). In this case we denote the ruin probability by I/J”(U), with the superscript s E [O, 1) indicating the time of the initial season. It can be written as V(u)
=l+;(
;+‘i
)Bj,
(6.1)
where for s E (0, 11, {M”(t)} is the nonhomogeneous Poisson process with intensity function p’(s + t). Note that for such a process the average arrival rate is /? = l@‘(t) dt. The model has been studied, for example, in Dassios and Embrechts (1989) and Asmussen and Rolski (1991). In the queueing literature, a few papers are devoted to queues with periodic input stream; see Harrison and Lemoine (1977), Lemoine (1981, 1989), Afanas’eva (1985), Rolski (1981, 1987, 1989) and Asmussen and Thorisson (1987). A characteristic of the theory of periodic queues theory for a periodic arrival rate function /3*(t) is the stationary workload I”“’ at given time s E [O, 1) of the season which, following Lemoine (1981), has representation M,(I) VS)
= sl’:
c
i where {M,(t)] /3’(t) =P*(-t)
)
q-t
i=l
(6.2)
1
is the nonhomogeneous we see that
Poisson
process
with intensity
l/V(U) = P{v”-“‘>U}.
function
/3*(s - t). Thus
setting
(6.3)
Rolski (1987) proposed to approximate queueing characteristics of periodic queues by Markov modulated ones. In risk theoretical terms his result translates via (6.3) to the following. For each rr = 1, 2, . . . define the nth approximating risk process as the Markov modulated risk process with the environment governed by a Markov process having states 0,. . . , n - 1 and intensity matrix A given by ‘---II 0 \ n
...
n -n . .
0 n . .
**. . . . .
()
(j
.::
0 ’ 0 (6.4) -n,
Arrivals occur at rate pi = p’((i each approximating risk process function for the corresponding continuous intensity functions /I’ Proposition
6.1.
Zf i/n
1)/n) if the governing Markov process is in state i. The claim sizes at are the same as in the original risk process. Let +,!“‘(u) be the risk nth risk process. The following holds for continuous or piecewise [see Rolski (1987)]:
+ s E [O, l), then I/J:“‘(U) -+ I/J”(~).
Note that it seems intuitively plausible that the above result can be extended to the more general model with claim size distribution B”(t) depending also on the season s E [O, 1). For a proof, see Asmussen and Rolski (1991).
S. Asmussen, T. Rolski / Computational methods in risk theory
271
....................
n =
-_-_-mm.
n= 16
8
_-----
0.00,
_
_
_
I
-
I
-
I
-
-
I
I 25’EU
-
-
I
-
-
I
I
-
-
I 50’EU
Fig. 7.
We performed p’(t)
computations
for the intensity function
= p( 1 + sin 27rt)
(6.5)
and exponentially distributed claim sizes with mean claim size l/6, say. The risk function I/#‘(U) for the approximating Markov modulated models is plotted in Figure 7 for n = 8, 16, 24. Here we chose p = p and S = 1 to have safety loading n = 10%. To demonstrate seasonal variability, we plotted in Figure 8 the function 4”(O), (0 IS < 1) approximated by I@~‘(O) (i = 0,. . , , 23), taking p = 2 or 10 and S determined by condition that n = 10%. The greater variability for the case p = 10 is explained by the fact that p is the average number of arriving claims per unit interval and that increasing fi makes the claim sizes smaller but more numerous. Thus, the stochastic variability in the claim sizes and the claim arrival times matters less compared to the seasonal variability of P’(t). Because of the general properties of the Markov-modulated model, it is not surprising that the curves again are almost exponential. However, under some regularity conditions, the periodic model also allows a version of the Cramer-Lundberg approximation. It is even possible to find adjustment coefficients y,, for the Markov modulated models as well as y for the periodic case. More precisely, it can be seen from
1.o m II0 l
. . .
0.9
.O
~!a.
. l *
0.6
Fig. 8.
p=2 p=10
272
S. Asmussen, T. Rolski / Computational methods in risk theory
Asmussen (1989) that for the nth approximating n
I + Yn+ (&Y,)
- l)P’((i
- 1)/n)
model ‘y, is the positive root of the equation = I
n
i=l ni
7
(6.6)
1 [here g(a) = /r e”“B(dx)] and y is the positive root of the equation Y =P(&,) - I)> (6.7) that is, the same as for the corresponding standard compound Poisson model with constant arrival intensity p. From (6.6) and (6.7) it may concluded that yn + y and moreover one can gain some information concerning of the speed of convergence. For details of these facts, see Asmussen and Rolski (1991).
7. Appendix We discuss here the basic procedures
used in our computations.
1. The exponential of a matrix.
There are various methods available for computing of matrix exponentials eK, see for a review Moler and Van Loan (1978). If K is a subintensity matrix, one further alternative not treated in Moler and Van Loan (1978) is popular in applied probability, namely uniformization (the procedure is similar in spirit to the method called scaling and squaring, which was recently used in Sengupta (1989) for similar purposes). For a subintensity matrix K = (kij), let n = 11K II = max,,j I kij I. Then IIZ + 7-lK II I 1 and
iu
,,,k
=c k=O
k!
e-"(I+
T-‘K)~
+Error,
(7.1)
where IIError II 5 C~=i,+l(~k/k!) e-7. This method was used to compute AIK] = eK/P for the Sparre Andersen process with deterministic arrivals, cf. Section 4. However, for the ruin probability function +(u) the numerical computations in this paper are based on the solution of a linear system of differential equations, for which we used the standard fourth-order Runge-Kutta method [see e.g. Press et al. (1988)]. More precisely, we need to evaluate functions of the form r eKth, where r is a row vector and h is a column vector. We note that the vector z(t) = eK’h fulfills the linear system of differential equations z’ = Kz, with the initial condition z(O) = h. This procedure appears particularly efficient whenever we need to compute the function r eKth = PZ for several values of t, whereas the before mentioned method is more convenient for computing a single or a few values. 2. Kronecker calculus. We recall that if A(‘) is a k, X m, and A’*’ a k, X m2 matrix, then the Kronecker Further(tensor) product A”’ C%J A’*’ is the (k, x k,) x (m, x m,) matrix with (i,i,X j, j,)th entry a$,a$. more, if k, = m, and k, = m2, then the Kronecker sum is defined by
A(‘) @A(*) = A(‘) @Zk*+ Zk, C3A’*‘.
(7.2)
See, e.g., Graham (1981). The following formula is basic for this paper: Suppose that A and B are both square such that CY+ /? # 0 whenever LYis an eigenvalue of A Lemma 7.1. and p is an eigenvalue of B. Let further 7r, v be any row vectors and h, k any column vectors. Then /0
XsreA’h~veB’kdt=(sr8v)(A@B)-‘(eA”BX-Z)(h@k).
(7.3)
S. Asmussen, T. Rolski / Computational methods in risk theory
273
If moreover (Y + p < 0 for all eigenvalues, then mreAthm?tk I0 Proof.
dt=
To prove (7.3) it suffices to differentiate
r eAth . v esrk = (T @v) eAQBx( h 8 k), see Graham (1981).
(7.4)
-(v8v)(A$B)-‘(h&k). both sides and verifjl the identity
P-5)
0
References Afanas’eva, L.L. (1985). On periodic distribution of waiting-time process. Stability Problems for Stochastic Models (Uzhgorod, 1984) l-20. Lecture Notes in Math. 1155. Springer, Berlin. Andersen, E. Sparre (1957). On the collective theory of risk in case of contagion between the claims. Transactions XVth International Congress of Actuaries., New York, II, 219-229. Asmussen, S. (1984). Approximations for the probability of ruin within finite time. Stand. Actuarial. J. 31-57; ibid. (1985) 64. Asmussen, S. (1987). Applied Probability and Queues. Wiley, Chichester. Asmussen, S. (1989). Risk theory in a Markovian environment. Stand. Actuarial J., 69-100. Asmussen, S. (1991). Ladder heights and the Markov-modulated M/G/l queue. Stoch. Proc. Appl. 39, no. 2. Asmussen, S. (1992a). Phase-type representations in random walk and queueing problems. The Annals of Probability 20, forthcoming. Asmussen, S. (1992b). Fundamentals of Ruin Probability Theory. Incomplete book manuscript. Asmussen. S. and 0. Nerman (1990). The EM Algoritm for phase-type distributions viewed as incompletely observed exponential families, in preparation. Preliminary version published in: K. Vest Nielsen, ed., Symposium i Anvendt Statistik (1991). UNI-C, Copenhagen, 335-346. Asmussen, S. and T. Rolski (1991). Risk theory in a periodic environment: The Cram&r-Lundberg approximation and Lundberg’s inequality. Submitted for publication. Asmussen, S. and H. Thorisson (1987). A Markov chain approach to periodic queues. J. Appl. Probab. 24, 215-225. Benckert, L.-G. and J. Jung (1974). Statistical models of claim distributions in fire insurance. ASTIN Bull. VIII, l-25. Bjiirk, T. and J. Grandell (1988). Exponential inequalities for ruin probabilities in the Cox case. Stand. Actuarial .I., 77-111. Bobbio, A. and A. Cumani (1990). ML estimation of the parameters of a PH distribution in triangular canonical form, Technical report. Burman, D.Y. and D.R. Smith (1986). An asymptotic analysis of a queueing system with Markov-modulated arrivals. Opns.Res. 34, 105-119. Bux, W. and U. Herzog (1977). The phase concept: approximation of measured data and performance analysis. In: K.M. Chandy and M. Reiser, eds., Computer Performance. North-Holland, Amsterdam, 23-38. Dassios, A. and P. Embrechts (1989). Martingales and insurance risk. Stochastic Models 5, 181-218. Embrechts, P. and N. Veraverbeke (1982). Estimates of the probability of ruin with special emphasis on the possibility of large claims. Insurance: Mathematics and Economics 1, 55-72. Falin, G.I. (1989). Periodic queues in heavy traffic. Adu. Appl. Prob. 21, 485-487. Gerber, H.U. (1979). An Introduction to Mathematical Risk Theory. S.S. Huebner Foundation Monographs, University of Pennsylvania, Philadelphia, PA. Graham, A. (1981). Kronecker Products and Matrix Calculus with Applications. Ellis Horwood, Chichester. Grandell, J. (1977). A class of approximations of ruin probabilities. Stand. Actuarial J., 37-52. Harrison, J.M. and A.J. Lemoine (1977). Limit theorems for periodic queues. J. Appl. Probab. 14, 566-576. Janssen, J. (1980). Some transient results on the M/SM/l special semi-Markov model in risk and queueing theories. ASTIN Bull. 1141-51. Janssen, J. and J.M. Reinhard (1985). Probabilit6s de ruine pour une classe de modeles de risque semi-Markoviens, ASTIN Bull. 15, 123-133. Johnson, M.A. and M.R. Taaffe (1989). Matching moments to phase distributions: mixtures of Erlang distributions of common order. Stochastic Models 6, 2.59-281. Lemoine, A.J. (1981). On queues with periodic Poisson input. J. Appl. Prob. 18, 889-900. Lemoine, A.J. (1989). Waiting time and workload in queues with periodic Poisson input. J Appl. Probab. 26, 390-397. Lucantoni, D.M. (1990). New results on the single server queue with a batch Markovian arrival process. Stoch. Models 6. Moler, C. and C. Van Loan (1978). Nineteen dubious ways to compute the exponential of a matrix. SIAM Review 20, 801-836. Neuts, M.F. (1981). Matrix-Geometric Solutions in Stochastic Models. Johns Hopkins University Press, Baltimore, MD. Neuts, M.F. (1989). Structured Stochastic Matrices of the M/G/I type and their Applications. Marcel Dekker, New York.
274
S. Asmussen, T. Rolski / Computational methods in risk theory
Press, W.H., B.P. Flannery, S.A. Teukalsky and W.T. Valtering (1988). Numerical Recipes. The Art of Scientific Computing. Cambridge University Press, Cambridge. Ramaswami, V. (1986). Nonlinear matrix equations in applied probability. SIAM Review 30, 256-263. Ramaswami, V. (1990). From the matrix-geometric to the matrix-exponential, Queueing Systems. Regterschot, G.J.K. and J.H.A. de Smit (1986). The queue M/G/l with Markov-modulated arrivals and services. Math. Oper. Res. 11, 456-483. Reinhard, J.M. (1984). On a class of semi-Markov risk models obtained as classical risk models in a Markovian environment. ASTIN Bull. XIV, 23-43. Rolski, T. (1981). Queues with non-stationary input stream: Ross’s conjecture. Adu. Appl. Prob. 13, 603-618. Rolski, T. (19871. Approximation of periodic queues. J. Appl. Probab. 19, 691-707. Rolski, T. (19891. Relationships between characteristics in periodic Poisson queues. Queueing Systems 4, 17-26. Rolski, T. (19901. Queues with nonstationary inputs. Queueing Systems 5, 113-130. Sengupta, B. (1989). Markov processes whose steady-state distribution is matrix-exponential with an application to the GZ/PH/l queue Adu. Appl. Probab. 21, 159-180. Sengupta, B. (1990). The semi-Markovian queue: Theory and applications. Stochastic Models. Thorin, 0. (19701. Some remarks on the ruin problem in case the epochs of claims form a renewal process. Stand. Actuarial. J., 29-50. Thorin, 0. (19711. Further remarks on the ruin problem in case the epochs of claims for a renewal process. Stand. Actuarial. J. 14-38, 121-142. Thorin, 0. (19741. On the asymptotic behaviour of the ruin probability for an infinite period when epochs of claims form a renewal process. Stand. Actuarial. .I., 81-99. Thorin, 0. (1975). Stationarity aspects of the Sparre Andersen risk process and the corresponding ruin probabilities. Stand. Actuarial. .I., 87-98. Thorin, 0. and N. Wikstad (1973). Numerical evaluation of ruin probabilities. ASTIN Bull. VII, 137-153. Thorin, 0. and N. Wikstad (1974). Numerical evaluation of ruin probabilities when the claim distribution is lognormal. ASTIN Bull. VIII, 231-245.
and Economics
259
10 (1991) 2.59-274
Computational methods in risk theory: A matrix-algorithmic approach Soren Asmussen Institute of Electronic Systems, Aalborg Unicersity, DK-9220 Aalborg, Denmark
Tomasz Rolski * Mathematical Institute, Wroclaw University, 50-370 Wroctaw, Poland Received
April
1991
Abstract: This paper is concerned with the numerical computation of the probability t)(u) of ruin with initial reserve u. The basic assumption states that the claim size distribution is phase-type in the sense of Neuts. The models considered are: the classical compound Poisson risk process, the Sparre Anderse process and varying environments which are either governed by a Markov process or exhibit periodic fluctuations. The computational steps involve the iterative solution of a non-linear matrix equation numerical examples are presented. Q = cY(Q) as well as the evaluation of matrix-exponentials e QU.A number of worked-out Keywords: Ruin probability, ment. Periodic environment.
Phase-type
distribution,
Matriu-algorithmic
methods,
Non-linear
matrix
iteration,
Markovian
environ-
1. Introduction
This paper is concerned with the numerical computation
of infinite horizon ruin probabilities
(1-l) where (S,} is the claim surplus process, u is the initial reserve and T.(U)= inf{t L 0: S, > u} is the time of ruin (first passage time to level u). Thus N
s, =
c
l_J-pt,
i=l
where IV, is the number of claims up to time t, the Ui are claim sizes and p is the premium rate. For example, in the classical compound Poisson model, {IV,}is a homogenous Poisson process and the Ui are independent of {IV,) and i.i.d. with common distribution B; even in this simple case, the numerical computation is non-trivial. At the core of our discussion is the fact that risk theory and queueing theory are mathematically closely related areas. For the compound Poisson model with p = 1 the relation is t/l(u) =
qw>u>,
(1.2)
where W is the stationary waiting time in the M/G/l queue with the same arrival process and service times distributed as the claim sizes. The quantities asked for in queueing applications are, however, * The second
author
0167-6687/92/$05.00
was supported
in part by The Philip Lundberg
0 1992 - Elsevier
Science
Publishers
Foundation
B.V. All rights reserved
S. Asmussen,
260
T. Rolski
/ Computational
methods
in risk theory
typically more simple than the whole distribution of W as in (1.2). Most often only the mean (El+’or the first few moments are required, and in the M/G/l case these are explicit in terms of moments of U. The GI/M/l queue is one more example, where explicit solutions are available, but more general systems present difficulties. In queueing theory, the development of computational methods has, however, been rapidly progressing during the last one or two decades, due not least to the application potential in areas like computer networks. Of particular importance for our discussion is the matrix-algorithmical method initiated by Neuts [see Neuts (1981, 1989)], and the extensive bibliographies therein. Briefly speaking the purpose of this paper is to explain that this approach is relevant in risk theory as well, to carry out the necessary adaptation, and to present some worked-out numerical examples. Much of the literature on ruin probabilities deals with approximations, for example the celebratedCram&-Lundberg approximation
(1.3) Our concern is here, however, the computation of the exact values. The basis assumption, which makes this possible, is restricting the claims to be phase-type; the crucial fact is that phase-type assumptions on the claims again provide a phase-type structure of the ruin probabilities. The phase-type set-up covers for example finite mixtures of exponential distributions and Gamma distributions with integer parameter. Whereas many papers in the literature deal with such examples, the general class of phase-type distributions have hardly been exploited at all in risk theory. Essentially, the distributions in the literature which do not fit into the framework are heavy-tailed ones like the lognormal [see Embrechts and Veraverbeke (1982)]. We start the paper by a brief summary of phase-type distributions in Section 2. Section 3 then deals with the classical risk model with Poisson arrivals of claims. For this case, the relevant expression for e(u) is a translation of a queueing result given in Neuts (1981) and, despite the simplicity of the solution and that it is well-known in queueing theory, it covers many special cases treated separately in the literature on ruin probabilities. The rest of the paper then relies on some more recent advances in algorithmical queueing theory, see in particular Asmussen (1991, 1992) Sengupta (1989, 1990), Ramaswami (1990) and Lucantoni (1990). The essential feature is that whereas largely the original matrix-geometric method as surveyed in Neuts (1981) is best suited for studying discrete processes (like the numbers of customers present in the queue), the papers just listed are able to deal directly with continuous-state problems. In Section 3, we study the risk process with renewal arrival of claims (sometimes called the Sparre Andersen process); the material is a translation of Asmussen (1992) which in turn has some features similar to Sengupta (1989). Section 4 deals with risk processes in a Markovian environment as studied by Asmussen (1989); compared to that paper which deals largely with approximations, the class of models for which we can present an explicit (or at least computationally tractable) exact solution is considerably widened. Finally, in Section 5 we give some remarks on models with a periodically varying Poisson intensity and how they can be approximated by a Markovian environment model. We remark at the end, that from the point of risk theory it may not feel all that satisfying to derive the results and algorithms by translation of the queueing literature. However, at least within the framework of this paper, we see both areas as special cases of random walks and their extensions (e.g. Markov additive processes); the approach that we present just happened to be initiated in queueing theory rather than in risk theory. A more self-contained (and therefore also more lengthy) exposition will be given elsewhere [Asmussen (1992)l. In the rest of the paper we scale for simplicity the claim surplus process to have premium rate p equal 1 and suppose that the safety loading 77 is positive; this implies that the surplus process tends to --OOwith probability 1. The traffic intensity p of queueing theory is related to 77 by p = l/(1
+ 7~)
or
77 = (l/p)
- 1.
(1.4)
S. Asmussen, T. Rolski / Computational methods in risk theory
2. Phase-type
261
distributions
That the distribution B of the generic claim U is of phase-type with representation that the d x d matrix T is the restriction to 11,. . . , d} of the intensity matrix
(r, T, d) means
of a Markov jump process on IO,. . . , d) with 0 as absorbing state, and that U is distributed as the time to absorption in 0 given the initial distribution r (written as row vector in the following). In particular, the vector t, of exit rates is t, = - Te where e is the column vector of ones. We assume that T is a proper subintensity matrix. For example a hyperexponential density (i.e. a finite mixture of exponential densities) Cfrjai epalx is phase-type with T = (a,, . . ., adjdiae, and an Erlang density adxd-‘/Cd l)! e- uX is phase-type with representation I-,
(Y --(Y
0 p=(l
Fundamental 1 -B(x)
0
..a
o),
. .
T=
0 G? ‘.
.
... ...
0 0
0 0
.
(2.1)
(j
(j
0
..:
_La
0
0
0
*.*
0
b(x)
=B’(x)
CY
-a
analytical identities are =lmB(dy) X
=reTxe,
=reTXtO,
(2.2)
T)-‘to,
(2.3)
-1
e=r(
pL(g)= on”B(dn) /
= ( - l)“n!rT-“e.
-sZ-
(2.4)
See Neuts (1981) or Asmussen (1987, Ch. III 6) for more detail. A crucial fact of these formulae is their amenability to numerical computations: all that is needed is (a) matrix inverses [the standard algorithm is Gauss-Jordan elimination with full pivoting, cf. Press et al. (198811 and (b) matrix exponentials for which a variety of computational recipes are available as well, see the appendix. A further important property of phase-type distributions is denseness: any given distribution on (0, w) can be approximated arbitrarily close by a phase-type distribution. This point is illustrated in Figures 1 and 2, where we have plotted phase-type fits to two truncated lognormal distributions, more precisely the distributions of eUU-” I2 A 2 where U is standard normal. These examples are the ones to be used in the rest of the paper, and thus some further comments are in order. Example 2.1 (Figure 1). The choice u = 1.8 in Figure 1 is from Benckert and Jung (1974) who found this to provide an adequate fit to some Swedish fire insurance data, and the distribution appears again in Thorin and Wikstad (1973); the truncation at 2 could correspond, e.g., to reinsurance. The parameters of the fitted phase-type distribution B, (say) in Figure 1 are d = 2,
.H = ( 0.5614
0.4386))
T=
- 8.640 0.101
(2.5)
The mean is the same as for the truncated log-normal, viz. 0.6015 (note that the unrestricted log-normal distribution has mean 1 and thus the truncation matters more than a brief look at Figure 1 might
262
S. Asmussen, T Rolski / Computational methods in risk theory
i
3
-
- -
- -
lognormal;
(T =
1.8
Bl
2
1 :\\
-4
--__
I
-----mm-__
I 2
1 Fig. 1.
indicate). The squared coefficient of variation (s.c.v.) is 1.8988. The standard comparison is the exponential distribution, where the S.C.V.is 1, and since a rough rule of thumb asserts that the risk increases with any added stochastic variability, we thus expect larger ruin probabilities than for exponentially distributed claim sizes with same mean.
Example 2.2 (Figure 2). The choice u = 0.8 in Figure 2 was made to obtain a distribution with a different shape than in Figure 1, but is otherwise arbitrary. The parameters of the distribution B, on Figure 2 are d = 4, m= (0.9731 T=
i
0.0152
- 28.648 0.102 0.113 0.100
0.0106
28.532 - 8.255 0.107 0.102
0.00 1
0.089 8.063 - 5.807 0.111
(2.6)
3-
--
2-
-
- -
lognormal;
0 = 0.8
B2
1
Fig. 2.
2
263
S. Asmusen, T. Rolski / Computational methods in risk theory
The mean of l3, is 0.8849 and the S.C.V.is 0.5132, and thus we can expect smaller for the exponential case.
ruin probabilities
than
Some basic references for fitting phase-type distributions are Bux and Herzog (1977), Johnson and Taaffe (1989) and Bobbio and Cumani (1990); the distributions B,, B, were obtained by the method of Asmussen and Nerman (1990). It should be noted, however, that our aims here are not - to survey just a few _ to discuss - to neglect distribution.
statistical aspects of phase-type fitting, only to point out that phase-type distributions with phases may catch the shape of a given (reasonably smooth) claim size distribution very well; whether phase-type distributions are appropriate for fitting heavy-tailed distributions; the problem that ruin probabilities are sensitive to even small changes in the claim size
Nevertheless, we may argue that if one wants to distinguish between lognormal in the figures, one needs a substantial amount of statistical data, presumably most other parametric models would fail a goodness of fit test anyway. With to explain how to calculate ruin probabilities using phase-type assumptions, purposes in the rest of the paper through (2.5), (2.6). For our computational the following structure of the ruin probabilites 4(u) will be found for a variety claims:
or phase-type distribution so much that these and these remarks, we proceed exemplified for numerical approach it is crucial that of models with phase-type
(u, S, d) if for some d-dimensional Definition 2.1. A function I,%(U)is of phase-type with representation subprobability vector u and some d x d subintensity matrix S we have I++(U)= v esue.
3. The classical compound
Poisson model
For the classical compound Poisson model with claim size distribution B of phase-type and arrival intensity p, the corresponding queue is M/PH/l and was considered by Neuts (1981, pp. 52-57). In risk theoretical terms the result is the following: Proposition
3.1.
T, d),
the ruin probability is phase-type with representation CT+, Q, d), where and
Q=T+tp+,
T+= -@rT-‘.
Proof. This follows directly by combining follows. First note that by standard matrix -rT-r
=
mrer* /0
(3.1) (1.2) and Neuts calculus
(1981). A short self-contained
dx
proof goes as
(3.2)
[the convergence in (3.2) is ensured by the Perron-Frobenius theorem because T is a proper subintensity matrix] so that the components of rr+ are nonnegative and the total mass IIr+ II = r+e of r+ is I(rr+ 1)= IIpuma
eTx dx II = p&)
Next, by the Pollaczek-Khinchine P(W4u) where
pB, P
= (1 -p)
PB
formula
[see e.g. Asmussen
(1987)1, (3.3)
5 p”Bo*“(u), n=O
is the (defective)
1 -B(x) (1)
= p < 1.
distribution
= /?P eTxe = p( -rT-‘)
with density e”(
- Te) = r+ eTXtO.
S. Asmussen, T. Rolski / Computational methods in risk theory
264
_______.
B
1
. . . . . . . . . . . . . . . . . . . . exponential
25’EU
WEU
Fig. 3.
the That is, pB, is defective phase-type with representation (P,, T, d). In (3.3) we can thus represent distribution of W by a geometric-compound of a (P,, T, d) phase-type distribution which is phase-type q (a,, T + c,,r+, d) as is readily seen by a probabilistic argument. Example 3.1. We shall reinspect an example in Gerber arrival intensity is p = 3 and the claim size density is
(1979) in view of Proposition
3.1. Here
the
b(~)=+.3e-~~++.7e-~~. Thus
B is hyperexponential
r+=
The eigenvalues
Q(U)
=
in accordance
with r = (+ i), T = (-3
(3
eQ"e =
.
Q=T+t,,r+=
&),
of Q are - 1 and
r+
- 7)diag, so that
- 6, and standard
diagonalization
techniques
yield
$ edu + & ep6’
with Gerber
(1979).
For numerical illustration, we fix from now on the safety loading 77 to 10% (so that p = 0.9091) and consider $(u) in the range 0 I u _<50 p, where p is the average claim size (for this simple model, we have just p = [EU). The ruin probability @l(u) is plotted on a logarithmic scale in Figure 3 for Poisson 2.1, 2.2 and the arrivals and three different claim size distributions, namely B,, B, of Examples exponential distribution with mean one (note that the arrival rate /? is determined by either of n, p sot that for B, we get p = p/p = 1.496 and for B, we get p = 1.017). The basic parameters of (3.1) were calculated as d = 2, ?r+= ( 0.10757
0.80152))
Q = ( - ;:;;;;
_;‘;;;;)
S. Asmussen, T. Rolski / Computational methods in risk theory
265
for Example 2.1, and as d = 4, r+=
(0.038
0.143
( -28.648
Q= \
0.102 0.124 0.171
0.210
28 532 - 8.254 0.149 0.368
0.518), 0.089 8.064 - 5.746 0.503
0.027 0.088 5.446 - 1.210
for Example 2.2. It is seen from Figure 3 that our expectations of which claim size distribution has the higher risk were confirmed. The almost linear appearence of the curves on Figure 3 and in later figures of the paper reflect the fact that the ruin probabilities are close to being exponential (the scale is logaritmical). There are essentially two reasons for this. The Cramer-Lundberg approximation (1.3) implies that the curves have linear asymptotes; thus possible deviations from linearity have to be sought for small values of U. However, also here the curves look almost linear, which we would explain by the fact that n = 10% is rather small so that we are in the domain of the diffusion approximation, cf. Grandell (1977).
4. The Sparre Andersen process In this section we suppose that the number of claims up to time t is a renewal process with interarrival distribution A and that the claim sizes are i.i.d. and independent of IN,] with common distribution B of phase-type with representation (w, T, d). The renewal assumption was introduced by Sparre Andersen (1957) and extensively studied by Thorin (1970, 1971, 1974, 1975) and Thorin and Wikstad (1973, 1974). We consider here two cases. In the first one we suppose that the process starts at a claim arrival (the arrivals form a non-delayed renewal process) and in the second case [Thorin (1975)] that the arrival process is a stationary renewal process. The ruin probability for the non-delayed case is denoted by I/J(U) and for the stationary case by 4*(u). The corresponding queue is GI/PH/l and was studied by Asmussen (1992) [see also Sengupta (1989)]. Before we recall the result from Asmussen (1992), we introduce some notation. For a subintensity matrix K (this provides the convergence), define AIK] = 10”eKXA(dx), and let V, be the nonlinear matrix operator defined by V,(K)
= T-t&&K],
(4.1)
where t, = - Te. Proposition 4.1. For the renewal model with claim size distribution of phase-type, the rain probability $(u) is of phase-type CT+, Q, d), where the subintensity matrix Q solves the non-linear matrix equation
Q = qd’l(e)
(4.2)
and rr+ = e’(Q - T)/e’t,; the ruin probability I/I*(u> is of phase-type (p y, Q, d), where y is the distribution soluing y(T + t,lr) = 0, ye = 1. Proof. The result concerning $(u) follows from Asmussen (1992) and the common interpretation of $(u) and p(W > U) as P(M > U) where M is the maximum of the random walk with generic increments distributed as the independent difference between a claim size (service time) and an interarrival time. For ICI*(u), one notes similarly that e*(u) = P(V> u), where V is the virtual waiting time of the GI/PH/l queue. 0
It is important
to notice that Q may be computed
Q n+l = !P,(Q,>, (n = 0, 1, . . .I.
as the limit Q = lim, +_Q,,, where Q, = T and
Two interarrival distributions, A,, A, were considered, one phase-type and one deterministic, respectively. We shall explain next how to carry out the crucial steps of calculating AIK] in (4.1). Consider first the phase-type case. Then (see the appendix for the Kronecker notation):
266
S. Asmussen, T. Rolski / Computational methods in risk theory
Proposition
d[K]
4.2.
Zf A is phase-type with representation (Y, S, q), then
= (z@v)(-K6Ls~‘(zc3s”),
(4.3)
where sg = -Se. Proof.
We have eKXveS”s,=(Z~Y)(eKx~esx)(Z~Sg)
(4.4)
and hence using eKx @ eSXe(“*‘jX we get
a[K]
=
irn eKxv eSXsOdx
= (ZC3v)im(eKr
@esX)
dx(Z@ss,)
= (Z@pY)( -KfBS)-‘(z@ss,).
0
(4.5)
In the deterministic case, where A, is degenerate at l/p (say), we have AlK] = eKla; to compute the matrix exponential, we used uniformization, see the appendix. We performed computations for the claim size distribution B, and the following interarrival distributions: A, - the hyperexponential distribution with density 0.4 .5c~ e-5a + 0.6. (yepa; A, - the distribution degenerate at l/p; and, for the sake of comparison, also the exponential distribution (which is already included in Figure 3). The scaling factors a, p were determined by 77= 10%. The illustration is in Figure 4. Again, the comparison of the risks involved in the different cases correspond nicely to the fact that the hyperexponential distribution has a S.C.V. which is higher than for the exponential distribution which is in turn higher than for the deterministic case. The basic parameters of Proposition 4.1 are d = 2, r+=
( 0.159
for the hyperexponential 7r+= (0.013
Q= ( -;$“6
0.780))
inter-arrival
-
t I
times, whereas for the deterministic
Q= ( -;$;
0.814),
0.001
_;‘i;;)
-
I
-
I
_;‘;;;).
-
I
case we got
-
I
-
-
I
25’EU Fig. 4.
__-_____
B
....................
exponential
-
I
-
I
-
I
1
_
-
I
50’EU
267
S. Asmussen, T. Rolski / Computational methods in risk theory
5. Markov modulated
arrivals
In this section we suppose that the environment is governed by an irreducible Markov process with states i= l,..., p and intensity matrix A = {hii}, say. This means that if the governing process (environment) is in state i, then claims are Poisson with intensity pi and claim sizes U”’ have common distribution C”’ of phase-type with representation (rci), Sci), 9”‘) (i = 1,. . . , p). For such a process the average arrival rate is p = Zip_, yipi, where y is the probability solution of the stationarity equation r/i = 0, the average claim size is
(5.1) and the average amount of claims per unit time is p = pp. Thus, as before, (1 -p)/p, and letting pi = p,lEU”‘, ni = l/p, - 1, we can rewrite (5.1) as P =
5
YiPi
,)./ =
i=l
5
ESi.
i=l
P
the safety loading is
(5.2)
This formula expresses overall characteristics of the Markov-modulated process in terms of the p standard compound Poisson risk processes that we obtain by fixing the parameters of the ith to be pi, Cc’). Note, however, that whereas we need to assume that the overall safety loading n is > 0, it is possible that vi < 0 for one or more i; such states i represent environmental conditions which are unfavourable in the sense that when present, a loss is suffered. The model has been studied in risk theory by Janssen and Reinhard (1985), Janssen (1980), Reinhard (1984) and Asmussen (1989), whereas some main early queueing references are Burman and Smith (1986) and Regterschot and de Smit (1986). We follow here Asmussen (1991). Note that the connection between the queue and the risk process is slightly more complicated than in (1.2) by involving a time-reversion of A; however, the basic fact is that both the waiting distribution for the queue and the ruin probabilities can be expressed in terms of maxima of a certain continuous-time random walk with Markov-dependent increments. We denote by J/i(~) the ruin probability subject to the initial condition that the reserve is u and the environment is i(i= l,..., p). Proposition 5.1. For a risk process in a Markovian environment with phase-type claim size distributions, the ruin probability *i(U) is phase type with representation CT?), Q, E) where: p, j’=l,...,
E={(j,j’):j=l,...,
q(j)},
and p~j;=Bj(el~rr(j))(-K~S(j))-‘(ej~Z), where the K is the solution of the nonlinear matrix equation K+A-(P)diag+
W(K).
(5.3)
Here W(K) denotes the matrix with the i&h element
The result is proved in Asmussen (1991) with W(K) = lFeKXs(dX), where S(dx) denotes integration of the jth column with respect to /3jC”‘(dx); the further calculations are as in Proposition 4.2. q Proof.
268
S. Asmussen,
T. Rolski / Computational
methods in risk theory ecn = 24 . . . . . . . . . . . . . . . . . . . . ecn = 16 I
I
O.Ol-------__-___
0.001
-
-
-
I
-
-
-
I
I
-
I
-
I 25’EU
-
I
-
I
_
I
_
_
I 50’EU
Fig. 5.
We performed computations for the following (thus 4 (I) = 2 , q(*) = 4) and the intensity matrix
data: p = 2, PI = 3.325, & = 0.781, C(r) = B,, C’*’ = B,
The value of & is chosen to obtain 7, = - 50% (or equivalently pr = p, pr = 2) and p2 is then determined from n = 10% by equation (5.2). The role of the parameter A is as a measure of the degree of Markov-modulation. As A + co, the Markov-modulated process tends to the homogeneous Poissonian claim arrival stream [this was noted in Asmussen (1989) and proved formally in Rolski (199O)l whereas for A + 0, the initial state becomes predominant. The three curves of the ruin function in Figure 5 correspond to three different values of h, which in turn were determined by letting the mean number ecn = - C?= ryipiAz; ’ of claims till completion of an environmental cycle be ecn = 4, 16, 64, respectively. Thus, the case of eta = 64 should correspond to an almost homogeneous Poissonian input, whereas the case ecn = 4 should be almost equivalent to the classical compound Poisson model with parameters given by the initial condition. The basic parameters of Proposition 5.1 computed for case ecn = 4 were K= P:)=
P$?=
Q=
0.511 2.554
-0.246 1.229 1 ’
( TT$~, : j = 1, 2, j’ = 1,. . . , q(j)) = (0.173
0.575
0.002
0.018
0.042
0.142),
(0.013
0.238
0.029
0.105
0.151
0.365 ),
1 -7.489 0.273 0.000 0.000 0.004 \ 0.024
5.813 - 0.525 0.000 0.001 0.069 0.443
0.012 0.012 - 28.648 0.102 0.121 0.153
0.122 0.018 28 532 - 8.255 0.137 0.297
0.278 0.042 0.089 8.064 - 5.763 0.393
0.945 0.141 0.027 0.087 5.402 - 1.495
269
S. Asmussen, T. Rolski / Computational methods in risk theory
The ordering environmental
of the states in E is lexicographical, so that e.g. the fifth element state 2 and phase 3 of B,. For ecn = 16 we have
K= ( 0.318 1.598
0.709
0.001
0.006
0.015
0.056),
P?‘= (0.009
0.211
0.029
0.107
0.157
0.383),
’ -7.380
Q=
0.289
6.709 - 0.390
0.000 0.003 0.018
0.001 0.061 0.392
0.000
\
to
-0.070 0.347 1’
(0.190
4!)=
corresponds
0 .ooo
0.004 0.001 - 28.648 0.102 0.121 0.154
0.042 0.006 28.532 - 8.255 0.138 0.301
0.375 \ 0.056 0.027 , 0.088 5.407 - 1.464 /
0.100 0.015 0.089 8.064 - 5.761 0.403
and for ecn = 64, K=
- 1.317 0.263
rY)= (0.195
0.774
0.000
0.002
0.004
0.017))
@=
(0.008
0.193
0.029
0.108
0.159
0.391),
Q=
l-7.342 0.295 0 .ooo 0 .ooo 0.002 \ 0.015
7.138 - 0.326 0.000 0.001 0.057 0.368
0.012 0.002 28.532 - 8.255 0.138 0.303
0.001 0.000 - 28.648 0.102 0.121 0.154
0.029 0.004 0.089 8.064 - 5.761 0.407
0.112 ’ 0.017 0.027 0.088 . 5.410 -1.4491
In Figure 5, the initial condition is state 1 of the environment, whereas in Figure 6, we have plotted @r(u) as well as I/J&U) for the intermediate case ecn = 16. Note that Figure 6 gives an example where the non-exponentiality of the ruin probabilities for small u is apparent; the version of the Cram&-Lundberg
ecn =I6 .. . . . . . . . . . . . . . . . . . .
0.1
-
____-------
0.01
-
___--------
0.00,
_
_
_
I
-
I
-
I
-
-
I
I 25’EU Fig. 6.
-
-
I
-
I
-
I
y
-
-
I 50’EU
270
S. Asmussen, T. Rolski / Computational methods in risk theory
approximation given in Asmussen (1989) implies that the asymptotic linear slope of the two curves is the same. Figure 5 reconfirms the findings of Asmussen (1989) (and much of the literature on similar models): the risk is considerably increased as the environment becomes more variable.
6. Periodic environment The claim arrival intensity may depend on some seasonal effect, in which case we can model the stream of claims by a periodic Poisson process. We scale for simplicity the model to have period 1 and assume that the claim sizes (Ui} are i.i.d. with common distribution B. Then, if we ask for the ruin probability starting at a given time of season S, N, = N”(t) is a Poisson process with a periodic arrival rate p’(s + .). In this case we denote the ruin probability by I/J”(U), with the superscript s E [O, 1) indicating the time of the initial season. It can be written as V(u)
=l+;(
;+‘i
)Bj,
(6.1)
where for s E (0, 11, {M”(t)} is the nonhomogeneous Poisson process with intensity function p’(s + t). Note that for such a process the average arrival rate is /? = l@‘(t) dt. The model has been studied, for example, in Dassios and Embrechts (1989) and Asmussen and Rolski (1991). In the queueing literature, a few papers are devoted to queues with periodic input stream; see Harrison and Lemoine (1977), Lemoine (1981, 1989), Afanas’eva (1985), Rolski (1981, 1987, 1989) and Asmussen and Thorisson (1987). A characteristic of the theory of periodic queues theory for a periodic arrival rate function /3*(t) is the stationary workload I”“’ at given time s E [O, 1) of the season which, following Lemoine (1981), has representation M,(I) VS)
= sl’:
c
i where {M,(t)] /3’(t) =P*(-t)
)
q-t
i=l
(6.2)
1
is the nonhomogeneous we see that
Poisson
process
with intensity
l/V(U) = P{v”-“‘>U}.
function
/3*(s - t). Thus
setting
(6.3)
Rolski (1987) proposed to approximate queueing characteristics of periodic queues by Markov modulated ones. In risk theoretical terms his result translates via (6.3) to the following. For each rr = 1, 2, . . . define the nth approximating risk process as the Markov modulated risk process with the environment governed by a Markov process having states 0,. . . , n - 1 and intensity matrix A given by ‘---II 0 \ n
...
n -n . .
0 n . .
**. . . . .
()
(j
.::
0 ’ 0 (6.4) -n,
Arrivals occur at rate pi = p’((i each approximating risk process function for the corresponding continuous intensity functions /I’ Proposition
6.1.
Zf i/n
1)/n) if the governing Markov process is in state i. The claim sizes at are the same as in the original risk process. Let +,!“‘(u) be the risk nth risk process. The following holds for continuous or piecewise [see Rolski (1987)]:
+ s E [O, l), then I/J:“‘(U) -+ I/J”(~).
Note that it seems intuitively plausible that the above result can be extended to the more general model with claim size distribution B”(t) depending also on the season s E [O, 1). For a proof, see Asmussen and Rolski (1991).
S. Asmussen, T. Rolski / Computational methods in risk theory
271
....................
n =
-_-_-mm.
n= 16
8
_-----
0.00,
_
_
_
I
-
I
-
I
-
-
I
I 25’EU
-
-
I
-
-
I
I
-
-
I 50’EU
Fig. 7.
We performed p’(t)
computations
for the intensity function
= p( 1 + sin 27rt)
(6.5)
and exponentially distributed claim sizes with mean claim size l/6, say. The risk function I/#‘(U) for the approximating Markov modulated models is plotted in Figure 7 for n = 8, 16, 24. Here we chose p = p and S = 1 to have safety loading n = 10%. To demonstrate seasonal variability, we plotted in Figure 8 the function 4”(O), (0 IS < 1) approximated by I@~‘(O) (i = 0,. . , , 23), taking p = 2 or 10 and S determined by condition that n = 10%. The greater variability for the case p = 10 is explained by the fact that p is the average number of arriving claims per unit interval and that increasing fi makes the claim sizes smaller but more numerous. Thus, the stochastic variability in the claim sizes and the claim arrival times matters less compared to the seasonal variability of P’(t). Because of the general properties of the Markov-modulated model, it is not surprising that the curves again are almost exponential. However, under some regularity conditions, the periodic model also allows a version of the Cramer-Lundberg approximation. It is even possible to find adjustment coefficients y,, for the Markov modulated models as well as y for the periodic case. More precisely, it can be seen from
1.o m II0 l
. . .
0.9
.O
~!a.
. l *
0.6
Fig. 8.
p=2 p=10
272
S. Asmussen, T. Rolski / Computational methods in risk theory
Asmussen (1989) that for the nth approximating n
I + Yn+ (&Y,)
- l)P’((i
- 1)/n)
model ‘y, is the positive root of the equation = I
n
i=l ni
7
(6.6)
1 [here g(a) = /r e”“B(dx)] and y is the positive root of the equation Y =P(&,) - I)> (6.7) that is, the same as for the corresponding standard compound Poisson model with constant arrival intensity p. From (6.6) and (6.7) it may concluded that yn + y and moreover one can gain some information concerning of the speed of convergence. For details of these facts, see Asmussen and Rolski (1991).
7. Appendix We discuss here the basic procedures
used in our computations.
1. The exponential of a matrix.
There are various methods available for computing of matrix exponentials eK, see for a review Moler and Van Loan (1978). If K is a subintensity matrix, one further alternative not treated in Moler and Van Loan (1978) is popular in applied probability, namely uniformization (the procedure is similar in spirit to the method called scaling and squaring, which was recently used in Sengupta (1989) for similar purposes). For a subintensity matrix K = (kij), let n = 11K II = max,,j I kij I. Then IIZ + 7-lK II I 1 and
iu
,,,k
=c k=O
k!
e-"(I+
T-‘K)~
+Error,
(7.1)
where IIError II 5 C~=i,+l(~k/k!) e-7. This method was used to compute AIK] = eK/P for the Sparre Andersen process with deterministic arrivals, cf. Section 4. However, for the ruin probability function +(u) the numerical computations in this paper are based on the solution of a linear system of differential equations, for which we used the standard fourth-order Runge-Kutta method [see e.g. Press et al. (1988)]. More precisely, we need to evaluate functions of the form r eKth, where r is a row vector and h is a column vector. We note that the vector z(t) = eK’h fulfills the linear system of differential equations z’ = Kz, with the initial condition z(O) = h. This procedure appears particularly efficient whenever we need to compute the function r eKth = PZ for several values of t, whereas the before mentioned method is more convenient for computing a single or a few values. 2. Kronecker calculus. We recall that if A(‘) is a k, X m, and A’*’ a k, X m2 matrix, then the Kronecker Further(tensor) product A”’ C%J A’*’ is the (k, x k,) x (m, x m,) matrix with (i,i,X j, j,)th entry a$,a$. more, if k, = m, and k, = m2, then the Kronecker sum is defined by
A(‘) @A(*) = A(‘) @Zk*+ Zk, C3A’*‘.
(7.2)
See, e.g., Graham (1981). The following formula is basic for this paper: Suppose that A and B are both square such that CY+ /? # 0 whenever LYis an eigenvalue of A Lemma 7.1. and p is an eigenvalue of B. Let further 7r, v be any row vectors and h, k any column vectors. Then /0
XsreA’h~veB’kdt=(sr8v)(A@B)-‘(eA”BX-Z)(h@k).
(7.3)
S. Asmussen, T. Rolski / Computational methods in risk theory
273
If moreover (Y + p < 0 for all eigenvalues, then mreAthm?tk I0 Proof.
dt=
To prove (7.3) it suffices to differentiate
r eAth . v esrk = (T @v) eAQBx( h 8 k), see Graham (1981).
(7.4)
-(v8v)(A$B)-‘(h&k). both sides and verifjl the identity
P-5)
0
References Afanas’eva, L.L. (1985). On periodic distribution of waiting-time process. Stability Problems for Stochastic Models (Uzhgorod, 1984) l-20. Lecture Notes in Math. 1155. Springer, Berlin. Andersen, E. Sparre (1957). On the collective theory of risk in case of contagion between the claims. Transactions XVth International Congress of Actuaries., New York, II, 219-229. Asmussen, S. (1984). Approximations for the probability of ruin within finite time. Stand. Actuarial. J. 31-57; ibid. (1985) 64. Asmussen, S. (1987). Applied Probability and Queues. Wiley, Chichester. Asmussen, S. (1989). Risk theory in a Markovian environment. Stand. Actuarial J., 69-100. Asmussen, S. (1991). Ladder heights and the Markov-modulated M/G/l queue. Stoch. Proc. Appl. 39, no. 2. Asmussen, S. (1992a). Phase-type representations in random walk and queueing problems. The Annals of Probability 20, forthcoming. Asmussen, S. (1992b). Fundamentals of Ruin Probability Theory. Incomplete book manuscript. Asmussen. S. and 0. Nerman (1990). The EM Algoritm for phase-type distributions viewed as incompletely observed exponential families, in preparation. Preliminary version published in: K. Vest Nielsen, ed., Symposium i Anvendt Statistik (1991). UNI-C, Copenhagen, 335-346. Asmussen, S. and T. Rolski (1991). Risk theory in a periodic environment: The Cram&r-Lundberg approximation and Lundberg’s inequality. Submitted for publication. Asmussen, S. and H. Thorisson (1987). A Markov chain approach to periodic queues. J. Appl. Probab. 24, 215-225. Benckert, L.-G. and J. Jung (1974). Statistical models of claim distributions in fire insurance. ASTIN Bull. VIII, l-25. Bjiirk, T. and J. Grandell (1988). Exponential inequalities for ruin probabilities in the Cox case. Stand. Actuarial .I., 77-111. Bobbio, A. and A. Cumani (1990). ML estimation of the parameters of a PH distribution in triangular canonical form, Technical report. Burman, D.Y. and D.R. Smith (1986). An asymptotic analysis of a queueing system with Markov-modulated arrivals. Opns.Res. 34, 105-119. Bux, W. and U. Herzog (1977). The phase concept: approximation of measured data and performance analysis. In: K.M. Chandy and M. Reiser, eds., Computer Performance. North-Holland, Amsterdam, 23-38. Dassios, A. and P. Embrechts (1989). Martingales and insurance risk. Stochastic Models 5, 181-218. Embrechts, P. and N. Veraverbeke (1982). Estimates of the probability of ruin with special emphasis on the possibility of large claims. Insurance: Mathematics and Economics 1, 55-72. Falin, G.I. (1989). Periodic queues in heavy traffic. Adu. Appl. Prob. 21, 485-487. Gerber, H.U. (1979). An Introduction to Mathematical Risk Theory. S.S. Huebner Foundation Monographs, University of Pennsylvania, Philadelphia, PA. Graham, A. (1981). Kronecker Products and Matrix Calculus with Applications. Ellis Horwood, Chichester. Grandell, J. (1977). A class of approximations of ruin probabilities. Stand. Actuarial J., 37-52. Harrison, J.M. and A.J. Lemoine (1977). Limit theorems for periodic queues. J. Appl. Probab. 14, 566-576. Janssen, J. (1980). Some transient results on the M/SM/l special semi-Markov model in risk and queueing theories. ASTIN Bull. 1141-51. Janssen, J. and J.M. Reinhard (1985). Probabilit6s de ruine pour une classe de modeles de risque semi-Markoviens, ASTIN Bull. 15, 123-133. Johnson, M.A. and M.R. Taaffe (1989). Matching moments to phase distributions: mixtures of Erlang distributions of common order. Stochastic Models 6, 2.59-281. Lemoine, A.J. (1981). On queues with periodic Poisson input. J. Appl. Prob. 18, 889-900. Lemoine, A.J. (1989). Waiting time and workload in queues with periodic Poisson input. J Appl. Probab. 26, 390-397. Lucantoni, D.M. (1990). New results on the single server queue with a batch Markovian arrival process. Stoch. Models 6. Moler, C. and C. Van Loan (1978). Nineteen dubious ways to compute the exponential of a matrix. SIAM Review 20, 801-836. Neuts, M.F. (1981). Matrix-Geometric Solutions in Stochastic Models. Johns Hopkins University Press, Baltimore, MD. Neuts, M.F. (1989). Structured Stochastic Matrices of the M/G/I type and their Applications. Marcel Dekker, New York.
274
S. Asmussen, T. Rolski / Computational methods in risk theory
Press, W.H., B.P. Flannery, S.A. Teukalsky and W.T. Valtering (1988). Numerical Recipes. The Art of Scientific Computing. Cambridge University Press, Cambridge. Ramaswami, V. (1986). Nonlinear matrix equations in applied probability. SIAM Review 30, 256-263. Ramaswami, V. (1990). From the matrix-geometric to the matrix-exponential, Queueing Systems. Regterschot, G.J.K. and J.H.A. de Smit (1986). The queue M/G/l with Markov-modulated arrivals and services. Math. Oper. Res. 11, 456-483. Reinhard, J.M. (1984). On a class of semi-Markov risk models obtained as classical risk models in a Markovian environment. ASTIN Bull. XIV, 23-43. Rolski, T. (1981). Queues with non-stationary input stream: Ross’s conjecture. Adu. Appl. Prob. 13, 603-618. Rolski, T. (19871. Approximation of periodic queues. J. Appl. Probab. 19, 691-707. Rolski, T. (19891. Relationships between characteristics in periodic Poisson queues. Queueing Systems 4, 17-26. Rolski, T. (19901. Queues with nonstationary inputs. Queueing Systems 5, 113-130. Sengupta, B. (1989). Markov processes whose steady-state distribution is matrix-exponential with an application to the GZ/PH/l queue Adu. Appl. Probab. 21, 159-180. Sengupta, B. (1990). The semi-Markovian queue: Theory and applications. Stochastic Models. Thorin, 0. (19701. Some remarks on the ruin problem in case the epochs of claims form a renewal process. Stand. Actuarial. J., 29-50. Thorin, 0. (19711. Further remarks on the ruin problem in case the epochs of claims for a renewal process. Stand. Actuarial. J. 14-38, 121-142. Thorin, 0. (19741. On the asymptotic behaviour of the ruin probability for an infinite period when epochs of claims form a renewal process. Stand. Actuarial. .I., 81-99. Thorin, 0. (1975). Stationarity aspects of the Sparre Andersen risk process and the corresponding ruin probabilities. Stand. Actuarial. .I., 87-98. Thorin, 0. and N. Wikstad (1973). Numerical evaluation of ruin probabilities. ASTIN Bull. VII, 137-153. Thorin, 0. and N. Wikstad (1974). Numerical evaluation of ruin probabilities when the claim distribution is lognormal. ASTIN Bull. VIII, 231-245.