FACULTY OF MATHEMATICS AND COMPUTER SCIENCE PhD Thesis Stochastic models in actuarial science. Theoretical contributions and applications Scienti c… [614330]
UNIVERSITY OF BUCHAREST
FACULTY OF MATHEMATICS AND COMPUTER SCIENCE
PhD Thesis
Stochastic models in actuarial
science. Theoretical contributions
and applications
Scienti c supersvisor:
Prof. Univ. Dr. Preda Vasile
PhD Student: [anonimizat]-Voinea Elena-Grat iela
2016
Dedicated to my missing father, Costin ROBE !
Contents
Published Papers vi
Notations and abbreviations ix
1 Introduction 1
1.1 The motivation and the importance of the research topic . . . . . . . . 1
1.2 Research objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 PhD thesis structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Preliminary notions 6
2.1 Collective models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.1 Univariate case . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.2 Bivariate case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.3 Multivariate case . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Techniques used to evaluate the compound distributions . . . . . . . . 8
2.2.1 Convolutions method . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.2 Panjer recursive formula . . . . . . . . . . . . . . . . . . . . . . 10
2.2.3 Inversion methods using the Fourier transform . . . . . . . . . . 15
2.2.4 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3 Recursions for compound multidimensional models 21
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Recursive evaluation: rst case, corresponding to model B . . . . . . . 23
3.3 Recursive evaluation: second case corresponding to model A . . . . . . 27
3.3.1 Method 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
iv
3.3.2 Method 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4 Alternative methods 37
4.1 Fast Fourier Transform for multivariate claims . . . . . . . . . . . . . . 38
4.1.1 Characteristic function . . . . . . . . . . . . . . . . . . . . . . . 38
4.1.2 Multidimensional discrete Fourier transforms and FFT algorithm 39
4.1.3 Types of errors. Exponential tilting . . . . . . . . . . . . . . . . 40
4.1.4 Numerical study . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.2 Simulation for multivaritate claims . . . . . . . . . . . . . . . . . . . . 54
4.2.1 Comparison between simulation and recursive methods . . . . . 54
4.2.2 Simulation method-Matlab implementation . . . . . . . . . . . . 55
5 Monte Carlo method – a particular study case 59
5.1 Risk analysis basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 Risk management aspects . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3 Risk analysis methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.3.1 Qualitative methods . . . . . . . . . . . . . . . . . . . . . . . . 63
5.3.2 Quantitative methods . . . . . . . . . . . . . . . . . . . . . . . 64
5.4 Monte Carlo simulation for a ship design project . . . . . . . . . . . . . 65
5.4.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.4.2 Simulation 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.4.3 Simulation 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.4.4 Risk Register- Cost and Duration sensitivity . . . . . . . . . . . 73
6 Conclusions 77
6.1 General conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.2 Future research perspectives . . . . . . . . . . . . . . . . . . . . . . . . 78
Bibliography 80
v
Published Papers
I. Published Papers between 2013-2016, in journals ISI Thomson Reuters
indexed from the country or abroad:
1.Elena-Gratiela Robe-Voinea , Raluca Vernic, On a multivariate aggregate
claims model with multivariate Poisson counting distribution, Proceedings of the
Romanian Academy, Series A , ISSN: 1454-9069, accepted. (cotat CNCSIS A, ID
=789, and ISI indexed), 2015-IF=1.658 , accepted
2.Elena-Gratiela Robe-Voinea , Raluca Vernic, Another appproach to the eval-
uation of a certain multivariate compound distribution Analele Universitatii
"Ovidius" Constanta, seria Matematica , accepted, ISSN 1224-1784, E-ISSN 1844-
0835. (cotat CNCSIS A, ID =31, ISI indexed), 2015-IF= 0.383 , accepted.
3.Elena-Gratiela Robe-Voinea , Raluca Vernic, Fast Fourier Transform for mul-
tivariate aggregate claims, Computational and Applied Mathematics, Springer In-
ternational Publishing , 2016, DOI 10.1007/s40314-016-0336-6, Print ISSN 0101-
8205, Online ISSN 1807-0302. (ISI indexed), 2016-IF=0.802 .
4.Elena-Gratiela Robe-Voinea , Raluca Vernic, On the recursive evaluation of
a certain multivariate compound distribution, Acta Mathematicae Applicatae
Sinica, English Series, Springer International Publishing , 2016,. (ISI indexed),
2016-IF=0.250 , accepted.
vi
II. National and international conferences:
1.Elena-Gratiela Robe-Voinea , Raluca Vernic, Fast Fourier Transforms for Bi-
variate aggregate losses. A Matlab application, The 16th Conference of Roma-
nian Society of Statistics and Probabilities, SPSR , April 26, 2013, The Bucharest
University of Economic Studies, Bucharest, Romania.
2.Elena-Gratiela Robe-Voinea , Raluca Vernic, Fast Fourier transform for Mul-
tivariate aggregate losses, The 17th Conference of Romanian Society of Statistics
and Probabilities, SPSR , April 25, 2014, The University of Bucharest, Romania.
3.Elena-Gratiela Robe-Voinea , Exponential tilting for computing multivariate
compound distributions using the Fast Fourier Transform, 22nd Conference on
Applied and Industrial Mathematics – CAIM 2014 , September 18-22, 2014, Uni-
versity of Bac au, Romania, Abstracts, pp. 38.
4.Elena-Gratiela Robe-Voinea , Raluca Vernic, What if di erent types of lcaims
simultaneously a ect an insurance portofolio?, The 18th Conference of Roma-
nian Society of Statistics and Probabilities,SPSR , May 8, 2015, University of
Bucharest, Romania.
5.Elena-Gratiela Robe-Voinea , Raluca Vernic, On the recursive evaluation of
a certain multivariate compound distribution, The 8th Congress of Romanian
Mathematicians , June 28 – July 1, 2015, "Alexandru Ioan Cuza" University of
Iasi, Romania, Abstracts, pp. 120.
6.Elena-Gratiela Robe-Voinea , On the recursive evaluation of a certain multi-
variate compund distribution , Scienti c PhD Students Session of the PhD School
of Mathematics , June 22, 2015, University of Bucharest, Romania.
7.Elena-Gratiela Robe-Voinea , Raluca Vernic, A recursive procedure for a com-
pound risk model with compound-type severity distributions, The 12th Balkan
Conference on Operational Research BALCOR , September 9-13, 2015, "Roma-
nian Naval Academy", Constanta, Romania.
vii
8.Elena-Gratiela Robe-Voinea , Raluca Vernic, Risk analysis based on the Monte
Carlo method for a ship design project, The 19th Conference of Romanian Soci-
ety of Statistics and Probabilities, SPSR , May 27, 2016, Technical University of
Civil engineering, Bucharest, Romania.
9.Elena-Gratiela Robe-Voinea ,Scienti c PhD Students Session of the PhD
School of Mathematics , June 21, 2016, University of Bucharest, Romania.
10.Elena-Gratiela Robe-Voinea , Raluca Vernic, Multivariate aggregate claims
evaluation using the Fast Fourier Transform, 13 eme Colloque Franco-Roumain
de Math ematiques Appliqu ees] , August 25-29, 2016, "Alexandru Ioan Cuza" Uni-
versity of Iasi, Romania.
viii
Notations and abbreviations
r:v: = random variabiles;
d:f: = distribution function;
p:g:f: = probability generating function;
m:g:f: = moment generating function;
c:f: = characteristic function;
c:d:f: = cumulative distribution function;
i:i:d = independent, identically distributed;
g= p.g.f. of a r.v.;
DFT =Discrete Fourier Transform ;
FFT =FastFourier Transform;
fS= the p.f. of S;
fj= the p.f. of Uj;
fL= the p.f.of L;
g= the p.g.f. of a r.v. indexed with the r.v.'s name;
n= (n0;:::;nm);t= (t1;:::;tm);0= (0;:::;0);x= (x1;:::;xm);
Po() = Poisson's distribution with parameter >0;
Bin(t;) = binomial distribution with parameters t2N;2(0;1);
NB( ;) = negative binomial with parameters >0;2(0;1;
MPom+1(; 0;:::;m ) = multinomial Poisson of dimension m+ 1;
Mnomm +1(n;p0:::;pm;) = multinomial of dimension m+ 1;
E[X] = expected value of the r.v. X;
ix
Chapter 1
Introduction
1.1 The motivation and the importance of the re-
search topic
The original results obtained in this research refer rst of all to how to study and model
the aggregate claims distributions of a portfolio insurance. The importance of aggregate
claims distributions is related to the main objectives of actuarial science: calculations
of premiums, evaluation of the needed reserves in order to avoid ruin probability, study
of the e ect of changing the deductibles or of the composition of the portfolio, study of
possible reinsurance covers, etc. There are two main models used for aggregate claims:
the individual and the collective models, see the books Klugman [49], Zbaganu [117],
[118]. To de ne the individual model, [84] assumes that the portfolio consists of a
xed number of independent policies, which implies that the aggregate claims of the
portfolio becomes the sum of the aggregate claim of each of the policies, hence the
claims of the portfolio are related to the individual policies. It is important to mention
that the individual model has been studied mainly by De Pril [16],[17], Dufresne [22]
Dhaene and Vandebroek [20]. In the same manner, [84] arms that in regarding the
collective model, we do not consider the individual policies, instead we must consider
the individual claims without relating them to the policies and assume that the amount
of these claims of the portfolio are still i.i.d. Generally speaking, the aggregate claims
of the portfolio is still a sum of i.i.d. random variables, the di erence between these
1
two models being the signi cation of the number of such variables: in the individual
model, the number of variables is the number of policies, while in the collective model
it is the number of claims, see book [49].
These two models can be extended to the multivariate setting. In an individual model, a
natural assumption could be that each policy holder has some policies, among which we
believe that there might be dependence. For such cases, we can extend the individual
model by replacing the one-dimensional aggregate claims of each policyholder with
a vector of elements representing his aggregate claims in, e.g.,home insurance, his
aggregate claims in automobile insurance, etc. We can consider a situation where
a claim event can cause claims on more than one policy. A typical example can be
considered extreme weather or earthquakes. In such situations, we could use a collective
model where one claim event produces a claim amount vector whose jth element is the
claim amount caused tot the jth policy by this claim event. We can also combine the
two settings and we shall present such a situation in Chapter 3.
Additionally, we will make a study case to show how an approximate technique such
as the Monte Carlo method can be used in ship design industry for risk analysis, see
e.g. Hullet and Avalon [40], J.R. van Dorp; M. R. Du ey [23], Vose [98].
1.2 Research objectives
The main objective of the research realized in the doctoral school program was
the study of some extension stochastic models concerning aggregate claims and their
implementation through several practical applications.
We extend such models from the univariate case to the multivariate case, models in
which we include some dependencies in order to help the science to make a more
realistic analysis of insurance scenarios.
In my opinion, the study of such models is a starting point in the theory of the aggregate
claims distributions and opens a variety of ways of research that other young future
PhD students could follow.
The operational objectives of my research were:
The extension to a multivariate settings of the bivariate collective model intro-
2
duced by Jin and Ren [43], to study aggregate claims when di erent types of
claims simultaneously a ect an insurance portfolio. The simpler multivatiate
collective models were studied by Sundt [94].
The description of several techniques used to evaluate the compound distribu-
tion such as: convolution, recursions, simulation, calculation with approximate
distribution and inversion methods, which also includes the FFT, and moreover
the study of how these methods interact with our new collective models in order
to obtain the distribution of the aggregate amount of claims occurred in an in
insurance portfolio within a given period of time. For details on these methods
see Klugmann [49], while for the FFT method see B uhlmann [12], Grubel and
Hermesmeier [31] and Embrechts et al. [26].
A comparison between these techniques mentioned above in order to nd the
optimal way to be closer to the reality;
We also focused on the so called alternative methods in order to simplify calcu-
lations and reduce the computing time, here we detailed a lot of aspects such as:
a variety of errors that these methods imply but also a strongly comparison be-
tween simulation and recursive methods. For details on such errors see Klugman
[49], Grubel and Hermesmeier [31], [32], and Sundt and Vernic(2009) [96].
A study case using the Monte Carlo method in order to demonstrate the utility
and the application of this method in the eld of engineering, more precisely in
the shipbuilding area. Fore more details regarding risk in shipbuilding see Dorp
and Du ey [23] and Kolic and Calic [50], and fore mode generic details see also
Storch and Lin, [92], and Winston [114] and Dlugokecki et al. [21].
1.3 PhD thesis structure
This PhD thesis is structured as it follows.
The rst chapter contains the motivation and the importance of the research topic,
the research targets and some comments regarding the bibliography.
3
The second chapter , Preliminary notions, as its name says, aims to familiarize the
reader with some general aspects such as: collective model – what it is and which are
the conditions to obtain such a model, together with a classi cation of such models
according to their dimension. Finally the chapter ends with a detailed review of the
techniques used to evaluate the compound distributions of stochastic models mentioned
above such as: convolutions method, Panjer recursive formula, inversion and simula-
tion.
In the third chapter , Recursions for compound multidimensional models, we extend
to a multivariate setting the two bivariate models introduced by [43]. The motivation
comes from the fact that some stochastic sources such as: res,
oods, trac acci-
dents or earthquakes may cause di erent types of dependent claims; for this reason
the recursions were made under the major assumption that the models include some
dependencies, together with other supplementary assumptions corresponding to each
model. The original results of this chapter are presented in [76] and [80].
If in the third chapter we present the method to obtain an exact recursive formula
for the probability function of the multivariate compound distribution, in Chapter 4
we deal with some alternative methods applied on the same models mentioned before.
We studied methods such as FFT and simulation in order to reduce rst of all the
computing time. There are also other issues that must be taken into consideration
when it comes to talking about these alternative methods, especially if we use the FFT
method, issues such as : errors. For that reason this chapter is also focused on the
related errors and their reduction. It is important to mention that all the numerical
calculations involved in this research were made in the Matlab software program, a
dedicated tool for such calculations, especially if we take into consideration the FFT.
The original results of this chapter can be found in [79] and [78].
The fth chapter may be considered the surprise of this thesis, due to the fact that
is it integrally a risk study case from the engineering eld, namely the industry of ship
design. Even though at a rst sight it may be extracurricular, it remains anchored in
the subject, because it is about how to perform risk analysis using the so called Monte
Carlo simulation. This study case was accepted to be presented at the 6th Interna-
tional Maritime Conference on Design for Safety, this autumn in Hamburg.
4
Finally, Chapter 6 , Conclusion, aims to bring together the importance of the subjects
studied in this thesis and the new perspective of research that may be developed in the
future in accordance with it.
5
Chapter 2
Preliminary notions
For insurance companies, studying the aggregate claims distribution, that is, the dis-
tribution of the aggregate amount of claims occurring in an insurance portfolio within
a given period, is crucial. Such distributions are usually modeled using individual or
collective models, which have already been introduced in Section 1.1. In the following
we concentrate on the collective models, which are the main subject of this thesis.
In our study we take into account the aspect that a probabilistic or a stochastic model
will always be an approximation to reality where one has to nd a balance between
realism and mathematichal simplicity, as stated by [84].
2.1 Collective models
2.1.1 Univariate case
The collective model in the univariate case can be written as
S=NX
j=0Uj; (2.1)
whereSrepresents the random variable (r.v.) aggregate claims, Nthe r.v. number
of claims occurring during the given time period, U0= 0 and (Uj)j1are independent
and identically distributed (i.i.d.) claim sizes positive r.v.s, which are also independent
ofNsee [35].
The distribution of Sis called compound, while the distribution of Nis called counting
6
distribution.
2.1.2 Bivariate case
The collective model in the bivariate case, can be written as
(S1;S2) = N1X
l=1Ul;N2X
j=1Vj!
; (2.2)
whereS1represents the aggregate losses of type I, S2the aggregate losses f type II, Ul
the cost of type I losses, (i.i.d.), Vjthe cost of type II losses, (i.i.d.), N1the number of
type I losses, N2the number of type II losses. Here N1andUlare independent while
N2;Vjare also independent.
Recently, [43] introduced a bivariate collective model to study aggregate claims in the
case when di erent types of claims simultaneously a ect an insurance portfolio (e.g.,
oods, storms or earthquakes).
Their model is
(S1;S2) = N1X
i=1Ui+N3X
k=1Lk;N2X
j=1Vj+N3X
k=1Qk!
; (2.3)
whereN1denotes the number of accidents that cause only type one claims, N2denotes
the number of accidents that cause only type two claims and N3denotes the number of
accidents that cause both types of claims. The claim number vector ( N1;N2;N3) has
probability function p(n1;n2;n3) =P(N1=n1;N2=n2;N3=n3); (Ui)i1and (Vj)j1
are mutually independent and independent of the claim numbers ( N1;N2;N3) and claim
sizes (Lk;Qk)k1, in the same manner claim sizes vectors ( Lk;Qk)k1are mutually
independent, identically distributed and independent of claim numbers ( N1;N2;N3)
and claim sizes UiandVj.
2.1.3 Multivariate case
The multivariate form for the collective model can be written as
(S1;:::;Sm) = N1X
i=1U1;:::;NmX
i=1Umi!
(2.4)
whereNkdenotes the number of claims of type kand (Uki)i1their corresponding
costs, which are i.i.d. and independent of the number of claims. [43] obtained re-
7
cursions for the bivariate form (2.3) under three di erent assumptions related to the
dependency structure of the claim numbers; the resulting models were named A, B and
C, corresponding to the similar ones from [35]. They also used FFT as an alternative
method.
Remark 2.1 We recall that a distribution belongs to the R1(a;b)class if its probability
function (p.f.) satis es the recursion
Pr (N=n) =
a+b
n
Pr (N=n 1);8n1;
for some constants a;b2R(for details on the Rkclasses see, e.g., [93] or [96]).
2.2 Techniques used to evaluate the compound dis-
tributions
The techniques used to evaluate the distribution of the aggregate amount of claims
occured in an insurance portfolio within a given period are the following: convolutions,
recursions, simulation, calculation with approximate distributions and inversion meth-
ods, which includes the Fast Fourier Transform (FFT)see [49].
From the beginning of this chapter it is very important to note that it can be shown
that in some cases, a portfolio heterogeneity can be treated mathematically as uniform,
so throughout this chapter, we will consider all risks of a portfolio identical.
Thus, we assume that we study an insurance portfolio within the time [0, t],tbeing
xed.
2.2.1 Convolutions method
We will start with the Convolutions Method, which is an obvious but a very laborious
technique. Let us denote with FUthe d.f. of a r.v. U, and also it is important to
mention the necessary formulas for a good description of the convolutions method,
namely:
8
The Convolution of order 0 is by de nition:
F0
U(x) =8
<
:0; x0
1; x> 0
The Convolution of order n>0;F0
U, is the d.f. of the sum U1+U2+:::+Unand
is given by
F1
U=FU (2.5)
Fn
U(x) =Zx
0F(n 1)
I (x y)dFU(y); n> 1 (2.6)
where it was taken into account that P(U > 0) = 1.
Remark 2.2 IfFUis absolutely continue and admits the density fU, then also Fn
Uis
absolutely continue,with the density fn
Ugiven by
f1
U=fU (2.7)
fn
U(x) =Zx
0f(n 1)
U (x y)fU(y)dy; n> 1 (2.8)
Remark 2.3 IfUis a discrete r.v., U2Nand if we denote
qx=P(U=x);x2Nthan the formula (2.6) becomes
q1
x=qx (2.9)
qn
x=P(U1+U2+:::;Un=x) =xX
y=1q(n 1)
x yqy; n> 1 (2.10)
whereq1
x=8
<
:1; x= 0
0;otherwise(2.11)
Proposition 2.4 Let beFSthe exact d.f. of S. Then it is given by
FS(x) =1X
n=0P(N=n)Fn
U(x);x> 0 (2.12)
Proof. The proof can be found in [49]
Remark 2.5 IfFUis absolutely continuous, it does not generally result that FSis also
absolutely continuous. In this case the distribution of Sis a mixed one:
FS(x)(2:12)=P(N= 0) +1X
n=1P(N=n)Zx
0fn
U(y)dy=
=p0+ (1 p(t))Zx
01X
n=1pn
1 p0fn
U(y)dy;
9
wherepn:=P(N=n).
We can observe that the formula of FSis a weighted sum between a degenerate discrete
distribution,with the value of 1, having the weight p0and an absolutely continuous
distribution with the weight (1 p0)and the density
1X
n=1pn
1 p0fn
U(y)
In conclusion, in order to have FS(x) absolutely continuous, it is necessary and sucient
thatp0= 0, in which case the density is
fS(x) =1X
n=1pnfn
U(x)
Remark 2.6 From the formula (2.12) we can observe that to calculate exactly the dis-
tribution of Sit is necessary to evaluate the convolutions of the r.v. U. Unfortunately,
these convolutions are generally dicult to calculate, even on the fastest computers.
Moreover, if FUis continuous or mixed, its convolutions require numerical procedures.
IfFUis discrete, convolution can be performed directly on the computer, but may re-
quire a lot of time, especially if F(x)must be calculated for a large number of values.
As a result, one may use speci c methods applicable in particular cases.
2.2.2 Panjer recursive formula
Compound distributions such as compound Poisson or Negative binomial are widely
used in modeling the distribution of cumulative claims that can occur in a xed period
of time.
The usual method of evaluating distribution functions requires many convolutions of
the claims distribution conditioned by the fact that the claim occurs. If the number
of alleged damage is very high, the calculation can become cumbersome even using
computers.
As long as the distribution of the claims number belongs to a speci c family (i.e. Pois-
son, binomial and negative binomial), an algorithm proposed by Panjer allows for the
recursive calculation of the distribution S, avoiding the calculation of convolutions.
In 1981, Panjer introduced a recursive formula to calculate the exact distribution of
10
the cumulative value of the losses in a particular case of the collective model. Since
then, a growing number of papers studied generalizations of this formula.
Recall the collective model (3.1) for a portfolio of insurance policies:
Srepresents the random variable (r.v.) aggregate claims, Nthe number of claims
andU0= 0;U1;U2;:::the costs of this claims, view as a selection under Uwith
P(U > 0) = 1, hence S=NP
j=0Uj.
We are in the case when Uis a discrete r.v. and suppose that its values can be arranged
so to t into the set f1;2;:::;rg;r2N. Let us denote pn=P(N=n);n2Nand
qj=P(U=j);j=1;r, with the condition qr>0. It is clear that the r.v. Sis also
dicrete, taking naturals values, so we have to evaluate P(S=x);x2N.
2.2.2.1 Panjer's Family
De nition 2.7 The family of probability distributions (pn)n2Nsatisfying:
p0>0 (2.13)
pn=
a+b
n
pn 1 n2N;a;b2R (2.14)
is called Panjer's Family.
Remark 2.8 pnit is the probability of exactly nclaims in a xed period of time.
Remark 2.9 The only distributions satisfying the relationship above are: Poisson,
Binomial and Negativ Binomial.
1. Poisson distribution ( )
pn=n
n!e ; n = 0;1;2;:::
p0=e ;pn
pn 1=
n=)a= 0;b= 0:
11
2. Binomial distribution ( t;)
pn=8
><
>:Cn
t(1 )t n; n = 0;1;2;:::;t
0; n =t+ 1;t+ 2;:::
p0= (1 )t;pn
pn 1=(t n+ 1)
n(1 )=)a=
( 1);b=(t+ 1)
(1 ):
3. Negative Binomial distribution ( NB( ;))
pn=C +n 1
nn(1 ) ; n = 0;1;2;:::
p0= (1 ) ;pn
pn 1=( +n 1)
n=)a=;b=( 1):
2.2.2.2 The recursive Panjer formula – discret case
Theorem 2.10 The recursive Panjer formula – discret case
LetShave a compound distribution for which the r.v. cost is taking non-negative inte-
gers values, with the distribution function f(x);x= 0;1;2;:::; and the probability nto
havepnclaims satisfy the recurrence (2.14).
The following formulas for the probability of the aggregate claims, s, hold:
P(S= 0) =g(0) =8
<
:P(N= 0); iff(0) = 0
GNf(0); iff(0)>0(2.15)
P(S=s) =g(s) =1
1 a(0)xX
y=1
a+by
x
f(y)g(x y); x = 1;2;:::(2.16)
whereGNis the p.g.f. of N.
Remark 2.11 The relationship (2.15), (2.16) for the three distribution: Poisson, Bi-
nomial and Negativ Binomial becomes:
1. Poisson ( ) witha= 0,b=
g(0) =e (1 f(0))
g(s) =1
sxX
y=1yf(y)g(x y)
12
2. Binomial with a=
1,b=
1 (t+ 1), (results that a<0)
g(x) =
(1 f(0)) 1xX
y=1
1 (t+ 1)y
x
f(y)g(x y); x = 1;2;:::
3. Negativ Binomial ( ;)ib=( 1), from where it results that 0< a < 1,
a+b>0and
g(x) =
1 f(0) 1xX
y=1
1 +( 1)y
x
f(y)g(x y); x = 1;2;::
2.2.2.3 The recursive Panjer Formula-continuous case
Theorem 2.12 The recursive Panjer Formula-continuous case
LetShave a compound distribution for which the r.v. cost is absolutely continuous,
with the density f taking values in the range of (0;1), and the probability pnof having
n claims satisfy the recurrence relation (2.14). Then the following relation holds:
g(x) =8
><
>:p0; ifx= 0
p1f(x) +xR
0
a+by
x
f(y)g(x y)dy ifx>0(2.17)
Remark 2.13 The integral equation (2.17)is called of Volterra type. It can be solved
with numerical methods, but in practice is easier to approximate the continuous distri-
bution of costs with the discrete one, using eventually an approximation which gives a
minimal or maximal limit for the exact distribution.
2.2.2.4 Extension of Panjer's recursion
In this section we shall brie
y present several literature references that extend Panjer's
recursions to other types of distributions.
In this sense, Panjer and Willmot [67] and Hesselager [34] discussed recursive evaluation
of compound distributions with counting distribution satisfying
p(n) =Pt
i=0c(i)ni
Pt
i=0d(i)nip(n 1);(n= 1;2;:::) (2.18)
Further on, in [2] it is discussed the recursion for a class of distributions pa;bthat satisfy
a relation in the form
pa;b=
u(a;b) +v(a;b)
n
pa+b;b(n 1):(n=l+ 1;l+ 2;:::) (2.19)
13
The special case where pa;bsatis es the relation
pa;b(n) =a
a+b
a+b
n
pa+b;b(n 1);(n= 1;2;:::) (2.20)
was treated by [3]. Moreover, [36] deduced recursions for a compound Lagrange distri-
butions and a compound shifted Lagrange distribution with kernel in the Panjer class.
Recursions in connection with Lagrange distributions have also been studied by [29],
[3], and [88].
Sundt [93] extended Panjer's class of counting distributions to
p(n) =nX
i=1
a(i) +b(i)
n
p(n 1): (2.21)
In [48] the generalised Charlier distribution has been studied p, which satis es a recur-
sion on the form
p(n) =
a+b
n
p(n 1) +
c+d
n+f
n 1
p(n 2) (n= 2;3;:::) (2.22)
The recursion
f(x) =1
1 (h(0)) mX
i=1 (i)hi(x) +xX
y=1f(x y)mX
i=1 (i)hi(y)!
(x= 1;2;:::)
was deduced by Eisele in [25] within the framework of phase-type distributions, in
the same manner [109] presented an algorithm for recursive evaluation of compound
distributions and counting distribution that satis es a recursion in the following form
p(n) =kX
i=1Pt
j=0c(i;j)nj
Pt
j=0d(j)njp(n i);(n=l+ 1;l+ 2;:::) (2.23)
Other recursive formulas concern the class of compound mixed Poisson distributions,
for more details see [85].
In the multivariate case, we mention the papers [94], [4],[101], [102],[103], [104], [ ?],[111],
[95], [38], [59], [60] and the overview [85]. More recent approaches are Wang et al. [110],
Minkova Balakrishnan [62], Anastasiadis Chukova [5],Lin, P., Neil Fenton [55] and
Liu Cheung [57].
14
2.2.3 Inversion methods using the Fourier transform
There are several inversion methods which are used to obtain in a numerical way com-
pound distributions or other elements derived from them, starting from the well known
expression as a transformed such as p.g.f., m.g.f. or c.f.. Compound distributions nat-
urally lead to such an approach, because their transforms are compound functions easy
to evaluate when both distribution of the claim frequency and their costs are known.
For example p.g.f. and c.f. of the aggregate claims are:
gS(u) =E(uS) =gN(gU(u));
'S(z) =E(eizX) =gN('U(z)).
The characteristic function is always unique and reciprocal, for a given c.f. there is
always a unique proper distribution. The scope of the following method is to obtain in
a numerical way a distribution starting from the c.f.
2.2.3.1 Fast Fourier Transform, FFT
De nition 2.14 For a certain function f:R !Rthe Fourier Transform is:
ef(z) =Z1
1f(x)eizxdx:
Remark 2.15 The starting function fmay be rebuilt from its Fourier transform through
f(x) =1
2Z1
1ef(z)e izxdz:
Remark 2.16 Iffis a density, then its Fourier transfom is its c.f.
De nition 2.17 Consider the vector (f0;f1;:::;fn 1)made ofnreal values, n2N.
Then we can de ne the discrete Fourier transform (DFT)attached to the vector as the
seriesef(k);k=0;n 1where
ef(k) =n 1X
j=0fje2i
njk; k=0;n 1 (2.24)
15
The transform ( f0;:::;fn 1)DFT !(ef0;:::;efn 1) is bijective, the inverse transform being
given by
fj=1
nn 1X
k=0ef(k)e 2i
njk; j=0;n 1 (2.25)
From (2.24) we can observe that in order to evaluate a value efkwe have to do the sum
ofnterms, so to evaluate all the nvaluesef0we have to sum and calculate n2terms
such asfje2i
njk. The calculation complexity will be of the order O(n2).
The algorithm FFT reduce the number of calculations at O(nlog 2n), and the reduction
becomes signi cant for big n. The algorithm is based to the property of the DFT of
lengthnto be written as a sum of two DFT having each of them the length n=2, the
rst being made of even points, and the second with the odd ones, thus:
ef(k) =efa
k+e2i
nkefb
k (2.26)
whereefa
kis the DFT coresponding to the vector ( f0;f2;:::;f 2[n=2]), andefb
kis DFT proper
of the vector ( f1;f3;:::;f 2[(n+1)=2] 1). On their turn, efa
kandefb
kcan be written as a sum
of two DFT with the length n=22, etc. In order that this numbers to be integers, the
algorithm FFT will begin with a vector of length n= 2r. The succesive writing of the
transforms in DFT having half of their length will conduct after rsteps at transforms
of length 1.
Knowing the transforms of length 1, we calculate successively the ones of length
2;22;23;:::;2rbased on the formula (2.26).
2.2.4 Simulation
As discussed by [49], the analytic methods presented before have two common charac-
teristics. First of all, they are exact until a speci c level of approximation, although
the errors can be reduced to almost 0. Secondly, both recursive and inverse methods
suppose that the aggregate claims can be written in the form S=U1+:::+UN, where
N;U 1;:::;UNare independent and U1;:::;UNare identically distributed. If the rst
characteristic doesn't concern us too much (the errors can be taken theoretically small
16
enough), the second one is restrictive and may separate the model from the reality. We
will present in the following a method which does not require restrictions like indepen-
dence or identically distribution, namely, we present simulation.
There are two typical situations in which the independence assumptions and identically
distributions are not valid. The rst one requires consideration of the time factor (in
particular, the change in value of money over time), and the second is due to changes
in deductibles, and it may also involve the time factor.
The beauty of simulation is that once the model is created, distribution of the aggre-
gated value will result without other inconvenience. However, to perform an ecient
simulation, it is necessary to have more creativity. In the following we present an
approach quite "raw", which works, but might require more running time on the com-
puter. The whole process can be summarized in three steps, the aim being to determine
the distribution function of the random variable S. Generally speaking, assuming that
we have to determine the distribution function for a random variable S, simulation
technique consists of the following steps:
1. It builds a model for Sby considering the r.v.s Y;Z;V;:::; whose distribution and
dependencies are known.
2. Fori= 1;:::;n pseudo-random values yi;zi;vi;:::;are generated yielding the val-
uessias resulting from the model in Step 1.
3. The distribution function of Sis approximated by the empirical distribution
function obtained from the values s1;:::;sn.
In particular, to obtain the distribution function for a compound random variable
S=NP
i=1Uithrough simulation we do the following:
Simulate a value nofN, of which distribution it is assumed to be known, for
examplePoisson ();
Simulatenvalues from the random variable claim U, whose distribution is known.
Sum ofnvalues ofUand obtain a value sofS;
Repeat Steps 1-2 mtimes, in order to obtain the values s1;:::;sm;
17
Evaluate the empirical distribution function
~F(a) =jfsi=si<a;i =1:mgj
m:
We have to answer two questions: what does it mean to generate pseudo-random values
and how large nshould be ?
Regarding the rst question, let Ube a r.v. with the distribution function F. This is
a real r.v. caused by a phenomenon of interest. For example, it can occur as a result
of an experiment like "collect random health care costs arising from car accidents".
Assume that the distribution of Uis known, for example Pareto with d.f. :
F(x) = 1 1000
1000 +x3
;x> 0
If we have a sample ( x1;:::;xn) extracted from the above Pareto r.v., it would be in-
distinguishable from a sample extracted from some other process, but with the same
Pareto distributions? However, to obtain a Pareto distributed sample exactly is dif-
cult, so we have to use another method; we can make progress by making some
concessions. Thus, we accept the replacement of a sample of Uwith a series of num-
bers (x
1;:::;x
n), which is not a random sample, but a sequence of numbers which may
not even be independent, but have been generated by a known process related to r.v.
U. Such a sequence of numbers is called pseudo-random since its origin is generally
not known, but can not be distinguished from a random sample of U. Such a sequence
will be sucient for our purpose.
The eld generating pseudo-random number sequences is well developed. One thing
that simpli es this generation is the ability to generate uniformly distributed numbers
in (0; 1). This is because we have the following property: "if VUniform (0; 1), than
the r.v.Y=F 1(V) will have the d.f. F". Therefore, it is enough to generate the
pseudo-randomly numbers ( u
1;:::;u
n) and then calculate x
i=F 1(u
i).
Let us consider now the second question : what value should we use for n? It is known
that an estimate is closer to the actual value when the sample size grows. With a little
e ort, we could determine the value of nthat will get us as close as possible to the real
value (i.e. with a speci ed probability).
An example of simulation for the collective model
18
Next we present the simulation algorithm for a collective model of claims, in the uni-
variate case:
Algorithm:
Simulates a value nofN, choosing from (Poisson, Binomial, Binomial Negative);
Simulates the cost Uichoosing between LogNormal and Exponential distribu-
tions; simulates in total nvalues forU;
CalculateSby summing the values u1+:::+un;
To obtain the empirical c.d.f. of Sevaluated in a, check if the value obtained for
Sis smaller than the value a;
Repeat steps 1-5 mtimes.
Use the formula ~F(a) =jfsi=si<a;i =1:mgj
m.
2.2.4.1 Monte Carlo method
The simulation algorithms presented above are based on the Monte Carlo method, see
e.g. Binder [8], Dunn and Shultis [24], Glasserman [28], Kroese et al. [51], Robert
[81], Rubinstein and Kroes [83]. Even if we use an analytic method or a simulation,
rst of all wee need to build a mathematical model of a system. Once we create it the
next step is to better understand and control it as eciently as possible. Such com-
plex stochastic systems should come from eld such as : telecommunications networks,
insurance, weather forecasting, and why not engineering. It is well known that an ana-
litical solution gives us the exact solution and requires far less work than simulation,
but if the scenario becomes more complex, simulation is the most indicated choice.
For that reason in Chapter 4 we present a study case from the naval engineering eld,
more exactly a risk analysis, which was made in a dedicated program software, Primav-
era Risk Analysis, program which has inside the Monte Carlo simulation method. The
rst question may be: how does it function? But once familiarized with the concept
of the classical simulation the answer may be very simple. In other words the Monte
Carlo method is the computerized mathematical technique of the classical simulation
19
mentioned above, dedicated especially to performing risk analysis by building models
of possible results by substituting a range of values with a probability distribution,
for example: Normal, LogNormal, Triangular, Uniform, etc., for any factor that has
inherent uncertainty. Next it calculates results over and over, each time using a dif-
ferent set of random values from the probability functions. Depending on the number
of uncertainties and the ranges speci ed for them, a Monte Carlo simulation could
involve thousands or tens of thousands of recalculations before it is complete. The
Monte Carlo simulation produces distributions of possible outcome values.
20
Chapter 3
Recursions for compound
multidimensional models
It is well known that for the evaluation of premiums and ruin probabilities, actuaries
study aggregate claims distributions corresponding to insurance portfolios. These dis-
tributions result from individual or collective models. In the following, we deal with
the collective model in a multivariate setting motivated by the fact that some stochas-
tic sources (like, e.g., res,
oods, trac accidents, earthquakes) may cause di erent
types of dependent claims. We start by extending to a multivariate setting the two
bivariate models A and B introduced by [43] to model insurance aggregate claims that
simultaneously a ect an insurance portfolio. Our goal was to obtain an exact recursive
formula for the probability function of the multivariate compound distribution corre-
sponding to both models A and B mentioned before.
These recursions were made under the major assumptions that the conditional multi-
variate counting distribution is multivariate Poisson for the model B and,conditioned
by the total number of claims, is multinomial for the model A. Moreover, regarding
the model A, we develop an alternative proof for the exact recursion for the evaluation
of the same distribution. We should mention that other supplementary assumptions
were used, but these ones will be presented in the following.
21
3.1 Introduction
Therefore, we consider the following multivariate aggregate claims model
(S1;:::;Sm) = N1X
i=0U1l+N0X
k=0L1k;:::;NmX
l=0Uml+N0X
k=0Lmk!
; (3.1)
wherem2 is the number of di erent types of claims a ecting the portfolio, Skdenotes
the aggregate claims of type k,Nkthe number of claims of only type k,N0the number of
common claims (e.g., accidents that causes all mtypes of di erent claims). Each set of
claim sizes ( Ujl)l1 are positive, independent and identically distributed (i.i.d.) as the
generic random variable (r.v.) Uj;1jm, independent of the claim numbers and of
the other claim sizes, including ( L1k;:::;Lmk). The random vectors ( L1k;:::;Lmk)k1are
non-negative i.i.d. as the generic ( L1;:::;Lm), and independent of the claim numbers.
Clearly,Uj0=Lj0= 0;8j. In the following, by a bold faced letter we denote a vector,
i.e.,X= (X1;:::;Xm) orx= (x1;:::;xm). We shall work with discrete r.v.s and if
the claim sizes distributions are continuous, they should be discretized using, e.g., the
rounding method, see [49]. If fis a probability function (p.f.), we denote by fnits
n-fold convolution corresponding to the distribution of the sum of ni.i.d. r.v.s having
p.f.f, note thatf1=fand, by convention, f0(x) =8
><
>:1x= 0
0x6= 0. LetfSdenote the
p.f. of S,fjthe p.f. ofUj,j=1;m,f0the p.f. of L, andpthe p.f. of N= (N0;:::;Nm).
Then from [49], we easily obtain
fS(x) =X
n00;:::;nm0p(n0;:::;nm)xX
k=0mY
j=1fnj
j(kj)fn0
0(x k);x0; (3.2)
where 0= (0;:::;0) , while the inequality x0and the di erence x-kare componen-
twise. Note that due to the convolutions in (3.2), fScould be dicult to evaluate
and time-consuming; therefore, alternative methods have been developed, from which
recursions play an important role in actuarial mathematics. We recall the fact that
the distribution of Nis also called counting distribution, while the distribution of Sis
called compound. Let us denote the p.f. of a discrete random vector XbyfXand its
probability generating function (p.g.f.) by gX; we recall that
gX(t)def=E"mY
j=1tXj
j#
:
22
Moreover, as a property of the p.g.f., it holds that
gX(t) =X
x0fX(X)mY
i=1txi
i (3.3)
and clearly, gx(0) =fx(0).
Proposition 3.1 Under the assumptions of model 3.1, the p.g.f. of Sis given by
gS(t) =gN(gL(t);gU1(t1);:::;gUm(tm)): (3.4)
Proof. We have that
gS(t) =E"mY
j=1tSj
j#
=E"
E"mY
j=1tPNj
l=0Ujl+PN0
k=0Ljk
j N0;:::;Nm##
=E2
4mY
j=1Eh
tUj
jiNjE"mY
j=1tLj
j#N03
5
=E"
gN0
L(t)mY
j=1gNj
Uj(tj)#
;
which immediately yields formula (3.4) and completes the proof.
3.2 Recursive evaluation: rst case, corresponding
to model B
This section is based on paper [76].
In this next section, we shall present a recursion to evaluate for model (1) when the
number of claims Nfollows a multivariate Poisson distribution. In the bivariate setting,
whenm= 2, a recursion for this case has been recently presented in [43]; our recursion
extends this one to a general m. In the simpler univariate case ( m= 1), the history
of recursions involving Poisson counting distributions starts in insurance with [68]
and [112], and continues with more complex recursions for compound mixed Poisson
distributions discussed in [113], [35] among others, or, in the multivariate case, in [84],
etc. In particular, as we mentioned from the beginnig if Nfollows the multivaritate
Poisson distribution MPom+1; 0;:::;mwith>0;j>0;, then, from [45] we have
23
that
gN(t) = exp(
mY
j=0tj 1!
+mX
j=0j(tj 1))
: (3.5)
Next we obtain a recursion for the p.f. for model (3.1) when Nfollows a multivariate
Poisson distribution.
Proposition 3.2 Under the assumptions of model (3.1), if Nfollows the multivariate
Poisson distribution MPom+1(; 0;:::;m), then the p.f. of Ssatis es the recursion
fS(X) =k
xkxkX
zk=0zkfk(zk)fS(x1;:::;xk 1;xk zk;xk+1;:::;xm+
+0
xkxX
z=0zkf0(z)fs(x z) +
xkxX
v=0 vX
u=0vkfi(ui)f0(v u)!
;xk1;1km;
with starting value
fS(0) = exp(
f0(0)mY
j=1fj(0) 1!
+mX
j=1j(fj(0) 1) +0(f0(0) 1))
:(3.6)
Proof. We start by inserting the p.g.f. (3.3) of Ninto (3.4) and obtain
gs(t) = exp(
gL(t)mY
j=1gUj(tj) 1!
+0(gL(t) 1))
:
Taking here t= 0 and using the properties of the p.g.f., we easily get the starting value
(3.6). To obtain the recursion, we recall from (3.3) that
gS(t) =X
x0fS(x)mY
l=1txl
l;gL(t) =X
x0f0(x)mY
l=1txl
l;gUk(tk) =1X
x=0fk(x)tx
k;
which yields
@gs(t)
@tk=X
x0fS(x)xktxk 1
kY
l=k
l6=ktxl
l;@gL(t)
@tk=X
x0f0(x)xktxk 1
kmY
l=1
l6=ktxl
l;g0
Uk(tk) =1X
x=0fk(x)txk 1
k:
(3.7)
We now consider the partial derivative of gswith respect to tk, 1, i.e.,
@gs(t)
@tk=gs(t)2
40
@@gL(t)
@tkm
l=1gUj(tj) +gL(t)dgUk(tk)
dtkm
j=1
j6=kgUj(tj)1
A+
+kdgUk(tk)
dtk+0@gL(t)
@tk
:
24
Inserting here the formulas (3.7) leads to
X
x0fS(x)xktxk 1
kY
l=k
l6=ktxl
l=T0(kT1+0T2+(T3+T4)); (3.8)
where
T0=gs(t) =X
y0fs(y)mY
l=1tyl
l;T1=g0
Uk(tk) =1X
zk=0fk(zk)zktzk 1
k;
T2=@gL(t)
@tk=X
z0f0(z)zktzk 1
kmY
l=1
l6=ktzl
l;T3=@gL(t)
@tkmY
j=1gUj(tj) =T2mY
j=11X
uj=0fj(uj)tuj
j;
T4=gL(t)g0
Uk(tk)mY
j=1
j6=kgUj(tj) =X
z0f0(z)mY
l=1tzl
l 1X
uk=0fk(uk)uktuk l
k!mY
j=1
j6=k1X
uj=0fj(uj)tuj
j:
We separately evaluate each term of(3.8). We start with
T0T1=X
y01X
zk=0fS(y)fk(zk)zktzk+yk 1
kmY
l=1
l6=ktyl
l;
where we change the variable xk=yk+zk)yk=xk zk, we interchange the sums
and obtain
T0T1=1X
yj=0
8j6=k1X
zk=01X
xk=zkfk(zk)zkfS(y1;:::;yk 1;xk zk;yk+1;:::;ym)txk 1
kmY
l=1
l6=ktyl
l
=1X
yj=0
8j6=k1X
xk=01X
zk=0zkfk(zk)fS(y1;:::;yk 1;xk zk;yk+1;:::;ym)txk 1
kmY
l=1
l6=ktyl
l
Similarly, changing variabiles xj=yj+zj;, and interchanging the sums, yields
T0T2=X
y0X
z0fs(y)f0(z)zktzk+yk 1
kY
l=1
l6=ktyl+zl
l=X
z0X
xzfS(x z)f0(z)zktxk 1
kmY
l=1
l6=ktxl
l
=X
x0xX
z=0fs(x z)f0(z)zktxk 1
kmY
l=1
l6=ktxl
l:
25
Also,
T3+T4=X
z0f0(z)0
B@zktxk 1
k 1X
uk=0fk(uk)tuk
k!mY
l=1
l6=k1X
ul=0fl(ul)tul+zl
l
+ 1X
uk=0ukfk(uk)tzl+ul 1
k!m
l=1
l6=k1X
ul=0fl(ul)tul+zl
l1
CA
=X
z0f0(z)0
B@mY
l=1
l6=k1X
ul=0fl(ul)tul+zl
l1
CA1X
uk=0fk(uk)tzl+ul 1
k (uk zk);
from where
T0(T3+T4) =X
y0fs(y)X
z0f0(z)0
B@mY
l=1
l6=k1X
ul=0fl(ul)tul+zl+yl
l1
CA1X
uk=0fk(uk)tzkuk+yk 1
k (uk+zk);
where we change variabiles vj=uj+zjandxj=vj+yj=uj+zj+yj;, which gives
T0(T3+T4) =X
u0X
vuX
xufs(x v)f0(v u)vkfk(uk)txk 1
kmY
l=1
l6=kfl(ul)txl
l
=X
x0xX
v=0vX
u=0fS(x v)f0(v u)vkfk(uk)txk 1
kmY
l=1
l6=kfl(ul)txl
l:
Inserting all these formulas into (3.8) and identifying the coecients of txk 1
kQm
l=1;l6=ktxl
l
yields
xkfs(x) =kxkX
zk=0zkfk(zk)fs(x1;:::;xk 1;xk zk;xk+1;:::;xm)
+0xX
z=0zkf0(z)fs(x z) +xX
v=0vX
u=0vk(m
l=1fl(ul))f0(v u)fs(x v);
from where, if xk1, we easily obtain formula (3.6).
Particular case. To see how the recursion presented in Proposition 3.2 works,
we illustrate it on the particular case m= 3. The evaluation of fsstarts with for-
mula (3.6), which gives fs(0;0;0); using then (3.6), we calculate fs(1;0;0) by tak-
ingk= 1,fs(0;1;0) by taking k= 2 andfs(0;0;1) by taking k= 3. For exam-
ple, we shall use fs(1;0;0) = [1f1(1) +0f0(1;0;0) +(f1(0)f2(0)f3(0)f0(1;0;0) +
26
f1(1)f2(0)f3(0)f0(0;0;0))]fs(0;0;0), etc.
Next step consits of calculating fs(1;1;0);fs(1;0;1);fs(0;1;1) followed by the evalua-
tion offs(1;1;1) and so on.
Remark 3.3 If, in particular, we take m= 2, formula (3.6) reduces to the particular
formulas obtained by Jin and Ren in [43] for their model B with the trivariate Poisson
as counting distribution. Moreover, using a similar reasoning as in the above proof, the
recursion from Proposition 3.2 can be extended after some tedious calculation to a more
general (than the multivariate Poisson) counting distribution, built like the one in the
general model B of [43] (which is in fact our model (3.1) for m= 2); more precisely,
the counting distribution results by taking Nj=Zj+Z;0with independent Zjs, also
independent of Z, all of them ( Zincluded) being distributed in the so-called Panjer
class, which consists of the Poisson, binomial and negative binomial distributions (for
more details on the Panjer class of distributions see [85]). In conclusion, Theorem 2.2
in [43] can be extended from the bivariate case considered in that work to .
3.3 Recursive evaluation: second case correspond-
ing to model A
This paper is based on [80].
In this case, we start by extending the recursion of model A from [43] in a multivariate
setting under the assumption that the distribution of the total number of claims N
belongs to the R1(a;b) class, while, conditionally on it, the distribution of Nis multi-
nomial.
Therefore, we derive the recursive formula corresponding to the multivariate model A,
and provide a numerical example in the trivariate case m= 3 making some remarks in
what regarding the computations order.
According to the model A in Jin and Ren (2014), we consider the following supplemen-
tary assumptions on model (3.1):
A1The rst one is related to the total number of claims N=N0+N1+:::+Nm;
27
whose probability function (p.f.) is assumed to satisfy the Panjer-type recursion
Pr (N=n) =
a+b
n
Pr (N=n 1);8n1;
for some constants a;b2R(for details on Panjer's class, see [66] or [96]);
A2The second one concerns the conditional distribution of NgivenN=n;which is
assumed to be multinomial Mnom (n;p1;:::;pm) with the parameters n2Nand
p1;:::;pm2(0;1) such that p0:= 1 Pm
i=1pi2(0;1):We recall that (see, e.g.,
[45]), withn=Pm
i=0ni;
Pr(N0=n0;N1=n1;:::;Nm=nmjN=n) =n!
mQ
i=0ni!mY
i=0pni
i:
Moreover, denoting by Ethe expected value operator, based on the p.g.f. formula of
this multinomial distribution (see, e.g.,[45]), the p.g.f. of Nbecomes
gN(n) =E"
E"mY
j=0nNj
j N##
=E2
4 mX
j=0pjnj!N3
5=gN mX
j=0pjnj!
:
From Proposition 2.1 in [76] we have that for the general model (3.1), the p.g.f. of S
is given by
gS(t) =gN(gL(t);gU1(t1);:::;gUm(tm));
and inserting the above formula of gNyields that for model A,
gS(t) =gN mX
j=1pjgUj(tj) +p0gL(t)!
: (3.9)
3.3.1 Method 1
Based on this formula and on the properties of the p.g.f., we shall now obtain a recursion
for the p.f. of S. In the following, we redenote fL:f0
Proposition 3.4 Under the assumptions of the above model A, the following starting
value and recursive formula hold:
fS(0) =gS(0) =gN mX
j=1pjfj(0) +p0fL(0)!
;
28
fS(x) =amX
j=1
j6=kpjxjX
zj=0fj(zj)fS(x1;:::;xj 1;xj zj;xj+1;:::;xm)
+pkxkX
zk=0
a+bzk
xk
fk(zk)fS(x1;:::;xk 1;xk zk;xk+1;:::;xm)
+p0xX
z=0
a+bzk
xk
fL(z)fS(x z); xk1: (3.10)
Proof. From the assumption N2R1(a;b) it is easy to obtain that (1 at)g0
N(t) =
(a+b)gN(t);or, equivalently,
g0
N(t) =a+b
1 atgN(t): (3.11)
By di erentiating (3.14) with respect to tk;1km;and inserting (3.11), we obtain
@gS(t)
@tk=@
@tkgN mX
j=1pjgUj(tj) +p0gL(t)!
=(a+b)gS(t)
1 a
mP
j=1pjgUj(tj) +p0gL(t)!@
@tk mX
j=1pjgUj(tj) +p0gL(t)!
=(a+b)gS(t)
1 a
mP
j=1pjgUj(tj) +p0gL(t)!
pkg0
Uk(tk) +p0@gL(t)
@tk
;
hence
@gS(t)
@tk"
1 a mX
j=1pjgUj(tj) +p0gL(t)!#
= (a+b)gS(t)
pkg0
Uk(tk) +p0@gL(t)
@tk
;
and rearranging,
@gS(t)
@tk=aT1+aT2+ (a+b) (T3+T4); (3.12)
where we used the following notation
T1=@gS(t)
@tkmX
j=1pjgUj(tj); T2=@gS(t)
@tkp0gL(t);
T3=gS(t)pkg0
Uk(tk); T4=gS(t)p0@gL(t)
@tk:
29
We recall the following formulas
gS(t) =X
y0fS(y)mY
j=1tyj
j;@gS(t)
@tk=X
y0fS(y)yktyk 1
kmY
j=1
j6=ktyj
j;
gL(t) =X
z0fL(z)mY
j=1tzj
j;@gL(t)
@tk=X
z0fL(z)zktzk 1
kmY
j=1
j6=ktzj
j;
gUj(tj) =1X
zj=0fj(zj)tzj
j; g0
Uk(tk) =1X
zk=0fk(zk)zktzk 1
k:
Note that to simplify the writing, in the above derivatives we did not restricted, e.g.,
zk1 as it should be, since the term zktzk 1
kvanishes for zk= 0. Inserting now these
intoT1, we obtain
T1=2
4mX
j=1pj0
@1X
zj=0fj(zj)tzj
j1
A3
50
B@X
y0fS(y)yktyk 1
kmY
l=1
l6=ktyl
l1
CA=X
y0fS(y)
2
664mX
j=1
j6=kpj1X
zj=0fj(zj)yktyk 1
k0
B@mY
l=1
l6=k;l6=jtyl
l1
CAtyj+zj
j +pk1X
zk=0fk(zk)yktyk+zk 1
k0
B@mY
l=1
l6=ktyl
l1
CA3
775:
We change variable xj=yj+zj;from where yj=xj zj;hencexjzjand
T1=mX
j=1
j6=kpj1X
xj=01X
yl=0
8l6=jxjX
zj=0yktyk 1
k0
B@mY
l=1
l6=k;l6=jtyl
l1
CAtxj
jfj(zj)fS(y1;:::;yj 1;xj zj;yj+1;:::;ym)
+pk1X
xk=01X
yj=0
8j6=kxkX
zk=0(xk zk)txk 1
k0
B@mY
l=1
l6=ktyl
l1
CAfk(zk)fS(y1;:::;yk 1;xk zk;yk+1;:::;ym)
=X
x0mX
j=1
j6=kpjxjX
zj=0xktxk 1
k0
B@mY
l=1
l6=ktxl
l1
CAfj(zj)fS(x1;:::;xj 1;xj zj;xj+1;:::;xm)
+pkX
x0xkX
zk=0(xk zk)txk 1
k0
B@mY
l=1
l6=ktxl
l1
CAfk(zk)fS(x1;:::;xk 1;xk zk;xk+1;:::;xm):
Similarly,
T2=p0X
y0X
z0fS(y)fL(z)yktyk+zk 1
k0
B@mY
l=1
l6=ktyl+zl
l1
CA;
30
where we change variables xl=yl+zlandxk=yk+zk;yielding
T2=p0X
x0xX
z=0fL(z)fS(x z)(xk zk)txk 1
k0
B@mY
l=1
l6=ktxl
l1
CA:
Also,
T3=pkX
y01X
zk=0fS(y)0
B@mY
l=1
l6=ktyl
l1
CAzkfk(zk)tzk+yk 1
k;
and denoting xk=yk+zk;xl=yl;8l6=k, we obtain
T3=pkX
x0xkX
zk=0zktxk 1
k0
B@mY
l=1
l6=ktxl
l1
CAfk(zk)fS(x1;:::;xk 1;xk zk;xk+1;:::;xm);
while
T4=p0X
y0X
z0fS(y)fL(z)zktzk 1+yk
k0
B@mY
l=1
l6=ktzl+yl
l1
CA;
withxl=yl+zl;8l;gives
T4=p0X
x0xX
z=0zkfL(z)fS(x z)txk 1
k0
B@mY
l=1
l6=ktxl
l1
CA:
Inserting now the obtained formulas of T1;T2;T3;T4into (3.12) and identifying the
coecients of txk 1
kmQ
l=1
l6=ktxl
lyields
xkfS(x) =a2
664mX
j=1
j6=kpjxjX
zj=0xkfj(zj)fS(x1;:::;xj 1;xj zj;xj+1;:::;xm)
+pkxkX
zk=0(xk zk)fk(zk)fS(x1;:::;xk 1;xk zk;xk+1;:::;xm)
+p0xX
z=0(xk zk)fL(z)fS(x z)#
+ (a+b)"
p0xX
z=0zkfL(z)fS(x z)
+pkxkX
zk=0zkfk(zk)fS(x1;:::;xk 1;xk zk;xk+1;:::;xm)#
;
31
which we divide by xkand obtain
fS(x) =amX
j=1
j6=kpjxjX
zj=0fj(zj)fS(x1;:::;xj 1;xj zj;xj+1;:::;xm)
+pkxkX
zk=0a(xk zk) + (a+b)zk
xkfk(zk)fS(x1;:::;xk 1;xk zk;xk+1;:::;xm)
+p0xX
z=0a(xk zk) + (a+b)zk
xkfL(z)fS(x z);
from where the recursion (4.2) is immediate.
Remark 3.5 Takingm= 2 in (4.2) yields the formulas in Theorem 2.1 from [43].
Form= 3 andx1>0;(4.2) becomes
fS(x1;x2;x3)
=K"
ap2x2X
z=1f2(z)fS(x1;x2 z;x 3) +ap3x3X
z=1f3(z)fS(x1;x2;x3 z)
+p1x1X
z=0
a+bz
x1
f1(z)fS(x1 z;x 2;x3)
+p0x1X
z1=0x2X
z2=0x3X
z3=0
z1+z2+z3>0
a+bz1
x1
fL(z1;z2;z3)fS(x1 z1;x2 z2;x3 z3)3
775;(3.13)
whereK=
1 a P3
i=1pifi(0) +p0fL(0;0;0) 1:We obtain similar formulas for
x2>0andx3>0:
3.3.1.1 Numerical example
We shall now numerically illustrate the recursion in the trivariate case m= 3, by
assuming that the distribution of the total number of claims is Poisson, i.e., N
Po(); > 0. We recall that this distribution belongs to the R1(a;b) class with
a= 0;b=(see, e.g., Sundt and Vernic, [96]). Then formula (3.13) reduces for x1>0
to
fS(x1;x2;x3) =b
x1"
p1x1X
z=1zf1(z)fS(x1 z;x 2;x3)
+p0x1X
z1=1x2X
zz=0x3X
z3=0z1fL(z1;z2;z3)fS(x1 z1;x2 z2;x3 z3)#
;
32
and we can write similar formulas for x2>0 andx3>0:To obtain fSup to
fS(r1;r2;r3);we implemented these formulas in Matlab in the following way: rstly,
we evaluated the starting value
fS(0;0;0) = exp(
3X
j=1pjfj(0) +p0fL(0;0;0) 1!)
;
secondly, we evaluated fS(0;0;x3) andfS(0;x2;0) for 1xjrj;j= 1;2; then we
evaluatedfS(0;x2;x3) for 1xjrj;j= 1;2; and nally, all values fS(x1;x2;x3) for
1x1r1and 0xjrj;j= 1;2.
Numerically, we considered the following claim sizes distributions
f10
@0 1 2 3
0:3 0:2 0:3 0:21
A; f20
@0 1 2 3
0:4 0:1 0:3 0:21
A; f30
@0 1 2 3
0:2 0:3 0:4 0:11
A;
(fL(0;i;j))i;j=0;1=2
40:15 0:1
0:05 0:23
5;(fL(1;i;j))i;j=0;1=2
40:2 0:12
0:1 0:083
5;
while the multinomial parameters values are n= 5;p0= 0:25;p1= 0:25;p2= 0:3;p3=
0:2, and the Poisson parameter = 5. We obtained the following results, computed
top-down and left-right:
fS(0;0;0) = 0:0263
fS(0;0;1) = 0:0112fS(0;0;2) = 0:0129fS(0;0;3) = 0:0074
fS(0;1;0) = 0:0056fS(0;2;0) = 0:0124fS(0;3;0) = 0:0105
fS(0;1;1) = 0:0090fS(0;1;2) = 0:0055fS(0;1;3) = 0:0048
fS(0;2;1) = 0:0067fS(0;2;2) = 0:0075fS(0;2;3) = 0:0046
fS(0;3;1) = 0:0076fS(0;3;2) = 0:0066fS(0;3;3) = 0:0046
fS(1;0;0) = 0:0132fS(1;0;1) = 0:0095fS(1;0;2) = 0:0081fS(1;0;3) = 0:0057
fS(1;1;0) = 0:0061fS(1;1;1) = 0:0093fS(1;1;2) = 0:0068fS(1;1;3) = 0:0055
fS(1;2;0) = 0:0069fS(1;2;1) = 0:0069fS(1;2;2) = 0:0063fS(1;2;3) = 0:0046
… … … …
fS(3;3;0) = 0:0056fS(3;3;1) = 0:0055fS(3;3;2) = 0:0048fS(3;3;3) = 0:0036
etc.
33
3.3.2 Method 2
This section is based on [77].
As we mentioned in the beginning of this chapter we develop a second proof of the exact
recursive formula for the probability function of the multivariate compound distribution
developed in Section 3.3.1. The proof of the recursion presented in Proposition 3.1 is
based on properties of the probability generating function (p.g.f.).We shall now present
a shorter proof based on similar recursions already proved in [94]. In the same time,
we shall also obtain the new alternative recursive formula (3.16).
From [77], we have that for the general model (3.1) under the assumptions (A1-A2),
the p.g.f. of Sis given by
gS(t) =gN mX
j=1pjgUj(tj) +p0gL(t)!
: (3.14)
Using the properties of the p.g.f., the following recursive formula (4.2) has been proved
in [77].
Proposition 3.6 Under the assumptions (A1-A2) of model (3.1), the following start-
ing value and recursive formulas hold:
fs(0) =gS(0) =gN mX
j=1pjfj(0) +p0fL(0)!
;
fS(x) =K8
>><
>>:amX
j=1
j6=kpjxjX
yj=1fj(yj)fS(x1;:::;xj 1;xj yj;xj+1;:::;xm)
+pkxkX
yk=1
a+byk
xk
fk(yk)fS(x1;:::;xk 1;xk yk;xk+1;:::;xm)
+p0X
0<yx
a+byk
xk
fL(y)fS(x y))
;xk1;xj0;8j6=k;(3.15)
fS(x) =K8
<
:mX
j=1pjxjX
yj=1
a+byj
x+
fj(yj)fS(x1;:::;xj 1;xj yj;xj+1;:::;xm)
+p0X
0<yx
a+by+
x+
fL(y)fS(x y))
;x>0; (3.16)
34
whereK="
1 a
mP
j=1pjfj(0) +p0fL(0)!# 1
andx+=mP
i=1xi:
Proof. Under the assumptions (1-2), model (3.1) can be represented as
(S1;:::;Sm) =NX
k=0(C1k;:::;Cmk);
where, as before, N=N0+N1+:::+Nm;Cj0= 0;8j=1;m;while the random vectors
(C1k;:::;Cmk)k1are i.i.d. as the generic Chaving the p.f.
fC(y) =mX
j=1I(yi= 0;8i=1;m;i6=j)pjfj(yj) +p0fL(y): (3.17)
HereIdenotes the indicator function de ned by I(A) = 1 ifAis true, and 0 otherwise.
From [94] (see also [96], formulas (15.4) and (15.5), respectively), it holds that
fS(x) =1
1 afC(0)X
0<yx
a+byk
xk
fC(y)fS(x y); xk1; (3.18)
fS(x) =1
1 afC(0)X
0<yx
a+by+
x+
fC(y)fS(x y);x>0: (3.19)
By inserting (3.17) into (3.18) we obtain for xk1;
fS(x) =1
1 a
mP
j=1pjfj(0) +p0fL(0)!
(mX
j=1pjX
0<yx
a+byk
xk
I(yi= 0;8i=1;m;i6=j)fj(yj)fS(x y)
+p0X
0<yx
a+byk
xk
fL(y)fS(x y))
:
In the rst sum, we separately consider the cases j=kandj6=k.
– Whenj=k, combining I
yi= 0;8i=1;m;i6=k
withy>0yieldsyk1;hence
the middle term of (4.2);
– Whenj6=k, fromI
yi= 0;8i=1;m;i6=j
we clearly have yk= 0, therefore
X
0<yx
a+byk
xk
I(yi= 0;8i=1;m;i6=j)fj(yj)fS(x y)
becomesxjX
yj=1afj(yj)fS(x1;:::;xj 1;xj yj;xj+1;:::;xm);
35
from where formula (4.2) is immediate.
Similarly, by inserting (3.17) into (3.19) we obtain formula (3.16) for x>0as
follows
fS(x) =K(mX
j=1pjX
0<yx
a+by+
x+
I
yi= 0;8i=1;m;i6=j
fj(yj)fS(x y)
+p0X
0<yx
a+by+
x+
fL(y)fS(x y))
;
and
X
0<yx
a+by+
x+
I
yi= 0;8i=1;m;i6=j
fj(yj)fS(x y)
=xjX
yj=1
a+byj
x+
fj(yj)fS(x1;:::;xj 1;xj yj;xj+1;:::;xm);
which completes the proof.
36
Chapter 4
Alternative methods
Apart from the exact methods (like the recursive one presented above), some approxi-
mate techniques have also been proposed for evaluating aggregate claims distributions
(as we mentioned in Chapter 1), with the purpose to simplify calculations and re-
duce the computing time (for details on these methods see, e.g., [49] or [96]). The
Fast Fourier Transform and Simulation are such techniques that reduce the computing
time, especially when one needs to evaluate the tail of the distribution. Moreover, these
techniques can be applied to models for which there is no recursion available. For that
reason, the FFT for example, received special attention in the actuarial literature, see,
e.g., [12], [26], [43], [49], or [79]. In the following, we shall present the FFT algorithm
and a simulation one corresponding to our model and compare them with the recursive
method.
Futhermore, taking into consideration the fact that the FFT algorithm yields an ap-
proximate distribution, the speci c errors generated by the method are discussed, with
accent on the aliasing error . To reduce the aliasing error, which can be important es-
pecially for heavy-tailed claim sizes distributions, we present a second FFT algorithm
that includes the so-called exponential tilting change of measure. In this chapter, we
numerically illustrate the algorithms on two particular cases of model (3.1), presented
in Chapter 2, with accent on the related errors and their reduction.
37
4.1 Fast Fourier Transform for multivariate claims
This section is based on [79].
4.1.1 Characteristic function
LetEdenote the expected value operator; let 'denote the characteristic function of a
r.v. or vector, and gits probability generating function (p.g.f.), both indexed with the
variable name. We recall that
'S(t) =Eh
eiPm
j=1tjSji
;gN(u) =E"mY
j=1uNj
j#
:
Proposition 4.1 Under the assumptions of model (3.1), it holds that the characteristic
function of Sis given by
'S(t) =gN('L(t);'U1(t1);:::;'Um(tm)): (4.1)
Proof. Based on the independency assumptions on the claim sizes distributions, we
have
'S(t) =Eh
eiPm
j=1tjSji
=E"mY
j=1eitjSj#
=E2
4E2
4mY
j=1eitj NjP
l=0Ujl+N0P
k=0Ljk! N0;:::;Nm3
53
5
=E2
4mY
j=1E
eitjUjNjE"mY
j=1eitjLj#N03
5=E"
'N0
L1;:::;Lm(t1;:::;tm)mY
j=1'Nj
Uj(tj)#
=gN0;:::;Nm('L1;:::;Lm(t1;:::;tm);'U1(t1);:::;'Um(tm));
i.e., formula (4.1). This completes the proof.
In particular, we consider two classical distributions (see [45] for details on them) for
the random vector number of claims, N.
Particular case 1 . Let Nbe multinomial distributed Mnomm+1(n;p0;:::;pm) with
the p.g.f.gN(u) =Pm
j=0pjujn
;wheren2N rf0g;pj2(0;1);0jm;and
Pm
j=0pj= 1. Then, from Proposition 4.1, it holds that
'S(t) = mX
j=1pj'Uj(tj) +p0'L(t)!n
:
38
Particular case 2 . Let Nfollow am+ 1 multivariate Poisson distribution
MPom+1(; 0;:::;m) with>0;j>0;0jm. Its p.g.f. is given by
gN(u) =exp(
mY
j=0uj 1!
+mX
j=0j(uj 1))
:
From Proposition 4.1 we obtain that
'S(t) =exp(
'L(t)mY
j=1'Uj(tj) 1!
+mX
j=1j
'Uj(tj) 1
+0('L(t) 1))
:
For this second particular case, under the same independency assumptions, [76] presents
the following recursive formula
fS(x) =k
xkxkX
zk=0zkfk(zk)fS(x1;:::;xk 1;xk zk;xk+1;:::;xm) +0
xkxX
z=0zkf0(z)fS(x z)
+
xkxX
v=0vk vX
u=0mY
j=1fj(uj)f0(v u)!
fS(x v); xk1; (4.2)
with starting value
fS(0) = exp(
f0(0)mY
j=1fj(0) 1!
+mX
j=1j(fj(0) 1) +0(f0(0) 1))
:
4.1.2 Multidimensional discrete Fourier transforms and FFT
algorithm
Letf(x) be anm-variate function de ned on the integer values xj= 0;1;:::;rj 1;1
jm:Then its discrete Fourier transform (DFT) ~fcan de ned by (de nition used in
Matlab)
~f(c) =r1 1X
x1=0:::rm 1X
xm=0f(x) exp(
2imX
j=1xjcj
rj)
; cj= 0;:::;rj 1;1jm:
Its inverse mapping is given by
f(x) =1
mQ
j=1rjr1 1X
c1=0:::rm 1X
cm=0~f(c) exp(
2imX
j=1xjcj
rj)
; xj= 0;:::;rj 1;1jm:
A FFT is an algorithm that computes the DFT and its inverse extremely fast. There
are many di erent FFT algorithms, one of their most important applications being in
39
signal and image processing, but not only. In order to apply the FFT method, the
valuesrjmust be powers of two for all j:For more details on Fourier transforms and
their applications see, e.g., [10] or [115]; for the use of the FFT method for model (2.1),
see [49]. For model (3.1), the following algorithm based on the FFT and its inverse
(IFFT) can be used to obtain an approximate distribution of S.
Algorithm 1
Step 1. Set the truncation points for the r.v.s claim sizes Ujatrj;1jm;and for L
at (r1;:::;rm). The truncated claim size distributions result as fj=ffj(0);fj(1);:::;fj(rj 1)g
forUj;1jm;andf0= [f0(j1;:::;jm)]j1;:::;jmforL;where 0jlrl 1;1lm.
If necessary, the resulting vectors fjor the table f0can be padded with zeros to force
therjs to be powers of two.
Step 2. Apply the one-dimensional FFT to fjyielding the vector ~fj;1jm; then
apply the multidimensional FFT to f0, yielding the multidimensional table ~f0.
Step 3. Use formula (4.1) to obtain the discrete characteristic function
~'S(j) =gN(~f0(j);~f1(j1);:::;~fm(jm));0jlrl 1;1lm:
Step 4. Apply the multidimensional IFFT to ~ 'Sto obtain the p.f. of S.
Remark 4.2 To nd optimal rjs, one can gradually increase them (e.g., 64, 128, 256
etc.) until the di erences between the solutions obtained for the current values of the
rjs and the previous ones are no more signi cant. More details on errors are given in
the next section.
4.1.3 Types of errors. Exponential tilting
The errors generated by the use of the DFT (and, in particular, of the FFT) have been
investigated mainly in connection with harmonic analysis applications (i.e., images and
signal processing) see, e.g., [42] and [10]. In the insurance eld, even if the calculation
of aggregate claims distributions using the Fourier method starts in 1983 with [33],
only later on, [31] and [32] conducted a thorough study of the related errors and even
proposed an improved FFT procedure based on an exponential change of measure.
Concerning the same problem, [86] noted that \such typical errors can be crucial for
40
the nal result, especially when working with heavy-tailed distributions". In nancial
mathematics, another application of the Fourier transform is related to the valuation
of options. For a detailed study of the corresponding errors, which -depending on the
model- can be important even for distributions with exponentially decaying tails, see,
e.g., [9], [53] and the references therein. Moreover, a comparison of the errors generated
by the DFT method in the context of option pricing versus (vs.) the context of model
(2.1) is realized in Section 6 of [86]. Back to the above Algorithm 1, it essentially
generates two types of errors: discretization and aliasing errors. We shall discuss them
in the sequel.
4.1.3.1 Discretization (arithmetization) errors
For simplicity, we assumed from the very beginning that all the claim sizes have discrete
distributions, which is essential for both FFT and recursive methods. However, in
practical applications, the claims distributions are of continuous type, in which case
they must be rst discretized, i.e., given a discretization (or span) parameter h >
0, the (univariate) claims distribution must be concentrated on the lattice hN0=
f0;h;2h;:::g. The most used discretization procedure is the rounding method (see, e.g.,
[49]), which, in the univariate case, places at jhthe probability fj=F(jh+h=2)
F(jh h=2);j1;whilef0=F(h=2);whereFdenotes the cumulative distribution
function (cdf). The discretized distribution converges weakly to the original one as
h!0, and the same convergence holds for the corresponding compound distributions
of model (2.1) obtained for both continuous and discretized claim sizes.
To provide upper and lower bounds for the true compound distribution of (2.1), [26]
and [86] de ned the so called forward and backward di erences,
f+
j=F((j+ 1)h) F(jh); f
j=F(jh) F((j 1)h):
Starting from these, [26] discussed the following method to choose a suitable span (i.e.,
small enough): his reduced until the di erence between the upper and lower bounds
is smaller than a given error. However, [26] also noticed from numerical examples that
these two bounds are usually too pessimistic, and therefore, they suggested to rather
proceed by successively reducing huntil the (relative) improvement in the correspond-
41
ingly computed compound distributions is small enough.
A detailed study of such discretization errors has also been carried out in [32], which
also discusses some extrapolation techniques to reduce these errors.
Regarding our model (3.1), we might also face the problem of discretizing a multivari-
ate continuous distribution (corresponding to the common claims L). The rounding
method can also be applied in the multivariate case, and for our numerical study we dis-
play it in the trivariate case: letting ( hi)i=1;2;3be the spans, the discretized probability
placed at (j1h1;j2h2;j3h3);ji= 0;1;:::;i = 1;2;3;is
fj1;j2;j3=F(^|1;^|2;^|3) F ~j1;^|2;^|3
F
^|1;~j2;^|3
F
^|1;^|2;~j3
+F ~j1;~j2;^|3
+F ~j1;^|2;~j3
+F
^|1;~j2;~j3
F ~j1;~j2;~j3
;
where ^|i=jihi+hi=2 and ~ji= max (0;jihi hi=2);i= 1;2;3:Upper and lower bounds
can also be provided in this case in a similar manner, but we preferred choosing the
spans according to the second method of [26]. However, we intend to widen the study
of the multivariate rounding errors in another work.
4.1.3.2 Aliasing errors
De nition and upper bound. As requested by the DFT method, the Algorithm 1
starts by truncating the discrete/discretized claim sizes distributions. This truncation
is an important source of errors, generating the so-called aliasing error (AE) ,which
consists of placing below the truncation point the compound mass which lies beyond
this point (i.e., the wrap-around e ect, a DFT speci c error). Based on the approach
of [86], we shall now give a mathematical de nition of the AE associated with our
multivariate model (3.1), and obtain an upper bound for it. Since convolutions are
the main ingredient of (3.3), we start the study with the convolution of two m-variate
functions.
Letf= (fj1:::jm)0jlrl 1
1lmandg= (gj1:::jm)0jlrl 1
1lmbe twom-variate functions with
discrete support. Their convolution c=fgis de ned by
ck1:::km=k1X
j1=0:::kmX
jm=0fj1:::jmgk1 j1:::km jm;0klrl 1;1lm;
42
that we shortly rewrite as
ck=kX
j=0fjgk j;0kr 1:
We now apply the DFT to fgfollowed by its inverse, and denote the result by ^ c. Based
on the fact that the convolution becomes multiplication under the DFT transform, a
straightforward calculation leads to the formula
^ck=r 1X
j=0fjgk j;0kr 1;
where each time a negative index kl jl<0 appears at gk j, it must be replaced with
the indexrl+kl jl(according to the periodicity of the Fourier transform). Clearly,
^ck=ck+AEk(fg), whereAEk(fg) is the corresponding aliasing error. Then we
de ne the overall AE of fgby
AE(fg) =r 1X
k=0AEk(fg);
and note that (as in the univariate case studied in [86]) we can rewrite AE(fg) =
P
Dcfjgk;whereDcconsists of all the multidimensional indexes 0j;kr 1, with
j+kcontaining at least one component jl+klrl; hence,Dc={Dis indeed the com-
plementary of the domain D=f(j;k)j0j;kr 1;j+kr 1g. Now, if f;gare
obtained from two probability density functions f;gasfj1:::jm=f(j1;:::;jm);gj1:::jm=
g(j1;:::;jm), withFandGthe cdf's of fand, respectively, g, then the following
inequality holds
AE(fg)Z
{f(x;y)jx+yr 1gf(x)g(y)dxdy= 1 (FG) (r 1): (4.3)
Clearly, this inequality can be extended to n-fold convolutions. Based on these, we
apply the DFT procedure to our model (3.1). Using also (3.2), for the corresponding
AE we obtain
AE(S) =X
n00;:::;nm0p(n0;:::;nm)AE
G(n1:::nm)Fn0
0
;
whereG(n1:::nm)denotes the cdf of g(n1:::nm);whileF0is the cdf of f0:Then
AE(S)X
n00;:::;nm0p(n0;:::;nm)
1
G(n1:::nm)Fn0
0
(r 1)
;
43
i.e.,
AE(S)1 FS(r 1): (4.4)
This result is consistent with the inequalities obtained by [31] and [86] for the univariate
model (2.1). Still, for this univariate model, in the case when the distribution of the
number of claims has nite variance, [31] deduced a more elaborate formula of the AE,
which supports the conclusion that this error is smaller for larger truncation points,
and when the tail of the claims distribution decays fast.
Exponential tilting. To obtain an exponential decay of the distribution's tail that
reduces the AE (especially when working with heavy-tailed claims distributions), [31]
suggested using a suitable change of measure called exponential tilting (orexponential
window in signal analysis). Mathematically speaking, this method is based on the
fact that, cf. [31], \an exponential change of measure commutes with the compound
distribution functional".
For our multivariate model, with the notation from Algorithm 1, we de ne the tilting
operators
Ejfj=
e jlfj(l)
0lrj 1;1jm; (4.5)
E1;:::;mf0="
exp(
mX
j=1jlj)
f0(l1;:::;lm)#
0ljrj 1;1jm; (4.6)
wherej>0;1jm;are the tilting parameters. To ensure the commutation
with our multivariate compound distribution as discussed above, we must use in the
m-variate operator E1;:::;mthe same tilting parameters as in the univariate operators
(the proof in the bivariate case m= 2 is given in [43], while the general proof is very
similar). Then, including the exponential tilting in Algorithm 1, we obtain the follow-
ing modi ed algorithm:
Algorithm 2 Step 1.Same as in Algorithm 1.
Step 2. Tilt the vectors fj;1jm;and the multidimensional table f0;using (4.5)-
(4.6) (for details on the choice of the tilting parameters see next section).
Step 3. Apply FFT to obtain ]Ejfj;1jm;and^E1;:::;mf0.
Step 4. Use formula (4.1) to obtain the discrete tilted characteristic function of S,
44
which results into a multidimensional table of dimension r1:::rm.
Step 5. Apply IFFT to this multidimensional table, then untilt the result by E 1;:::; m
to obtain the p.f. of S.
We shall now see how this exponential tilting a ects the AE's upper bound (4.4). Since
this exponential tilting commutes with the compound distribution, the p.f. of Sob-
tained by tilting the claim sizes distributions as in (4.5)-(4.6) equals e Pm
j=1jxjfS(x):
Therefore, the overall AE corresponding to Algorithm 2, denoted AETilt
1;:::;m(S);results
in a similar way as AE(S);with the main di erence in inequality (4.3), where, e.g.,
in the univariate case,
AETilt
f2
Z
{f(x;y)jx+yr 1ge xf(x)e yf(y)dxdy =
=Z
f(x;y)jx+y>r 1ge (x+y)f(x)f(y)dxdye (r 1)
1 F2(r 1)
:
Since the exponential tilting also commutes with convolutions (see, e.g., the proofs of
[43]), we therefore obtain
AETilt
1;:::;m(S)e Pm
j=1j(rj 1)(1 FS(r 1)) =e Pm
j=1j(rj 1)AE(S):
This supports the idea that the exponential tilting can produce an important decrease
of the AE by increasing r. However, applying the exponential tilting might cause under-
or over
ow errors, see also the following discussion. Moreover, under- or over
ow errors
might be produced by the p.g.f. of Nin case of \large" frequencies. To solve this issue,
[26] suggested using a portfolio splitting technique.
Reducing the aliasing error. To reduce the AE, the following techniques are rec-
ommended:
– Reasonably increasing the truncation points; if for certain claims distributions there
is no available probability mass between some (smaller) point and the truncation point,
the empty range will be padded with zeros (see rst numerical example);
– Applying the exponential tilting with a careful choice of the tilting parameters. In
practice, the values of these parameters should be large enough, however with spe-
cial care to under- or over
ow errors that might appear at steps 2 and 5 of Algo-
rithm 2. In this sense, considering the truncation points rjand the min/max values
45
exp(
mP
j=1jrj)
, [26] and [31] recommended a rough value ofmP
j=1jrj'20 to avoid
other numerical errors. In the univariate case this gives = 20=r;in the bivariate
case [43] took j= 10=rj;while in our numerical part, we shall use j= 7=rjin the
trivariate case. More details on the choice of the tilting parameters can be found in
[86].
Remark 4.3 In connection with the discrete nature of the Fourier transform algo-
rithm, the numerical integration error must also be considered. Basically, this appears
by replacing the integralRxM
0eitxf(x)dxwith a discrete approximation (here xMis
the right endpoint of the discretization grid). However, since in our case the claims
distributions that we use as input for the FFT are already discretized, and since the
FFT itself does not introduce additional errors, the main error sources remain the ones
described above.
4.1.4 Numerical study
We shall illustrate the use of the FFT in the trivariate case, i.e., for m= 3. Paying
attention to the speci c errors described above, we implemented the method with
and without exponential tilting. In the rst example, we consider only discrete claim
sizes distributions and we illustrate the AE reduction by padding with zeros. The
distribution of the number of claims Nis taken to be the multivariate Poisson, hence,
for comparison, we implemented in Matlab both the exact recursion given in formula
(4.2) in the Particular case 2 of Section 4.1.1, and the FFT method (Matlab already
contains the function tn, which evaluates multidimensional FFTs). Note that the
computing time requested by the exact recursion proved to be so large that we preferred
to evaluate recursively only the values of fS(x) up to some xM= (24;24;24). Hence,
although with the FFT method we obtained much more values of fS, the comparisons
are realized in all examples only on the evaluation support of the recursive method.
46
Example 1 . (AE reduction by padding with zeros) Let
f10
@0 1 2 3
0:3 0:2 0:3 0:21
A; f20
@0 1 2 3
0:4 0:1 0:3 0:21
A; f30
@0 1 2 3
0:2 0:3 0:4 0:11
A;
(f0(0;i;j))i;j=0;1=2
40:15 0:05
0:2 0:13
5;(f0(1;i;j))i;j=0;1=2
40:1 0:2
0:12 0:083
5;
while the number of claims Nis Poisson distributed MPo 4(2;2;3;1;2):As already
mentioned, there is no reason to apply the exponential tilting in this case (more pre-
cisely, we noticed that the best results are for all j= 0), hence we applied only the
Algorithm 1, in which we varied the rjs (taking each time r1=r2=r3) to illustrate the
AE reduction by padding with zeros. In Table 3.1 we display the exact recursive value
FS(xM), the FFT one FFFT
S(xM), the AE and the maximum absolute error between
the exact and the FFT based p.f.s, de ned, respectively, by
AE =X
0xxM fS(x) fFFT
S(x) ;
Max:err = max
0xxM fS(x) fFFT
S(x) ;
where xMis the right endpoint of the support on which fSwas recursively evaluated
(i.e., whenrj= 16;xM= (15;15;15);while forrj32;xM= (24;24;24)). These values
show how increasing the limits rjs improves the results, until that -in this case- taking
rj= 128 yielded about the same results as for 64. The last two lines of Table 3.1
also contain the computing time (shortly named \time"), emphasizing the high speed
gained using the FFT.
Moreover, in Figure 1 we plotted the exact p.f. of the marginal S1and its FFT based
p.f. obtained for rj= 16, which shows the AE cumulated in the left side. For rj= 32,
the FFT based p.f. almost overlaps the exact one.
Since the use of exponential tilting is relevant when the claims distributions are
heavy-tailed, in the following examples we shall consider Pareto distributions of type
II for the claim sizes, having the decumulative distribution (or survival) functions:
– For the univariate Pareto distribution Pa1II( ;); ;> 0;
F(x) = 1 F(x) =
1 +x
;x> 0;
47
rj= 16 rj= 32 rj= 64
xM (15;15;15) (24;24;24) (24;24;24)
FS(xM) 0:8848 0:996969 0:996969
FFFT
S(xM) 1 0:997049 0:996969
AE 0:1152 7:9910 51:2610 14
Max:err 2:7510 42:4910 71:1910 16
total FFT time 0:025sec. 0:104sec. 0:978sec.
recursive
time till xM31:82sec.19min.
Table 4.1: Comparison of recursive and FFT methods for Example 1
Figure 4.1: Recursive and FFT p.f.s of S1for the data in Example 1
48
The expected value of this distribution equals =( 1) and exists only if >1.
– For the trivariate Pareto distribution Pa3II
; (i)i=1;2;3
; ;i>0;i= 1;2;3;
F(x) =
1 +3X
i=1xi
i!
;x>0:
We discretized these distributions using the method of rounding as discussed in Sec-
tion 4.1.3.1. To nd an optimal span h(for simplicity, we used the same span for
all discretized distributions), we proceeded by successively reducing it until the dif-
ferences between the discretized distributions were no longer signi cant (e.g., of order
10 4). For the Pareto distributions considered in Examples 2-4, a good span would
beh= 0:01 (or smaller), but this created computing problems in the trivariate case:
to capture a signi cant mass of probability (over 90%) from the trivariate Pareto dis-
tribution, we needed to evaluate so many discretized values that the computing time
became tremendously long, and the use of the tnMatlab function generated the error
message "out of memory" (on an average computer, these problems appeared from
rj= 29). Unfortunately, these aspects restrain the applicability of the FFT method in
the multivariate setting. This is why in the following three examples we consider larger
values ofh, which, though not optimal, enables us to check at least the eciency of
the exponential tilting.
Regarding this last mentioned aspect, in Algorithm 2 we varied the jvalues and con-
cluded from the following examples that the results improve with the increasing of the
js up to a certain value, in this case, up to j= 7=rjas mentioned in the algorithm;
larger values like j= 9=rjbring no signi cant improvement, however, taking j= 5=rj
worsen the results (see comparisons in the examples). We also varied the (equal) trun-
cation points rjasrj= 64;rj= 128 and rj= 256;the value 128 yielding a good
enough accuracy compared with the exact recursive method up to xM= (24h;24h;24h)
for Examples 2 and 3. Since we did not recursively evaluate fS(x) for larger values of
xdue to the exaggerated computing time, we also present comparisons between the
distributions obtained when varying the rj's; we also used such comparisons as stop-
ping criterion for Example 4, where no recursive formula is available.
49
Example 2 . (AE reduction by exponential tilting) In this example, we considered
the same p.f. f0and the same multivariate Poisson distribution for Nas in Example
1, while the distributions of ( Ui)i=1;2;3were taken to be univariate Pareto, i.e., U1
Pa1II(1;1);U2Pa1II(2;2);U3Pa1II(3;1):Note that the expected value of U1
and the variance of U2do not exist, which enables us to see how the exponential tilting
performs in the heavy-tailed case. In this sense, in Table 3.2 we have a comparison
between the recursive and the two FFT algorithms, in which we varied the parameters
handrj(the formulas of AEandMax:err are the ones de ned in Example 1). All the
displayed cdf's (exact and FFT based) are evaluated till xM= (24h;24h;24h) because
of the recursive method. The results show the improvement brought by the exponential
tilting.
Concerning the computing time, the recursive method up to xMneeded about 20
min. compared with only about 40 sec. required by the FFT method, which evaluates
all values up to (128 h;128h;128h) (note that both computing times do not include the
discretization time); in this case, the exponential tilting did not signi cantly increase
the time compared with no tilting. Moreover, taking rj= 256;the FFT computing
time up to (256 h;256h;256h) increased to about 4 min.
In Table 3.3, we display the maximum probability mass FFFT
S(xM) captured by
the FFT method for h= 0:1 and di erent values of the parameters; note how small is
this value for this discretization span. In Table 3.4, we compare the FFT results when
varying di erent parameters (like the tilting parameters or when increasing the trun-
cation points from 128 to 256) by means of the AEandMax:err de ned in Example
1, but this time based on the p.f.s obtained by the FFT method with the parame-
ters indicated in column 1, evaluated till the smallest xM= (r1h;r2h;r3h). To get an
idea on the discretization accuracy, we also evaluated the FFT cdf for rj= 256 with
exponential tilting with j= 7=rj;by varying h2f0:1;1g. The maximum xMconsid-
ered is (12;12;12), for which FFFT
S;h=1(xM) FFFT
S;h=0:1(xM) = 0:4177 0:3968 = 0:0209;
quite large, meaning that a smaller span his recommended. However, due to time and
memory reasons (see discussion above), we did not try it in the current study.
50
h= 1;xM= (24;24;24) h= 0:1;xM= (2:4;2:4;2:4)
Exact FFT no tilt.FFT tilt.
j= 7=rjExact FFT no tilt.FFT tilt.
j= 7=rj
FS=FFFT
S(xM)0:6879 0:6932 0:6879 0:0206 0:0275 0:0206
AE 5:2710 34:7810 66:8110 35:2810 6
Max:err 7:5210 66:8410 93:9310 63:4110 9
Table 4.2: Comparison of recursive and FFT methods for Example 2; rj= 128
rj= 128;xM= (12:8;12:8;12:8);tilt.rj= 256;xM= (25:6;25:6;25:6)
j= 5=rjj= 7=rjj= 9=rj no tilt. tilt.j= 7=rj
FFFT
S(xM)0:4226 0:4215 0:4213 0:81555 0:6998
Table 4.3: Comparison of probability masses captured by FFT for Example 2; h=
0:1;xM= (r1h;r2h;r3h)
FFT parameters AEbetween FFT p.f.s Max:err
no tilt. vs. tilt. with j= 7=rj,rj= 256 0:1157 1:148110 6
rj= 256 vs.rj= 128;tilt. withj= 7=rj 1:076410 42:384310 9
tilt. withj= 5=rjvsj= 7=rj,rj= 128 4:069010 42:181210 8
tilt. withj= 7=rjvsj= 9=rj,rj= 128 5:498510 52:951010 9
Table 3.4: Comparison of FFT results when varying di erent parameters for Example
2;h= 0:1;xM= (r1h;r2h;r3h)
Example 3 . Considering the same multivariate Poisson distribution for Nand the
same univariate Pareto distributions of ( Ui)i=1;2;3as in Example 2, we now assume that
L= (L1;L2;L3) follows the trivariate Pareto distribution Pa3II(1:5; 2;2;2);which we
also discretized. We proceeded with the same comparisons as in Example 2, the results
being displayed in Tables 3.5-3.7, from where we see a clear improvement brought by
the exponential tilting. Also, we noticed again a strong reduction of the computing
time from about 26 min. for the recursion up to (24 h;24h;24h) to about 42 sec. for the
FFT method up to (128 h;128h;128h) (again, these do not include the discretization
51
time, which now became very long due to the trivariate Pareto discretization).
Regarding the discretization accuracy, we evaluated the FFT cdf for rj= 256
with exponential tilting with j= 7=rj;andh2f0:1;1g. For the maximum xM=
(12;12;12), we obtained FFFT
S;h=1(xM) FFFT
S;h=0:1(xM) = 0:1707 0:1559 = 0:0148;hence
a smaller span his recommended .
h= 1;xM= (24;24;24) h= 0:1;xM= (2:4;2:4;2:4)
Exact FFT no tilt.FFT tilt.
j= 7=rjExact FFT no tilt.FFT tilt.
j= 7=rj
FS=FFFT
S(xM)0:4287 0:4336 0:4287 7:6510 31:0910 27:6510 3
AE 4:8910 34:4210 63:2910 31:6110 6
Max:err 1:4610 61:3310 92:3610 61:9810 9
Table 3.5: Comparison of recursive and FFT methods for Example 3; rj= 128
h= 0:1;rj= 128 h= 1;tilt.j= 7=rj
no tilt. tilt.j= 5=rjj= 7=rjj= 9=rjrj= 64 rj= 128
FFFT
S(xM)0:4461 0:1721 0:1713 0:1712 0:8066 0:9262
xM (12:8;12:8;12:8) (64;64;64) (128;128;128)
Table 3.6: Comparison of probability masses captured by FFT for Example 3;
xM= (r1h;r2h;r3h)
FFT parameters AEbetween FFTs Max:err
no tilt. vs. tilt. with j= 7=rj,rj= 128;h= 1 0:0194 1:464010 6
rj= 64 vs.rj= 128;tilt. withj= 7=rj;h= 1 5:027710 54:994010 9
no tilt. vs. tilt. with j= 5=rj,rj= 128;h= 0:1 0:2739 2:346710 6
tilt. withj= 5=rjvsj= 7=rj,rj= 128;h= 0:1 7:496910 41:263610 8
tilt. withj= 7=rjvsj= 9=rj,rj= 128;h= 0:1 1:838410 51:709610 9
Table 3.7: Comparison of FFT results when varying di erent parameters for Example
3;xM= (r1h;r2h;r3h)
Example 4 . In this example, we considered the same Pareto distributions for
(Ui)i=1;2;3and for Las in Example 3. However, we changed the distribution of Nto be
52
multinomial Mnom 4(10; 0:25;0:3;0:25;0:2) as de ned in the rst particular case from
Section 2.1. Note that for this model we do not have an exact recursive formula, hence
we directly implemented Algorithm 2 for rj= 64 and then for rj= 128;j= 1;2;3;with
j= 7=rjin both cases. The di erences AEandMax:err based on the corresponding
p.f.s (and evaluated till the smallest xM= (r1h;r2h;r3h)) being very small (see Table
3.8), we conclude that, at least for these data, the choice rj= 128 would be enough.
More comparisons of the FFT results when varying di erent parameters are shown in
Tables 3.9-3.10. The computing time for rj= 64 was about 5 sec., while for rj=
128 about 41 sec. (discretization time not included). In this case also, an accuracy
improvement would probably be generated by a smaller span.
h= 1 h= 0:1
rj= 64 rj= 128 rj= 64 rj= 128
xM= (r1h;r2h;r3h)(64;64;64) (128;128;128) (6:4;6:4;6:4)(12:8;12:8;12:8)
FFFT
S(xM) 0:9019 0:9606 0:0803 0:3386
AEbetween FFTs 1:665910 55:817910 5
Max:err 8:731810 93:999410 10
Table 3.8: Comparison of FFT methods for Example 4 when varying rj;j= 7=rj
h= 0:1;rj= 128;xM= (12:8;12:8;12:8);tilt./no tilt.
no tilt.j= 5=rjj= 7=rjj= 9=rj
FFFT
S(xM)0:5863 0:3397 0:3386 0:3385
Table 3.9: Comparison of probability masses captured by FFT for Example 3
FFT parameters AEbetween FFTs Max:err
no tilt. vs. tilt. with j= 7=rj,rj= 128;h= 1 5:776310 32:779410 6
no tilt. vs. tilt. with j= 5=rj,rj= 128;h= 0:1 0:2465 4:652410 7
tilt. withj= 5=rjvsj= 7=rj,rj= 128;h= 0:1 1:072810 32:094910 9
tilt. withj= 7=rjvsj= 9=rj,rj= 128;h= 0:1 1:449210 42:831110 10
Table 3.10: Comparison of FFT results when varying other parameters for Example
4;xM= (r1h;r2h;r3h)
53
4.2 Simulation for multivaritate claims
4.2.1 Comparison between simulation and recursive methods
For comparison reasons, we also implemented the simulation and FFT methods for the
same data. Some results of the FFT method can be found in [78], where the main
advantage of the FFT -its high speed- is emphasized by presenting several computing
times (e.g., to recursively evaluate all the exact values up to fS(31;31;31);19.66 seconds
were needed, while the FFT method did the same in only 0.19 seconds, with a very
good accuracy).
Regarding the simulation method, we implemented in Matlab an algorithm that
simulates model (3.1) starting with the simulation of a Poisson value of N, followed by
the simulation of the multinomial values of N0;N1;N2;N3;and, based on these values,
we simulated corresponding numbers of claims distributed according to fL;f1;f2;f3;
then, by summing these claims values according to model (3.1), we obtained a simulated
value of S= (S1;S2;S3). By repeating this process Ktimes (Klarge enough), we
obtained the simulated values of fSthat we compared with the recursive exact ones by
calculating the largest absolute error value between the simulated and exact p.f.s till
the point r= (16;16;16), i.e.,
Err(r) := max
0xr fS(x) fSimulated
S (x) :
This error is displayed in Table 1 together with some values of the cumulative distri-
bution function of S,FS. Note that the error di erences between performing 5 105
simulations and 106simulations are not so important, hence we can conclude that for
these data, the simulation method based on 106works well enough.
Regarding the computing time, simulation with K= 106took about 6 min, while
withK= 5105about 3 min. This might seem long, but note that by simulation
we obtain approximate values for a very large domain of fS;and to obtain exact
recursive values for the same domain it would take an excessively long time (e.g., for
x= (64;64;46) it needed about 23 min).
54
Err(16;16;16)FS(3;3;3)FS(8;8;8)FS(16;16;16)
Exact values 0.4454 0.9658 0.99995
Number of K= 1062:228310 40.4456 0.9656 0.99994
simulations K= 51052:669710 40.4450 0.9660 0.99994
K= 1056:970610 40.4436 0.9660 0.99997
Table 3.11: Comparison between the recursive and simulation methods
4.2.2 Simulation method-Matlab implementation
1clc
2clear all
3M=100000; lbd =5; p1 =0.25; p2 =0.3; p3 =0.2;
4f=zeros (30 ,30 ,30) ;
5for i =1:M
6 n=poissrnd ( lbd ) ;
7 i fn==0
8 s=zeros (3) ;
9 else
10 nn=zeros (3) ; n0=0;
11 for j =1:n
12 u=rand ;
13 i fu<=p1
14 nn (1)=nn (1) +1;
15 elseif u<=p1+p2
16 nn (2)=nn (2) +1;
17 elseif u<=p1+p2+p3
18 nn (3)=nn (3) +1;
19 else n0=n0+1;
20 end
21 end
55
22 for j =1:3
23 su ( j )=sumaU(nn( j ) , j ) ;
24 end
25 sL=sumaL( n0 ) ;
26 for j =1:3
27 s ( j )=su ( j )+sL ( j ) ;
28 end
29 end
30 f ( s (1) +1, s (2) +1, s (3) +1)=f ( s (1) +1, s (2) +1, s (3) +1)+1/ M;
31end
32cdf =0;
33for i 1 =1:4
34 for i 2 =1:4
35 for i 3 =1:4
36 cdf=cdf+f ( i1 , i2 , i 3 ) ;
37 f ( i1 , i2 , i 3 )
38 end
39 end
40end
41cdf
1function s=sumaL(n)
2ss=zeros (3) ;
3for i =1:n
4 u=rand ;
5 i fu<=0.15
6 ss (1)=ss (1) +0;
7 elseif u<=0.25
8 ss (3)=ss (3) +1;
9 elseif u<=0.3
10 ss (2)=ss (2) +1;
56
11 elseif u<=0.5
12 ss (3)=ss (3) +1; ss (2)=ss (2) +1;
13 elseif u<=0.7
14 ss (1)=ss (1) +1;
15 elseif u<=0.82
16 ss (1)=ss (1) +1;
17 ss (3)=ss (3) +1;
18 elseif u<=0.92
19 ss (1)=ss (1) +1;
20 ss (2)=ss (2) +1;
21 else ss (1)=ss (1) +1;
22 ss (2)=ss (2) +1;
23 ss (3)=ss (3) +1;
24 end
25end
26 s=ss ;
1function s=sumaU(n , j )
2i fj==1
3 f j =[0.3 0.2 0.3 0 . 2 ] ;
4elseif j==2
5 f j =[0.4 0.1 0.3 0 . 2 ] ;
6else f j =[0.2 0.3 0.4 0 . 1 ] ;
7end
8ss =0;
9for i =1:n
10 u=rand ;
11 i fu<=f j (1)
12 ss=ss +0;
13 elseif u<=f j (1)+f j (2)
14 ss=ss +1;
57
15 elseif u<=f j (1)+f j (2)+f j (3)
16 ss=ss +2;
17 else ss=ss +3;
18 end
19end
20s=ss ;
58
Chapter 5
Monte Carlo method – a particular
study case
A well-known method such as Monte Carlo simulation is quite often used to analyze
the risks in a project development. This chapter aims to present the method using an
ongoing ship design project for a petroleum chemical tanker.
First of all it is important to know that the primary focus during the development of a
ship basic/detailed design is to pay attention to an important aspect called risk. Taking
into consideration that building a ship implies a lot of risks, an important duty is to
prevent them by maximizing the probability and consequences of positive events and,
in the same time, by minimizing the probability and consequences of adverse events
related to the project's objectives. Shipbulding industry is for centuries a worldwide
subject, see for example Hoving [39], Chida and Davies [15], Briggs Vernon [11] and
involves a lot of aspects such as those that we will present in the following. According
to studies regarding the evolution of naval constructions, evidenced in Market Research
Reports, the industry of naval constructions has registered a series of major
uctuations
lately. Ship orders have constantly risen during a four year period, from 2003 to
2007 and plummeted in 2008-2009, as a consequence of the global nancial crisis.
After a slight comeback in 2010, the naval industry starts experiencing problems again
in 2011. The slow evolution of ship building, worldwide speaking, re
ects the weak
economical foundations, the debt crisis within the euro zone and the lack of loans
59
in maritime transport. The ship building industry is faced with a dark perspective
due to the constant growth of naval shipyards eciency, modest investments in the
eld and also the market saturation. Building new ships means substantial costs, a
great consumption of raw matter and materials, signi cant human resources as well
as launching new dedicated techniques and technologies. Consequently, to start a new
project, a risk analysis is required to re
ect as accurate as possible the implications and
consequences of all the factors which are about to participate in the task achievement.
5.1 Risk analysis basis
Risk analysis is not speci c only to insurance, but it is also involved in many other elds
see for example Hillerbrand [37]; we also mention elds such as: insurance and business,
Grandell [30] and Bhlmann [13]; nance: Murphy [65] and Sortino [91]; engineering:
Ayyub [6], McCormick [61], Modares et al. [63]; ecology: Scholze et al. [87] and Lash
et al. [52]; medicine: Levitan et al. [54] and Vincent et al. [97]; agriculture: Radulescu
et al.,[74], shing: Radulescu et al.[73], queueing theory: Raducan et al. [72] and other
industry: Akintoye and MacLeod [1], Khan and Abbasi [47], Ramiro and Aisa [75], Wu
et al. [116].
Designing a new ship involves the identi cation and assuming multiple risks, dependent
on the environment, the mobility of the production process, long execution terms, the
large number of coordinated activities and the complexity of the entire process as
well as contractual relations with third parties. The risk is the possibility of an event
occurring, usually an unfavorable one, bearing consequences on the cost performances
(it usually leads to further expense), on the quality (not rising up to the necessary
quality speci cations) and on time (not meeting deadlines). In other words, risk can
be de ned as uncertainty associated with a desired result. Risk can appear in three
situations:
the event occurs, but the result is uncertain;
the e ect of an event is known, but the result is uncertain;
the event is uncertain as well as the result.
60
Risk analysis uses the theory of probabilities and it consists of studying all variables
of an objective which can contribute to the obtaining of acceptable performances, con-
sidering all the factors already discussed. This analysis can be made quantitatively,
through stochastic methods by quantifying the possible value dispersion as previsions
in terms of cost and/or time, and from a qualitative point of view, through a formal
investigation, based on judgments, experience and diagnostic methods using checklists.
The risk element in a project is represented by any element which possesses a measur-
able probability to deviate from a well-de ned plan. Any element of such a structure
can be used as a potential risk element. The more representative a structure is, the
more signi cant the risk elements will be taken into consideration. Thus, to achieve
a goal, a series of activities must be run. An activity, (a), can be considered a risk
element if two simultaneous conditions are met :
8
><
>:0<P(a)<1
L(a) = 0(5.1)
Here,P(a) represents the probability that the event ( a) will occur; E(a) represents
the e ect of event ( a) on the project; L(a) represents the nancial evaluation of E(a).
Generally, evaluating risks means covering the following stages (see [82] ):
sensitivity analysis;
critical variables probability distribution;
risk assessment;
evaluating acceptable levels of risk;
risk prevention.
5.2 Risk management aspects
Risk management is a cyclic process which unfolds throughout the entire project and
it consists of covering at least three distinct stages, as it can be observed in Table 1.
61
No.Risk management
stagesActivities
1. RISK IDENTIFICATION- identify the risk
– identify signi cant risks
– eliminate incongruous risks
2. RISK ANALYSIS- in-depth analysis of signi cant risks
– result assessment
-probability assessment
-expected value determination
3. REACTION TO RISK- risk reduction
– risk elimination
– risk distribution
Table 5.1: Risk management stages
In the risk identifying stage one must evaluate potentials perils, as well as e ects
and apparition probabilities in order to decide which speci c risks need prevention.
Actually, at this stage, one must identify all the elements which meet the criteria ( ??).
At the same time, one must eliminate incongruous risks that are those elements with
low probability of occurrence or with insigni cant consequences. This means that we
shall not take into consideration those elements whose P(a) or L(a) value tend to zero.
In the risk analysis stage one must take into consideration the risks primarily identi ed
and another in-depth analysis must be drawn up. To perform a risk analysis, one must
use various mathematical instruments, ranging from probability analysis to Monte
Carlo analysis. In the reaction to risk stage one takes action to eliminate, to reduce or
to distribute risks. These actions are carried out according to a risk management plan
which includes the procedures to be used in risk management, designating individuals
responsible for certain actions in risk domains, resources assigned to this purpose, as
well as assessing results obtained in risk management.
62
5.3 Risk analysis methods
5.3.1 Qualitative methods
The qualitative method of risk analysis implies using qualitative criteria, using various
categories for separating parameters, with de nitions establishing the range for every
category. Qualitative analysis is based on expert opinions, it doesnt require specialized
instruments and it is easily applied. Interpreting results is clear to most people, not
requiring special knowledge, and the conclusions are sucient. Qualitative risk anal-
ysis elements consist of risk probability, risk impact and risk assessment matrix (see
Ionescu 2010). The risk probability is de ned as the probability of that risk occurring.
According to the probabilities theory, the probability of the event (a) is de ned as the
ratio of the number of favorable results (m) of the event (a) and the total number of
the event results (n), considered equally possible, that is :
P(a) =m
n(5.2)
Risk impact points to the e ect that risk has upon the project objectives, if it occurs.
The risk probability and impact, I(a), are seen as very high, high, medium, low or very
low. Risk score, S(a) is determined using the formula
S(a) =P(a)(a) (5.3)
To determine the actual scores of every risk, one must draw up the risk score matrix,
as can be seen in Table 2. The signi cance of the data in Table 2 is as follows:
S(a)<0;05 = low impact (green);
0;05S(a)<0;15 = medium impact (yellow);
S(a)>0;15 = high impact (red).
The qualitative approach is subjective but it allows a higher generalization degree,
being less restrictive.
63
Risk impact
Probability
0,50 0,10 0,20 0,40 0,80
0,90 0,05 0,09 0,18 0,36 0,72
0,70 0,04 0,07 0,14 0,28 0,56
0,50 0,03 0,05 0,10 0,20 0,40
0,30 0,02 0,03 0,06 0,12 0,24
0,10 0,01 0,01 0,02 0,04 0,08
Table 5.2: Risk score matrix
5.3.2 Quantitative methods
The quantitative risk analysis implies the evaluation of the risk value through numerical
methods. Through quantitative risk analysis one takes into consideration the numerical
evaluation of probability and impact of each risk upon the entire project. To this end,
several quantitative methods are employed, such as decision trees analysis or the Monte
Carlo simulation. The quantitative analysis cannot be made without the qualitative
analysis to determine the initiating events. The method requires a high volume of
data, very well prepared professionals and expensive software. Decision trees method
represents decisions and random events, in the way they are viewed by the decision
makers. As represented in Fig. 1 ([41]), each probable future event (represented by
a circle) is matched with an action (a square) which can be taken by the decision
maker, resulting in a tree-like graph. The actual occurrence of a future event may
mean that one of more actions needs to be taken, and the sum of all probabilities
equal to unity. The signi cance of notations in Fig. 1 is as follows: D decision
point; E probable event; U pro t according to various uncertainty degrees. The
Monte Carlo analysis, elaborated in 1940, uses repeated random sampling to generate
simulated data to use with a mathematical model (see [82] and [51]). In this context,
the simulation is represented by the approximation process of a result through the
random and repetitive use of the algorithm. Monte Carlo simulation allows project
managers to incorporate uncertainty and risk in project planning unlike other risk
analysis methods which cannot quantify these aspects equally well. The Monte Carlo
64
Figure 5.1: Decision tree
simulation is frequently applied, due to its acknowledged advantages, both by the
practitioners (see [ ?]) and by the academic community (see [98]). The quantitative
approach is much more objective and accurate, but one must mention the fact that
results can be a ected by the precision and the validity of input parameters. Due to
this fact, the quantitative results mustnt be considered exact numbers, but estimative
ones, with a variable scale depending on the quality of the data used in the analysis.
5.4 Monte Carlo simulation for a ship design project
In this study case we consider a short part of a ship design project which was planned
to start on January 4, 2016 and scheduled to be completed on July 18, 2017 at a cost
of $120.000.
5.4.1 Preliminaries
The planning of the project was developed in the program software Primavera P6,
having initially a level 3 work breakdown structure (WBS) as in gure 2. After the
planning was made and the Baseline was assigned, the next step was to export our
project into a software able to help us to make a risk analysis. In consequence, we
imported it in Primavera Risk Analysis (Fig. 3). In order to start an analysis we
65
Figure 5.2: Basic design project schedule in Primavera P6
Figure 5.3: Basic Design project schedule imported to Primavera Risk Analysis
66
must assure that aspects such as: nish date, duration and cost are well completed
and imported from P6. First of all we will we evaluate the risk only under the e ect of
duration uncertainty using two distribution types: triangle and uniform. The triangular
distribution is commonly used to model the duration of a task. Its simple set of
parameters make it easy to relate to real life. It requires minimum, maximum and
most likely duration, having the following density
f(x) =8
>>><
>>>:2(x a)
(b a)(c a);ax<c
2(b x)
(b a)(b c);c<xb
0;otherwise(5.4)
This distribution is often skewed to the left. This is because a lot of tasks cannot
physically be completed in less than a certain duration, but all tasks can generally be
delayed for any number of reasons. This leads to the minimum duration being closer
to the most likely than the maximum duration. A long thin long tail on the triangular
distribution models a range of things that could go wrong but are unlikely. The Uniform
distribution is used when only the extremes of uncertainty can be speci ed and where
intermediate values have equal probability of occurrence (or where no inference can be
made as to the likely shape of the distribution). Its density is given by
f(x) =8
<
:1
(b a);ax<b
0;otherwise(5.5)
The software allows us to simulate all mentioned above using also other distribution
types such as: Normal, Beta Pert, Lognormal, Discrete, Trigged, cumulative, etc. In
the second part of the simulation the project will be constrained by a risk register. All
these issues will be detailed later.
5.4.2 Simulation 1
The percentage assigned to the minimum, likely and maximum are applied to the total
duration/ activity. After running the Monte Carlo simulation, we can observe (Fig.
5.4) that the possibility of nishing on time is 90%. Also, we can observe that the P-80
date is now July 04, 2017 for a Triangle distribution, compared with May 30, 2017 for
the Uniform distribution.
67
WBS NameDistribution
typeMin. Likely Max.
Structural design documents Triangle 25% 75% 100%
Ships equipment, out tting Triangle 70% 105% 120%
Accommodation Triangle 85% 110% 130%
Table 5.3: Set uncertainties to activities
Figure 5.4: Triangle distribution nish date
68
Figure 5.5: Uniform distribution nish
From Fig. 5.5, we can observe that if we choose a uniform distribution the proba-
bility to nish in time is more optimistic, 98%. This is due to the shape of the Uniform
distribution, but we would expect a Triangular distribution to be more realistic for
such a study.
Figure 5.6: Triangle distribution cost
69
Figure 5.7: Uniform distribution cost
The costs for Simulation 2 are almost similarly in both cases: even though we use
a Uniform or a Triangle distribution (Fig.5.6 and Fig.5.7), there are 94% chances for
the Uniform distribution and 96% chances in case we chose a Triangle distribution to
maintain the initial budget.
5.4.3 Simulation 2
In this part, as we mentioned before, we use a Risk Register (Fig. 5.9). Every risk is
de ned by four attributes such as: probability, schedule, cost, performance and score.
Every attribute has assigned a factor scale (low, medium, high) in accordance with its
project impact factor. The Matrix Risk Register (Fig. 5.10) shows the before Pre-
Mitigation period.
We assume risks such as:
Changing the input data
70
Figure 5.8: Risk Register
Leaving employees
Class certi cation
Logistic database
Figure 5.9: Matrix Risk Register
71
Figure 5.10: Risk result for Finish date Monte Carlo simulation
We can easily observe that in this case, the P-80 date is now August 18, 2017 or almost
more than one month from the scheduled date of July 18, 2017, and that the possibility
of nishing on time is in this scenario only 39% (Fig. 11). Regarding the costs, there
72
Figure 5.11: Cost – with the default Risk Register
are only 19% chances to nish the project with the initial budget ($120.000), and a
P-80 posibility to have a cost of $282.448, which means a $162.488 amount di erence.
(Fig. 5.12)
5.4.4 Risk Register- Cost and Duration sensitivity
We can also make a classi cation regarding the Cost (Fig. 5.13) and Duration (Fig.5.14)
sensitivity. It results that the Leaving employees risk is the most sensitive with a 85%
possibility to happen if we take into consideration the cost, while Class certi cation is
the most dangerous risk if we look at the duration issue. By adding default risk to
73
Figure 5.12: Cost – Cost
Figure 5.13: Cost – Duration
74
Item Initial Pre-mitigation
Duration 562 d 586 d
Cost $120.000 $383.283
Start Date 04/01/2016 04/01/2016
Finish Date 18/07/2017 11/08/2017
Table 5.4: Result of analysis
the task activities, the following aspects changed (see Table 4):
After a Pre-Mitigation process, a Post-Mitigated process should follow. Furthermore,
in order to improve the possibility of having a project succes, it is indicated – and the
experience has shown this- that it a risk analysis must exist periodically throughout
all stages of a project, namely: project planning, engineering, procurement, execution
and commissioning. Finally, in order to improve the perspective of a success, Hullet
and Avalon [40] recommended a plan which implies the following stages:
Create a summary schedule at Level 2 or Level 3 that is integrated, logically-
linked, and complies with scheduling best practices. The budget or cost esti-
mate is incorporated into the schedule with summary time-dependent and time-
independent resources assigned to activities.
Gather data about uncertainty and risk. Experience has shown that individual
con dential interviews of subject matter experts in the project teams and man-
agement, contractors, and other knowledgeable personnel in the organization, can
uncover risks that are either new or unpopular to discuss risk workshops.
Use a modern Monte Carlo simulation software package to model both uncertain-
ties to activity duration and to time-dependent resourcess burn rates or time-
dependent resources total material costs.
Prioritize the Risk Drivers to schedule and to cost (only schedule was illustrated
above).
Prepare a pre-mitigated risk result presentation.
75
Conduct a risk mitigation workshop using the prioritized risk contained in the
pre-mitigation presentation as a guide to the most important risks.
Gain organizational commitment to those mitigation actions to be adopted,
scheduled, budgeted, sta ed, and monitored.
Prepare a post-mitigation model and simulate it.
Prepare a post-mitigation report.
76
Chapter 6
Conclusions
6.1 General conclusions
In this PhD thesis original results have been presented in accordance with the study of
the aggregate claims distributions of a portfolio insurance, portfolio which may contain
policies for: res,
oods, earthquakes, trac accidents, and so on.
We focused rstly on the extension of two bivariate models to a multivariate setting,
and moreover we created our own multivariate aggregate claims model, with its own
dependencies.
Our goal was to obtain an exact recursive formula for the probability function of the
multivariate compound distribution of the both models.
Why did we do this, and did we aim by doing this? -the answer is obvious, we wanted to
contribute to the improvement of the actuarial science in order to assist the actuaries to
a better analysis of an insurance scenario, especially when a portfolio is simultaneously
a ected by di erent types of claims. It is important to mention here that all these
aspects were proved not only in classical mathematical ways but also with numerical
examples.
Furthermore, due to the fact that I personally activate as a planner engineer in the
shipbuilding industry I consider that it was a good idea to complete my thesis with a
study case from this eld of activity. I wanted to emphasize the strong link between
planning and the mathematical science, and arm on that way that : "Math is every-
77
where!".
I showed how important a simulation method such as Monte Carlo could be and how
this can help us make a relevant risk analysis of a ship design project.
In my opinion, I believe that making this mixture between the classical stochastic the-
ories and real life circumstances make my research interesting to read and undoubtedly
original.
It is very important to mention that the entire research of these thesis is the result
of a strong collaboration with my scienti c advisor, PhD Professor Vasile Preda and
with PhD Associate Professor Raluca Vernic who, due to their academic experience,
coordinated my ideas in the best way possible in order to make a good and ecient
work.
6.2 Future research perspectives
Actuarial science is a eld that give us many research opportunities. In other words
we could say that it challenges us to discover and develop many other aspects from
that perpective. For that reason in the future we are interested in approaching the
following issues:
Study how we can determine the distribution of aggregate claims with a so called
moment-based technique, in the literature; [44] studied this technique in the
univariate and the bivariate cases, for instance we propose to see if we could do
it in a multivariate way. Other moment-based approaches of distributions are:
Provost [71], Lindsay et al. [56], Nadarajah et al. [58], and [7].
Study in more detail the FFT method in the multivariate case, i.e., the discretiza-
tion method and the corresponding errors.
Introduce new multivariate models with new and more complex type of depen-
dencies between the claims.
Introduce dependency between the number of claims and the claim sizes.
78
We would also like to study some numerical aspects related to recursions, like
under
ow/over
ow and stability. Waldmann [99], Waldmann [100] and Wang and
Panjer [107] discussed the under
ow/over
ow question related to the recursive
evaluation of compound distributions. Even if the recursive method is an exact
one, the stability issue appears due to the rounding errors generated by the
computer. Stability of recursions for distributions has been considered by [106],
[108], [107], [109],[108] and [27]. Therefore, we would also like to investigate such
aspects for our multivariate models.
Comparative studies with other methods, like simulation.
Personally I am interested in developing in parallel a university career and for
that reason I intend to develop more study cases as the one presented in my
thesis and after that to implement it in a teaching way in order to persuade the
students that mathematics has applications in various elds of activity and once
they discover how to do that the problem becomes simpler.
79
Bibliography
[1] Akintoye, A. S., MacLeod, M. J. (1997). Risk analysis and management in con-
struction. International journal of project management, 15(1), 31-38.
[2] Ambagaspitiya, R.S. (1995). A family of discrete distributions. Insur. Math. Econ.
16, 107-127.
[3] Ambagaspitiya, R.S., Balakrishnan, N. (1994). On the compound generalized Pois-
son distribution. ASTIN Bull. 24, 255-263.
[4] Ambagaspitiya, R.S. (1999). On the distributions of two classes of correlated ag-
gregate claims, Insur. Math. Econ. 24, 301-308.
[5] Anastasiadis, S., Chukova, S. (2012). Multivariate insurance models: an overview.
Insurance: Mathematics and Economics, 51(1), 222-227.
[6] Ayyub, B. M. (2014). Risk analysis in engineering and economics. CRC Press.
[7] Bee, M. (2015). Density Approximations and VAR Computation for Compound
Poisson-lognormal Distributions. Communications in Statistics-Simulation and
Computation, (to appear).
[8] Binder, K. (1986). Introduction: Theory and technical aspects of Monte Carlo
simulations. In Monte Carlo Methods in Statistical Physics (pp. 1-45). Springer
Berlin Heidelberg.
[9] Boyarchenko, S. and Levendorskii, S. (2014). Ecient variations of the Fourier
transform in applications to option pricing. The Journal of Computational Finance
18 (2), 57-90.
80
[10] Briggs, W.L. and Henson, V.E. (1995). The DFT: An Owners' Manual for the
Discrete Fourier Transform. Siam.
[11] Briggs Vernon, L. (2015). History of Shipbuilding on North River, Plymouth
County, Massachusetts: With Genealogies of the Shipbuilders, and Accounts of
the Industries Upon Its Tributaries, 1640 to 1872 (Classic Reprint), Forgotten
Books.
[12] B uhlmann, H. (1984). Numerical evaluation of the compound Poisson distribution:
recursion or fast Fourier transform? Scand. Actuar. J., 116-126.
[13] B uhlmann, H. (2007). Mathematical methods in risk theory (Vol. 172). Springer
Science Business Media.
[14] Charpentier, A. (2015). Computational Acturial Science with R, CRC Press.
[15] Chida, T., Davies, O. (2013). The Japanese Shipping and Shipbuilding Industries:
A History of their Modern Growth,Bloomsbury Academic Collections.
[16] De Pril, N. (1986). On the exact computation of the aggregate claism distribution
in the individual life model, ASTIN Bull. 16, 109-112.
[17] De Pril, N. (1989). The aggregate claims distribution in the individual model with
arbitrary positive claims, ASTIN Bull. 19, 9-24.
[18] Dedu, S.,Ciumara, R. (2010). Restricted Optimal Retention in Stop-Loss Reinsur-
ance, Proceedings of The Romanian Academy, 11, 3, pp. 213-217
[19] Denuit, M., Dhaene, J., Goovaerts, M., Kaas, R. (2006). Actuarial Theory for
Dependent Risks: Measures, Orders and Models, John Wiley Sons.
[20] Dhaene, J., Vandebroek, M. (1995). Recursions for the individual model, Insur.
Math. Econ. 16, 31-38.
[21] Dlugokecki, V., Fanguy, D., Hepintstall, L.(2009). Leading the way for mid-tier
shipyards to implement design for production methodologies, Journal of Ship Pro-
duction, Vol. 25, No. 2, p. 99-108.
81
[22] Dufresne, F. (1996). An extension of Kornya's method with application to pension
funds, Bull. Swiss Assoc. Actuar., 171-181.
[23] Dorp van J.R. ; Du ey, M. R. (1997). Statistical dependence in risk analysis for
project networks using Monte Carlo methods, Elsevier International.
[24] Dunn, W. L., Shultis, J. K. (2011). Exploring Monte Carlo Methods. Elsevier.
[25] Eisele, K. T. (2006). Recursions for compound phase distributions. Insur. Math.
Econ. 38, 149-156.
[26] Embrechts,P., Gr ubel R. and Pitts, S. M. (1993). Some applications of the fast
Fourier transform algorithm in insurance mathematics. This paper is dedicated to
Professor WS Jewell on the occasion of his 60th birthday , Stat. Neerl.,59-75, 47.
[27] Gerhold, S., Schmock, U., Warnung, R.,(2008, in press). A generalization of Pan-
jer's recursion and numerically stable risk aggregation. Finance Stoch.
[28] Glasserman, P. (2003). Monte Carlo methods in nancial engineering (Vol. 53).
Springer Science Business Media.
[29] Goovaerts, M.J., Kaas, R. (1991). Evaluating compound generalized Poisson dis-
tributions recursively, ASTIN Bull. 21, 193-198.
[30] Grandell, J. (2012). Aspects of risk theory. Springer Science Business Media.
[31] Grubel, R. and Hermesmeier, R. (1999). Computation of compound distributions
I: aliasing errors and exponential tilting. ASTIN Bulletin 29 (2), 197-214.
[32] Grubel, R. and Hermesmeier, R. (2000). Computation of compound distributions
II: discretization errors and Richardson extrapolation. ASTIN Bulletin 30 (2),
309-331.
[33] Heckman, P. E. and Meyers, G. G. (1983). The calculation of aggregate loss dis-
tributions from claim severity and claim count distributions. Proceedings of the
Casualty Actuarial Society LXX, 22-61.
82
[34] Hesselager, O. (1994). A recursive procedure for calculation of some compound
distributions ASTIN Bull. 24, 19-32.
[35] Hesselager, O. (1996). A recursive procedure for calculation of some mixed com-
pound Poisson distributions. Scandinavian Actuarial Journal, 1, pp. 54-63.
[36] Hesselager, O. (1997). Recursions for a class of compound Lagrangian distribu-
tion. Bull. Swiss Assoc. Actuar., 95-101 Correction: Bull. Swiss Assoc. Actuar.,
135(1998).
[37] Hillerbrand, R., Sandin, P., Peterson, M., Roeser, S. (2012). Handbook of Risk
Theory: Epistemology, Decision Theory, Ethics, and Social Implications of Risk,
Springer-Verlag.
[38] Horn, R.A., Steutel, F.W.(1978). On multivariate in nitely divisible distributions.
Stoch. Process. Their Appl. 6, 139-151.
[39] Hoving, A.J.,Wildeman, D.(2012). Nicolaes Witsen and Shipbuilding in the Dutch
Golden Age,Texas AM University Press.
[40] Hulett,D.T., Avalon, A.(2015). Integrated Cost and Schedule Risk Analysis, Long
International Inc.
[41] Ionescu (Mandru),L.(2010). Theoretical research and case studies on risk assess-
ment in an integrated management system for quality-risk industrial companies.
PhD Thesis Summary, Transilvania University of Brasov.
[42] Jerri, A.(1992). Integral and discrete transforms with applications and error anal-
ysis (Vol. 162). CRC.
[43] Jin, T.,Ren, J. (2014). Recursions and fast Fourier transforms for a new bivariate
aggregate claims model. Scandinavian Actuarial Journal, 8, pp. 729-752.
[44] Jin, T., Provost, S.B.,Ren, J. (2016). Recursions and fast Fourier transforms for
a new bivariate aggregate claims model. Scandinavian Actuarial Journal, 3, pp.
216-245.
83
[45] Johnson, N.L., Kotz, S., Balakrishnan, N. (1997). Discrete Multivariate Distribu-
tions. Wiley, New York.
[46] Johnson, N.L., Kemp, A. W., Kotz, S. (2005). Univariate Discrete Distributions,
3rd edn. Wiley-Interscience, Hoboken.
[47] Khan, F. I., Abbasi, S. A.(2001). Risk analysis of a typical chemical industry using
ORA procedure. Journal of Loss Prevention in the Process Industries, 14(1), 43-59.
[48] Kitano, M., Shimizu, K., Ong, S.H. (2005). The generalized charlierseries distri-
bution as a distribution with two-step recursion. Stat. Probab. Lett. 75, 280-290.
[49] Klugman, S.A.,Panjer, H.H.,Willmot, G.E. (2004). Loss Models: From Data to
Decisions, 2nd edn, Wiley-Interscience. Hoboken.
[50] Kolic, D., Fafandjel, N., Calic, B.(2010).Determining how to apply the design for
production concept in shipyards through risk analysis, Engineering Review, 30-1,
63-72.
[51] Kroese, D.P., Taimre, T. Botev, Z. I. (2013). Handbook of Monte Carlo methods
(Vol. 706). John Wiley Sons.
[52] Lash, S., Szerszynski, B., Wynne, B. (Eds.)(1996). Risk, environment and moder-
nity: towards a new ecology (Vol. 40). Sage.
[53] Levendorskii, S.(2013). Pitfalls of the Fourier transform method in ane models,
and remedies. Available at SSRN 2367547.
[54] Levitan, N., Dowlati, A., Remick, S. C., Tahsildar, H. I., Sivinski, L. D., Beyth,
R., Rimm, A. A. (1999). Rates of initial and recurrent thromboembolic disease
among patients with malignancy versus those without malignancy: risk analysis
using Medicare claims data. Medicine, 78(5), 285-91.
[55] Lin, P., Neil, M., Fenton, N. (2014). Risk aggregation in the presence of discrete
causally connected random variables. Annals of Actuarial Science, 8(02), 298-319.
84
[56] Lindsay, B. G., Pilla, R. S., Basak, P. (2000). Moment-based approximations of
distributions using mixtures: Theory and applications. Annals of the Institute of
Statistical Mathematics, 52(2), 215-230.
[57] Liu, L., Cheung, E. C. (2015). On a bivariate risk process with a dividend barrier
strategy. Annals of Actuarial Science, 9(01), 3-35.
[58] Nadarajah, S., Chu, J., Jiang, X. (2016). On moment based density approxima-
tions for aggregate losses. Journal of Computational and Applied Mathematics,
298, 152-166.
[59] Mata, A.J.(2000). Pricing excess of loss reinsurance with reinstatements. ASTIN
Bull. 30, 349-368.
[60] Mata, A.J.(2003). Asymptotic dependence of reinsurance aggregate claim
amounts. ASTIN Bull. 33, 239-263.
[61] McCormick, N. J. (1981). Reliability and risk analysis: methods and nuclear power
applications (pp. 221-223). San Diego, CA:: Academic Press.
[62] Minkova, L. D., Balakrishnan, N. (2014). On a bivariate Plya-Aeppli distribution.
Communications in Statistics-Theory and Methods, 43(23), 5026-5038.
[63] Modarres, M., Kaminskiy, M. P., Krivtsov, V. (2009). Reliability engineering and
risk analysis: a practical guide. CRC press.
[64] Murthy, K. P. N.(2003). Theoretical Studies Section, Materials Science Division,
Indira Gandhi Centre for Atomic Research, Kalpakkam 603 102, Tamil Nadu
INDIA.
[65] Murphy, D. (2008). Understanding risk: The theory and practice of nancial risk
management. CRC Press.
[66] Panjer, H. H. (1981). Recursive evaluation of a family of compound distributions,
ASTIN Bull., 12(1), 22-26.
85
[67] Panjer, H. H., Willmot, G.E.(1982). Recursions for compound distributions.
ASTIN Bull. 13, 1-11.
[68] Panjer, H.H. (1990). The aggregate claims distribution and stop-loss reinsurance.
Transactions of the Society of Actuaries, 32, pp. 523-545
[69] Preda, V., Panaitescu, E.,Ciumara, R.(2011). The modi ed exponential-Poisson
distribution. Proceedings of the Romanian Academy, 12, 1, pp. 22-29.
[70] Promislow, S. D. (2014). Fundamentals of Actuarial Mahtematics, John Wiley
Sons.
[71] Provost, S. B.(2005). Moment-based density approximants. Mathematica Journal,
9(4), 727-756.
[72] Raducan, A. M., Lakatos, L. A. S. Z. L. O., Zbaganu, G. (2008). Computable
Lindley processes in queueing and risk theories. Rev. Roumaine Math. Pures Appl.,
53(23), 239-266.
[73] Radulescu, M., Radulescu, C. Z., Turek Rahoveanu, M., Zbaganu, G. (2010).
A portfolio theory approach to shery management. Studies in Informatics and
Control, 19(3), 285-294.
[74] Radulescu, M., Radulescu, C. Z. and Zbaganu, G. (2014). A portfolio theory
approach to crop planning under environmental constraints. Annals of Operations
Research, 219(1), 243-264.
[75] Ramiro, J. S., Asa, P. B. (2012). Risk analysis and reduction in the chemical
process industry. Springer Science Business Media.
[76] Robe-Voinea, E.G. and Vernic, R. (2015a). On a multivariate aggregate claims
model with multivariate Poisson counting distribution. To appear in Proceedings
of the Romanian Academy – Series A.
[77] Robe-Voinea, E.G. and Vernic, R.(2015b). On the recursive evaluation of a certain
multivariate compound distribution. 8th Congress of Romanian Mathematicians,
June 28-July 1, 2015, "Alexandru Ioan Cuza" University of Iasi.
86
[78] Robe-Voinea, E.G. and Vernic, R. (2015b). Another approach to the evaluation of
a certain multivariate compound distribution. To appear in Analele S tiinti ce ale
Universitatii Ovidius
[79] Robe-Voinea, E.G. and Vernic, R.,(2016). Fast Fourier Transform for multivariate
aggregate claims, Computational and Applied, Springer International Publishing ,
DOI 10.1007/s40314-016-0336-6, Print ISSN 0101-8205, Online ISSN 1807-0302.
[80] Robe-Voinea, E.G. and Vernic, R. (2015c). On the recursive evaluation of a certain
multivariate compound distribution. To appear in Acta Mathematicae Applicatae
Sinica (English Series), Springer International Publishing.
[81] Robert, C., Casella, G. (2009). Introducing Monte Carlo Methods with R.
Springer Science Business Media.
[82] Roman, M., Andreica, M. (2012). Paper no. 10- Developing risk analysis into the
cost-bene t analysis of projects nanced by the ERDF and CF, Contract "Ca-
pacity development for cost-bene t analysis", project co- nanced by the ERDF
through OPTA 2007-2013, Bucharest, University of Economics Studies.
[83] Rubinstein, R. Y., Kroese, D. P. (2011). Simulation and the Monte Carlo method
(Vol. 707). John Wiley Sons.
[84] Sundt, B., Vernic, R.(2004). Recursions for compound mixed multivariate Poisson
distributions. Bltter der DGVFM, 26, 4, pp. 665-691.
[85] Sundt, B.,Vernic, R.(2009). Recursions for Convolutions and Compound Distribu-
tions with Insurance Applications, EAA Lectures Notes. Springer-Verlag Berlin.
[86] Schaller, P. and Temnov, G. (2008). Ecient and precise computation of con-
volutions: applying t to heavy tailed distributions. Computational Methods in
Applied Mathematics 8 (2), 187-200.
[87] Scholze, M., Knorr, W., Arnell, N. W., Prentice, I. C. (2006). A climate-change
risk analysis for world ecosystems. Proceedings of the National Academy of Sci-
ences, 103(35), 13116-13120.
87
[88] Sharif, A.H., Panjer, H.H. (1995). An improved recursion for the compound gen-
eralised Poisson distribution. Bull. Swiss Assoc. Actuar., 93-98.
[89] Sharif, A.H. (1996). Generalized recursions for total claims distributions. Disser-
tation, University of Waterloo.
[90] Sharif, A.H., Panjer, H.H. (1998). Stepwiserecursions for a class of Lagrangian
distributions. Actuar. Res. Clearing house, 333-364.
[91] Sortino, F. A., Satchell, S. (Eds.). (2001). Managing downside risk in nancial
markets. Butterworth-Heinemann.
[92] Storch, R.L., Lim, S.(1999).Improving
ow to achieve lean manufacturing in ship-
building, Production Planning and Control, Vol. 10, No.2, p.127-137.
[93] Sundt, B. (1992). On some extensions of Panjer's class of counting distributions.
ASTIN Bulletin 22, 1, 61-80.
[94] Sundt, B.(1999). On multivariate Panjer recursions. ASTIN Bulletin, 29, 1, 29-45.
[95] Sundt, B.(2000). Multivariate compound Poisson distributions and in nite divisi-
bility. ASTIN Bull. 30, 305-308.
[96] Sundt, B. and Vernic, R.(2009). Recursions for convolutions and compound distri-
butions with insurance applications. EAA Lectures Notes, Springer-Verlag Berlin.
[97] Vincent, C., Taylor-Adams, S., Stanhope, N. (1998). Framework for analysing
risk and safety in clinical medicine. Bmj, 316(7138), 1154-1157.
[98] Vose, D. (2008). Risk Analysis: A Quantitative Guide, 2nd Edition, Wiley.
[99] Waldmann, K.H. (1994). On the exact calculation of the aggregate claims distri-
butions in the individual life model, ASTIN Bull. 24, 89-96.
[100] Waldmann, K.H. (1995). Exact calculation of the aggregate claims distribution
in the individual life model by use of an n-layer model. Bl atter Dtsch. Ges. ver-
sicher.math. 22, 279-287.
88
[101] Walhin, J.F., Paris, J. (2000). The e ect of excess of loss reinsurance with re-
instatements on the cedent's potfolio, Bl atter Dtsch. Ges. Versicher.math. 24,
615-627.
[102] Walhin, J.F., Paris, J. (2001). Excess of loss reinsurance with reinstatements:
premium calculation and ruin probability of the cedent, Bl atter Dtsch. Ges. Ver-
sicher.math. 25, 1-12.
[103] Walhin, J.F., Herfurth, L., De Longueville, P. (2001). The practical pricing of
excess of loss treaties: actuarial, nancial, economic and comercial aspects. Belg.
Actuar. Bull. 1, 40-57.
[104] Walhin, J.F. (2002a). On the optimality of multiline excess of loss covers. Belg.
Actuar. Bull. 2, 92-96.
[105] Walhin, J.F. (2002b). Some comments on the pricing of an exotic excess of loss
treaty. J. Actuar. Pract. 10, 175-191.
[106] Wang, S., Panjer, H.H.(1993). Critical starting points for stable evaluation of
mixed Poisson probabilities, Insur. Math. Econ. 13, 287-297.
[107] Wang, S., Panjer, H.H. (1994). Proportional convergence and tail-cutting tech-
niques in evaluating aggregate claim distributions, Insur. Math. Econ. 24, 161-166.
[108] Wang, S. (1995). On two side compound binomial distributions, Insur. Math.
Econ. 17, 35-41.
[109] Wang, S., Sobrero, M. (1994). Further results on Hessalager's recursive procedure
for calculation of some compound distributions. ASTIN Bull. 24, 161-166.
[110] Wang, R., Peng, L., Yang, J. (2015). CreditRisk+ Model with Dependent Risk
Factors. North American Actuarial Journal, 19(1), 24-40.
[111] Widouck, S., Walhin, J.F. (2004). Some remarks on thecedent's retention risk in
presence of an annual aggregate deductible or reinstatements. Bl atter Dtsch. Ges.
Versicher.-Finanzmath. 26, 461-481.
89
[112] Williams, R.E.(1980). Computing the probability function for aggregate claims.
Proceedings of the Canadian Institute of Actuaries, XI, pp. 39-47.
[113] Willmot, G.E.(1986). Mixed compound Poisson distributions. ASTIN Bulletin,
16, pp. 59-79.
[114] Winston, W.L.(2004) Introduction to Probability Models Operations Research,
Vol. 2, 4th edition, Thomson Learning, Canada.
[115] Wong,M. W. (2011). Discrete Fourier Analysis, Springer, Basel.
[116] Wu, T., Blackhurst, J., Chidambaram, V. (2006). A model for inbound supply
risk analysis. Computers in industry, 57(4), 350-365.
[117] Zbaganu, G. (2004). Metode matematice in teoria riscului si actuariat, Editura
Universitatii, Bucureti.
[118] Zbaganu,G. (2007). Elemente de teoria ruinei, Geometry Balkan Press,
Bucharest.
90
Copyright Notice
© Licențiada.org respectă drepturile de proprietate intelectuală și așteaptă ca toți utilizatorii să facă același lucru. Dacă consideri că un conținut de pe site încalcă drepturile tale de autor, te rugăm să trimiți o notificare DMCA.
Acest articol: FACULTY OF MATHEMATICS AND COMPUTER SCIENCE PhD Thesis Stochastic models in actuarial science. Theoretical contributions and applications Scienti c… [614330] (ID: 614330)
Dacă considerați că acest conținut vă încalcă drepturile de autor, vă rugăm să depuneți o cerere pe pagina noastră Copyright Takedown.