Rainer Kolisch
email: rainer.kolisch@wi.tum.de

Christian Heimerl
email: christian.heimerl@wi.tum.de

Work assignment to and qualification of multi-skilled human resources under knowledge depreciation and company skill level targets

Keywords: Project staffing, Resource assignment, Multi-skilled resources, Learning, Forgetting

Qualification of multi-skilled human resources under knowledge depreciation and company skill level targets

Assigning project work to human resources (henceforth called resources) taking into account resource-specific skills and efficiencies is a general planning task which has to be performed in any organization. It is of particular importance for service firms where the labour intensity is high (cf. e.g. Brusco [START_REF] Brusco | An exact algorithm for a workforce allocation problem with application to an analysis of cross-training policies[END_REF]). The multi-skill assignment problem is an important variant of this planning problem. It requires work packages with a certain demand for one out of several skills to be assigned to resources with multiple skills (cf. Heimerl and Kolisch [START_REF] Heimerl | Scheduling and staffing multiple projects with a multi-skilled workforce[END_REF]). Because of its importance considerable effort has been devoted to address this problem class, and various models and solution methods have been proposed in literature (e.g. cf. Alba and Chicano [START_REF] Alba | Software project management with GAs[END_REF], Ballou and Tayi [START_REF] Ballou | A decision aid for the selection and scheduling of software maintenance projects[END_REF], Bassett [START_REF] Bassett | Assigning projects to optimize the utilization of employees' time and expertise[END_REF], Bellenguez and Néron [START_REF] Bellenguez | Lower bounds for the multi-skill project scheduling problem with hierarchical levels of skills[END_REF], Gutjahr et al. [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF]).

Concerning assignment decisions previous research followed two different approaches: Several authors (e.g. [1, 6, 7, 9-11, 13, 15, 26, 27, 33-35, 39]) assume a discrete 0-1 assignment of resources to tasks. This approach is appropriate if e.g. transfer times between tasks (e.g. due to travelling) are relatively high and cannot be neglected (cf. Krüger and Scholl [START_REF] Krüger | Managing and modelling general resource transfers in (multi-)project scheduling[END_REF]), the lengths of the periods are relatively short, or the transfer between tasks must occur concurrently (e.g. due to a minimum number of resources per task constraint). Other publications (cf. e.g. [2, 4, 18-21, 23, 42]) allow (continuous or discrete) fractional assignments of resources to tasks. With this assumption it is possible to model multi-tasking of resources which is a realistic assumption for period lengths of weeks or months. While from a mathematical point of view, the continuous problem is easier solve, the interpretation of the fractional assignment results within a period remains as a subsequent problem to solve.

Usually there are company-and/or resource-specific goals that assignment problems want to achieve: Company-specific goals can be the maximization of the company's revenue [START_REF] Abboud | Manpower allocation using genetic annealing[END_REF][START_REF] Kwak | Human resource allocation in a CPA firm: A fuzzy set approach[END_REF], the maximization of the company's utility [START_REF] Campbell | Development and evaluation of an assignment heuristic for allocating cross-trained workers[END_REF][START_REF] Sayin | Assigning cross-trained workers to departments: A two-stage optimization model to maximize utility and skill improvement[END_REF], the maximization of the company's performance or productivity [START_REF] Chu | A manpower allocation model of job specialization[END_REF][START_REF] Kwak | Human resource allocation in a CPA firm: A fuzzy set approach[END_REF][START_REF] Nembhard | Learning and forgetting-based worker selection for tasks of varying complexity[END_REF], the minimization of costs for (external) resources [START_REF] Bassett | Assigning projects to optimize the utilization of employees' time and expertise[END_REF][START_REF] Bhatnagar | Contingent manpower planning in a high clock speed industry[END_REF][START_REF] Cai | A genetic algorithm for scheduling staff of mixed skills under multi-criteria[END_REF][START_REF] Eiselt | Employee positioning and workload allocation[END_REF][START_REF] Wu | A project scheduling and staff assignment model considering learning effect[END_REF][START_REF] Wu | On the manpower allocation within matrix organization: A fuzzy linear programming approach[END_REF] and/or overtime [START_REF] Bhatnagar | Contingent manpower planning in a high clock speed industry[END_REF][START_REF] Eiselt | Employee positioning and workload allocation[END_REF], the minimization of project duration [START_REF] Bellenguez | Lower bounds for the multi-skill project scheduling problem with hierarchical levels of skills[END_REF][START_REF] Vairaktarakis | The value of resource flexibility in the resource-constrained job assignment problem[END_REF] or the optimization of other performance indicators [START_REF] Hendriks | Human resource allocation in a multi-project R&D environment: Resource capacity allocation and project portfolio planning in practice[END_REF]. Typical resource-specific goals are the maximization of the resources' satisfaction [START_REF] Abboud | Manpower allocation using genetic annealing[END_REF][START_REF] Yoshimura | Decision-making support system for human resource allocation in product development projects[END_REF], the maximization of preferences or priorities for a task [START_REF] Corominas | Multi-objective allocation of multifunction workers with lower bounded capacity[END_REF], the maximization of appropriateness for a task [START_REF] Ballou | A decision aid for the selection and scheduling of software maintenance projects[END_REF], the minimization of the resources' boredom [START_REF] Eiselt | Employee positioning and workload allocation[END_REF], the maximization of skill improvement [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF][START_REF] Sayin | Assigning cross-trained workers to departments: A two-stage optimization model to maximize utility and skill improvement[END_REF], equalizing workload among resources [START_REF] Cai | A genetic algorithm for scheduling staff of mixed skills under multi-criteria[END_REF][START_REF] Eiselt | Employee positioning and workload allocation[END_REF] or compliance with soft-constraints like a limit on the maximum number of similar tasks per resource [START_REF] Corominas | Rotational allocation of tasks to multifunctional workers in a service industry[END_REF].

Only few of these aforementioned models consider learning effects of resources (cf. e.g. Gutjahr et al. [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF], Wu and Sun [START_REF] Wu | A project scheduling and staff assignment model considering learning effect[END_REF]). Instead, most models assume static and often homogeneous efficiencies of resources. However, efficiencies of resources are usually dynamic due to the acquisition of knowledge, such that more experienced resources have learned and therefore have higher efficiencies. Learning effects can be observed for repetitive tasks in manual, cognitive and knowledge-based work (cf. e.g. Arzi and Shtub [START_REF] Arzi | Learning and forgetting in mental and mechanical tasks: A comparative study[END_REF], Boh et al. [START_REF] Boh | Learning from experience in software development: A multilevel analysis[END_REF], Nembhard and Uzumeri [START_REF] Nembhard | Experiential learning and forgetting for manual and cognitive tasks[END_REF]). Boh et al. [START_REF] Boh | Learning from experience in software development: A multilevel analysis[END_REF] show that learning effects are present in special project environments with knowledgebase work, too. Wright [START_REF] Wright | Factors affecting the cost of airplanes[END_REF] introduced learning curves as a common way to describe learning processes and an extensive amount of literature deals with different types of learning curves (cf. e.g. Nembhard and Uzumeri [START_REF] Nembhard | Individual-based description of learning within an organization[END_REF], Yelle [START_REF] Yelle | The learning curve: Historical review and comprehensive survey[END_REF] for an overview). On the other hand there is depreciation of knowledge. This effect is considered by even fewer assignment models (cf. e.g. Gutjahr et al. [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF], Nembhard [START_REF] Nembhard | Heuristic approach for assigning workers to tasks based on individual learning rates[END_REF], Süer and Tummaluri [START_REF] Süer | Multi-period operator assignment considering skills, learning and forgetting in labour-intensive cells[END_REF]). Depreciation of knowledge can be caused by internal and/or external effects. Internal effects can depict e.g. forgetting, external effects can represent e.g. technological progress, which renders some older experience useless in the new environment or requires learning of new concepts in order to maintain efficiency (cf. Chen and Edgington [START_REF] Chen | Assessing value in organizational knowledge creation: Considerations for knowledge workers[END_REF]).

In this paper we address the problem of assigning project work to resources. We assume that project schedules are given and that project work can be aggregated across individual projects such that we obtain a resource demand for each skill and period of e.g. one month length. The aggregation enables the use of continuous assignment decisions. Learning and depreciation of knowledge of resources is taken into account and in this context we try to achieve two company-specific goals: The operative goal is to minimize costs for performing a given amount of project work during the planning horizon. On a strategic level, however, we need to decide who is applying and improving which skills in order to obtain given company skill level targets measured in production rates. The latter goal is expressed as a constraint. Of course, assignment problems can be applied to arbitrary-sized organizational units, e.g. companies, departments, or groups. W.l.o.g. we will refer to the organizational unit of a company for the remainder of this paper.

To the best of our knowledge, the aforementioned problem has only been scarcely discussed in the literature. Nembhard [START_REF] Nembhard | Heuristic approach for assigning workers to tasks based on individual learning rates[END_REF] solves the problem of assigning resources to skills with assignment policies but does not consider depreciation of knowledge or strategic goals. Wu and Sun [START_REF] Wu | A project scheduling and staff assignment model considering learning effect[END_REF] apply a meta-heuristic to the scheduling and staffing problem of projects w.r.t. learning but they neglect depreciation of knowledge and strategic goals, as well.

Our approach has many similarities with Gutjahr et al. [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF]. Both models use multi-skilled resources and apply the concept of learning and forgetting. However, Gutjahr et al. [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF] employ an objective function maximizing the weighted sum of strategic goals (learning perspective) and portfolio value while we minimize costs in the objective function and use constraints in order to accomplish the strategic goals. We neither consider project selection nor scheduling, but take into account outsourcing decisions. Finally, in the experimental studies we solve the model with non-linear learning curves as opposed to the linear approximations used in Gutjahr et al. [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF].

The remainder of this paper is organized as follows: The following Section 2 introduces an optimization model for the assignment of resources to project work. In Section 3 we describe implementation issues and our test setup. Section 4 presents computational results and gives managerial insights to the impact of different prob- lem parameters on the costs, outsourcing decisions and specialization measures of resources and the company. Finally, the paper concludes with a summary and an outlook in Section 5.

Model

We refer to a work package (s, t) with the demand r st as the aggregated amount of project work requiring skill s ∈ S in period t = 1, . . . , T , where T denotes the planning horizon. The pool of resources consists of internal resources R i and external resources R e , i.e. R = R i ∪ R e . The demand of work packages requiring skill s can be assigned to resources k ∈ R s = R i s ∪ R e s with cost rates c kt in period t. R s is the set of resources capable of performing project work which requires skill s. The set of skills that resource k possesses is denoted with S k . Working times of resources are limited by time-dependent availabilities R kt .

The projects are broken up into skill-specific work packages (s, t), which -in contrast to the projects itself -are expected to have a repetitive character. Thus, the related learning process is well-suited to be described by learning curves.

Learning curves are usually monotonically decreasing and convex functions. A learning curve f ks (z ks ) describes the unit production time, i.e. the amount of time that resource k requires to produce one additional unit after having produced z ks units using skill s. The number of units z ks produced by resource k in skill s can be interpreted as the experience of the resource w.r.t. the particular skill. We index the function f ks with resource k and skill s since the learning process depends on the type of work s and the abilities of resource k to adapt new knowledge. Furthermore, the argument z ks is also indexed with k and s, i.e. we neglect cross-skill or team learning effects.

The time F ks (z ks ) that resource k requires to produce the first z ks units using skill s can be expressed by the integral F ks (z ks ) = z ks z =0 f ks (z )dz . The time τ ks required to produce x ks units after having produced z ks -x ks units can then be calculated by τ ks = F ks (z ks ) -F ks (z ks -x ks ).

We define z kst as the experience level of resource k in skill s at the end of period t. Resource k's amount of project work performed using skill s in period t is represented by x kst . We model resource k's amount of depreciated knowledge β kst in skill s as a loss of experience at the beginning of period t. Given an initial experience level z ks0 of resource k in skill s at the end of period t = 0, z kst can be calculated by z kst = z ks(t-1) -β kst + x kst . As stated earlier β kst can depict e.g. forgetting and/or technological progress: For example release dates of new technologies which are known in advance can be modelled due to the index t. Using index s the impact on experience can be modelled skill-dependent. E.g. the release of a new programming language might affect programming skills but not database skills. With index k we can discriminate resources w.r.t. their forgetting properties. Note that z kst can become negative. However, w.l.o.g. z kst ≥ 0 can be assured by shifting the learning curve T t=1 β kst units to the right and setting z ks0 := z ks0 + T t=1 β kst . This approach is similar to dynamic inventory models (cf. e.g. Zipkin [START_REF] Zipkin | Foundations of Inventory Management[END_REF]) with z kst being the current inventory of knowledge, β kst being the demand or loss of knowledge and x kst being the ordered and instantaneously delivered (or produced) knowledge.

Strategic goals might require to keep or develop certain skills across the company. To define strategic goals we employ the production possibility frontier (PPF, also known as transformation curve). The PPF is a concave curve describing maximum production quantities for different skills given several (but limited) resources (cf. Varian [START_REF] Varian | Intermediate Microeconomics: A Modern Approach[END_REF]). Note that in our assignment problem the PPF is not static or given in advance but depends on the actual assignments as it is transformed over time by learning and depreciation of knowledge of human resources. In our model we will not employ requirements on the shape of the complete PPF but only on the values of its intercept points at the axes at the end of the planning horizon. We require the intercept point at axis s to be at least φ s . Thus, the level φ s is the guaranteed total amount of work per time unit requiring skill s that internal human resources could perform at the end of the planning horizon T if they all used that skill. Figure 1 gives an example for the possible result of the PPF at the end of the planning horizon. A different assignment decision would yield a different PPF. W.r.t. the given values of φ 1 and φ 2 this solution is feasible, since production rates for both skills could attain level φ s independently (but not concurrently). The intercept points of the PPF can be calculated as the sum of production rates of the internal resources with skill s at the end of the planning horizon T defined by Note that the production rate is the reciprocal of the unit production time calculated by f ks (•).

k∈R i s 1 f ks (z ksT ) . ( 1
We model the outlined problem as a non-linear program by employing the following decision variables: The amount of work done by resource k with skill s in period t is denoted with x kst . Resource k's experience in skill s cumulated up to period t considering depreciation of knowledge is depicted with z kst . Finally, τ kst is the time that resource k requires to process the amount of work x kst of skill s in period t.

Z = Min s∈S k∈R s T t=1 c kt • τ kst (2)
subject to

z kst = z ks(t-1) -β kst + x kst k ∈ R s s ∈ S t = 1, . . . , T (3)
τ kst = F ks (z kst ) -F ks (z kst -x kst ) k ∈ R s s ∈ S t = 1, . . . , T (4)
k∈R i s 1 f ks (z ksT ) ≥ φ s s ∈ S (5) k∈Rs x kst ≥ r st s ∈ S t = 1, . . . , T (6)
s∈S

τ kst ≤ R kt k ∈ R t = 1, . . . , T (7)
x kst , τ kst ≥ 0

k ∈ R s s ∈ S t = 1, . . . , T (8)
z kst ∈ R k ∈ R s s ∈ S t = 1, . . . , T (9)
The objective function (2) minimizes the costs for performing the work packages (s, t). Constraints (3) are the dynamic experience level constraints, where the experience level z kst of resource k in skill s at the end of period t is the experience level of the preceding period (t -1) decreased by depreciation of knowledge β kst and increased by the acquisition of knowledge x kst due to work assignments in period t. Constraints (4) calculate the time τ kst required to perform x kst units of work package (s, t). Constraints (5) enforce the company skill level targets φ s at the end of the planning horizon. Constraints [START_REF] Bellenguez | Lower bounds for the multi-skill project scheduling problem with hierarchical levels of skills[END_REF] ensure that the demand r st of work package (s, t) is assigned to the resources. Continuous fractional assignment of work packages to resources is allowed. Constraints [START_REF] Bhatnagar | Contingent manpower planning in a high clock speed industry[END_REF] restrict the availability of the resources. Constraints ( 8) and ( 9) define the decision variables.

If there is no learning and no forgetting, model ( 2)-( 9) becomes linear with the learning function f ks (z ks ) = 1 η ks , where η ks is the static efficiency of resource k w.r.t. skill s (cf. e.g. Heimerl and Kolisch [START_REF] Heimerl | Scheduling and staffing multiple projects with a multi-skilled workforce[END_REF]).

Although the learning functions f ks (z ks ) are assumed to be decreasing and F ks (z ks ) is therefore concave, the model is neither a convex nor a concave optimization problem since constraints (4) calculate the difference of two concave functions resulting in an unclassifiable function.

Note that in Gutjahr et al. [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF] the experience ("competence score" in [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF]) of resources increases due to the amount of time invested into a skill. In contrast, we assume that the amount of work done by a resource drives its experience which determines the production rate ("efficiency values" in [START_REF] Gutjahr | Competence-driven project portfolio selection, scheduling and staff assignment[END_REF]).

Implementation and Test setup

For our experimental study we employed an adaptation of the exponential learning function (cf. Nembhard and Uzumeri [START_REF] Nembhard | Individual-based description of learning within an organization[END_REF], Pendharkar and Subramanian [START_REF] Pendharkar | An empirical study of icase learning curves and probability bounds for software development effort[END_REF]), i.e.

f ks (z ks ) = a ks e -λ ks z ks + b ks [START_REF] Cai | A genetic algorithm for scheduling staff of mixed skills under multi-criteria[END_REF] and therefore

F ks (z ks ) = a ks λ ks 1 -e -λ ks z ks + b ks z ks . ( 11
)
b ks > 0 represents the steady state unit production time, λ ks ≥ 0 is the learning rate and a ks ≥ 0 the learning potential of resource k in skill s. Note that knowledge depreciation is depicted by β kst in constraints (3) and is therefore implicitly included in the learning curve. We chose this type of learning function due to the ability to depict steady state unit production times and its mathematical tractability. However, the implementation can easily be adapted to other learning functions.

The model was implemented in C++ using Ipopt 3.3.2 of the COIN-OR library (cf. [START_REF]Coin-OR[END_REF]). Ipopt uses a primal-dual interior point filter line search algorithm (cf. Wächter and Biegler [START_REF] Wächter | On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming[END_REF]) to return an (at least locally) optimal solution. We solved the model 50 times using random starting points in order to partially overcome the lack of only local optimality.

We generated test instances using the base case depicted in Table 1. The values are chosen to obtain realistic small size instances with period lengths of one month and efforts measured in man days. A set of external resources R e s was taken into account by adding one additional resource with unlimited availabilities for each skill, i.e. s . These numbers are based on the assumption that there will always be enough qualified external resources with state-of-the-art knowledge. The costs of external resources depend on the requested skill, and we refer to skill 1 as the cheapest and skill 4 as the most expensive skill.

All six internal resources are capable of performing two out of four skills and each resource has a unique combination of skills. A learning rate of λ ks = 0.012 together with a ks = 0.2 and b ks = 0.9 leads to a unit production time of 1 after approximately 60 work units. The parameters of |R i |, |S k |, r st , R kt were chosen such that the utilization of internal resources is 150% if no external resources are employed and all internal resources had static unit production times of 1, i.e. f ks (•) ≡ 1. Furthermore, despite learning effects, the utilization will always be higher than 100%. We varied four parameters to generate the 16 instances given in Table 2. The first group of instances (1-8) uses a flat learning curve (λ ks = 0.012), the other group (9-16) uses a steeper learning curve (λ ks = 0.02). The shape of the learning curve affects the speed of learning in the sense that resources with a steeper learning curve learn faster. Figure 2 depicts the plots of the learning curves. Within each group some instances (1-4, 9-12) neglect forgetting (β kst = 0), the remaining instances impose forgetting (β kst = 10). Some instances (1-2, 5-6, 9-10, 13-14) do not demand company skill level targets (φ s = 0), while all the other instances do (φ s > 0). Furthermore, we compare cases with internal costs c i k = 0 (odd instance numbers) and c i k = 500 (even instance numbers). For the latter instances internal resources are more expensive than external resources w.r.t. skill 1 (c i k = 500 vs. c e 1 = 400) in case both internal and external resources had equal and horizontal learning curves (i.e. f ks (•) ≡ 1).

T = 6 |S| = 4 |R i | = 6 R kt = 20 ∀t, k ∈ R i r st = 45 ∀s, t a ks = 0.2 ∀s, k ∈ R i s b ks = 0.9 ∀s, k ∈ R i s z ks0 = 0 ∀s, k ∈ R i

Computational results

The test instances led to a model with 96 variables and 70 constraints. They were solved on a Pentium 4 with 2.4 GHz and 1.5 GB RAM. We recorded the required time to solve 50 random starting point problems (cf. Table 3). For the instances given in Table 2 the solution time is 21 seconds on average and always less than 33 seconds. The instances with internal costs c i k = 500 are harder to solve than those with c i k = 0 (24.6 vs. 17.0 seconds). This is probably due to the fact that the full utilization of internal resources is not always an optimal solution in cases where c i k = 500 and thus constraints [START_REF] Bhatnagar | Contingent manpower planning in a high clock speed industry[END_REF] are not always binding. The benefit of using random starting points is measured using the average and maximum deviation of different starting point solutions to the best solution found. For 5 out of 16 instances (31%) the use of random starting points led to no improvements. On the other hand, the average distance from the locally optimal solutions to the best solutions found was 0.24% and the maximum distance was 8.3%. Thus, for some instances choosing a good starting point matters and random starting points are a straight-forward strategy to find a good solution.

# time [ 4. Obviously, even large instances with 100 resources and 50 skills can be solved within minutes using the 50 random starting points approach. However, the solutions did not benefit from the use of random starting points. x R ks be the amount of work done by internal (company) resources for skill s during the planning horizon, then

|R i |
ρ s = x C
s t r st is the fraction of work ρ s done by internal resources for skill s during the planning horizon. We use ρ s to analyse outsourcing decisions. Further let gives an indication about the specialization of the company: A high coefficient of variation indicates high specialization of the company, since the assignment of internal resources to work packages of different skills is very unequal, i.e. high for some skills, low for others. The values of ρ s and CV(X C ) for the solutions of the instances are listed in Table 5.

X R k = {x R k1 , x R k2 , x R k3 , . . .} and X C = {x C 1 , x C 2 ,
The specialization of individual resource k ∈ R i can be measured with

σ k = σ k σ max k
which uses the standard deviation

σ k = VAR(X R k ) = E (X R k ) 2 -(E (X R k )) 2 = s (x R ks ) 2 |S k | - s x R ks |S k | 2
and is normalized with the maximum possible standard deviation σ max k . σ max k can be determined by calculating σ k if only one skill is used as much as possible. If all skills are used equally by resource k, i.e. x R ks = x R ks ∀s, s , we obtain σ k = 0. If only one skill is used, σ k = 1 holds. The values for the specialization σ k are also listed in Table 5.

# ρ 1 ρ 2 ρ 3 ρ 4 CV(X C ) σ 1 σ 2 σ 3 σ 4 σ 5 σ 6
avg. 1 0.00 0.62 1.00 1.00 0.156 1.00 1.00 1.00 1.00 0.37 0.48 0.81 2 0.00 0.44 1.00 1.00 0.172 1.00 1.00 1.00 1.00 1.00 0.17 0.86 3 0.16 0.46 1.00 1.00 0.138 1.00 1.00 0.40 1.00 0.92 0.48 0.80 4 0.16 0.16 1.00 1.00 0.180 0.03 1.00 1.00 1.00 1.00 0.17 0.70 5 0.00 0.48 1.00 1.00 0.167 1.00 1.00 1.00 0.68 1.00 0.21 0.82 6 0.00 0.00 1.00 1.00 0.250 0.00 1.00 1.00 1.00 1.00 0.11 0.68 7 0.31 0.36 0.81 1.00 0.117 0.84 0.91 0.21 1.00 1.00 1.00 0.83 8 0.31 0.36 0.81 1.00 0.117 0.84 0.91 0.21 1.00 1.00 1.00 0.83 9 0.00 0.81 0.91 1.00 0.146 1.00 1.00 1.00 1.00 0.64 1.00 0.94 10 0.00 0.46 1.00 1.00 0.170 1.00 1.00 1.00 1.00 1.00 0.20 0.87 11 0.39 0.47 0.81 1.00 0.092 0.22 0.56 0.41 1.00 0.15 1.00 0.56 12 0.39 0.47 0.81 1.00 0.092 0.33 0.57 0.41 1.00 0.15 1.00 0.58 13 0.00 0.68 0.86 1.00 0.151 1.00 1.00 1.00 1.00 0.30 1.00 0.88 14 0.00 0.00 0.98 1.00 0.250 0.00 1.00 1.00 1.00 1.00 0.24 0.71 15 0.28 0.40 0.86 1.00 0.120 0.91 1.00 0.30 1.00 1.00 1.00 0.87 16 0.27 0.30 0.86 1.00 0.134 0.83 1.00 0.27 1.00 1.00 1.00 0.85 In addition, we visualize the assignment of work packages to internal resources over time on the company level as in Figure 3(a). The work done by all internal resources ( k∈R i x kst ) is shown on the ordinate while the time line is on the abscissa. The four skills are represented by four different curves. For example in Figure 3( skill 1 is not done internally at all while the amount of work that is done internally for skill 2 increases with each period. Note that the curves for skill 3 and 4 are super-imposed at a maximum of k∈R i x k3t = k∈R i x k4t = 45.

Outsourcing decisions

As stated in Section 3, due to the chosen problem parameters, the average utilization of the internal resources is higher than 100%. Therefore, it is not possible to do all the work internally. From Table 5 we can observe that, in general, those work packages are outsourced (low values of ρ 1 and ρ 2 ) which require cheap skills, i.e. a skill s that can be outsourced at relatively low costs c e s . On the other hand work packages which require expensive skills are (almost) completely done internally (ρ 3 and ρ 4 close or equal to 1). As expected, the inequality ρ s ≥ ρ s holds if c e s > c e s .

Impact of the shape of the learning curve

Impact on the company When comparing instances 1 and 5 with its counterparts 9 and 13 that have a steeper learning curve, we can see that a steeper learning curve leads to a broader qualification of the company, since CV(X C ) in Table 6 is lower for instances 9 and 13. This holds regardless of the depreciation of knowledge.

# λ ks β kst ρ 1 ρ 2 ρ 3 ρ 4 CV(X C ) 1 0.012 0 0.00 0.62 1.00 1.00 0.156 5 0.012 10 0.00 0.48 1.00 1.00 0.167 9 0.02 0 0.00 0.81 0.91 1.00 0.146 13 0.02 10 0.00 0.68 0.86 1.00 0.151 Table 6: Fraction ρ s of work done internally and specialization of the company CV(X C ) for instances 1, 5, 9 and 13 Figure 3(a) which graphically shows the solution for instance 1 illustrates this finding. Skill 3 and 4 are completely done by internal resources ( k∈R i x kst = r st = 45), skill 1 is completely outsourced ( k∈R i x kst = 0) and the amount of work that is done internally for skill 2 increases over time. The increase of internal work for skill 2 can be explained with the four diagrams shown in Figure 3(b). It stems from the efficiency gains of resources 3 and 6, that enable resource 5 to use more of its capacity performing skill 2. Furthermore, resources 5 and 1 are becoming more efficient in doing skill 2 as well. At first sight it seems that the following greedy heuristic would provide optimal results: i) Order work packages in decreasing order of external skill costs. ii) Beginning with the work package requiring the most expensive external skill, assign work packages to available internal resources as long as possible. iii) Outsource the remaining work. However, Figure 4(a) gives a counterexample. The figure shows the solution for instance 9 which, compared to instance 1, has a steeper learning curve. The difference compared to instance 1 is that some of the work for the more expensive skill 3 has been outsourced ( k∈R i x kst < r st = 45) in favour of skill 2 which is now done more intensively by internal resources and less outsourced. The reason for this phenomenon can be examined using Figure 4(b) which reveals that the full dedication of resources 2 and 4 to skill 3 and resource's 5 higher dedication to skill 2 leads to lower overall costs.

Impact on individual resources

The impact of the learning curve's shape on the individual resources can be examined with a comparison of σ k for instances 1 and 5 with its counterparts 9 and 13 in Table 7. Regarding these four instances a steeper learning curve seems to lead to more specialized individual resources since σ k is higher on average. These numbers are supported by the fact that there are 4 fully specialized resources (σ k = 1) in instance 1 while there are 5 of them in instance 9. Resource 5 (with skills 2 and 4) and resource 6 (with skills 3 and 4) are both capable of performing skill 4, while skill 2 (the other skill of resource 5) is cheaper than skill 3 (the other skill of resource 6) (cf. Table 1). Therefore, one would expect that resource 5 (with the cheaper other skill) is devoted to perform skill 4 rather than resource 6 (with the more expensive other skill). However, in instance 9 resource 5 is qualified for both skills while resource 6 is fully specialized for skill 4. This contradicts the intuition that a higher difference of skill values (as for resource 5) should lead to a stricter specialization for the more expensive skill.

# λ ks β kst σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 avg. 1 0.012 0 1.00 1.00 1.00 1.00 0.37 0.48 0.81 5 0.012 10 1.00 1.00 1.00 0.68 1.00 0.21 0.82 9 0.02 0 1.00 1.00 1.00 1.00 0.64 1.00 0.94 13 0.02 10 1.00 1.00 1.00 1.00 0.30 1.00 0.88 Table 7: Specialization σ k of internal resources for instances 1, 5, 9 and 13

The numbers disclose that the specialization measures on the individual and company level diverge. A steeper learning curve leads to more specialized individual resources in terms of σ k but a broader qualified company in terms of CV(X C ).

Impact of company skill level targets

Company skill level targets broaden the company's qualification. This can be seen by comparing the coefficients of variation CV(X C ) of instances 1 and 3 (0.156 5, respectively. The proposition is also supported by comparison of Figures 3(a) and 5(a). We can clearly see that in the latter figure the cheapest skill 1 is used and trained to reach the required company skill level targets of φ s = 2.8 for all skills, while it is completely outsourced in Figure 3(a). In order to reach this goal less work of skill 2 is done by internal resources which in turn leads to a higher amount of outsourced work demanding this skill. Of course, increasing company skill level targets φ s lead to higher overall costs since the solution space is reduced due to constraints [START_REF] Bassett | Assigning projects to optimize the utilization of employees' time and expertise[END_REF]. Due to the different shapes of the learning curves the company skill level targets φ s are hardly comparable in a quantitative manner. Therefore, Figure 6 depicts the dependence of the costs on the company skill level targets for further variations of instance 1 with a flat learning curve. For φ s ≥ 3 no feasible solutions could be found.

Impact of knowledge depreciation

We observed basically similar results comparing instances 1 and 5 as well as 9 and 13. No structural difference on the company level between the solutions of different levels of knowledge depreciation can be identified for instances with company skill level targets φ s = 0: As can be seen in Figure 3(a) and Figure 7 3). No particular impact of the influence of knowledge depreciation on specialization measures CV(X C ) and σ k can be observed.

We illustrate the combined effects of knowledge depreciation and synchronously imposed company skill level targets φ s = 2.5 by Figure 8(a) which shows the results for instance 7. As before, assignment in decreasing order of external skill costs is basically still reasonable. But as it has been observed for the case of company skill level targets without knowledge depreciation, the amount of work for skill 1 done by internal resources must be non-zero.

It is noteworthy that, in order to reach a company skill level target in constraints (5), only the cumulated amount of work done in that skill is important. This is due to constraints (3) in which depreciation of knowledge is a subtractive term. Thus, it is a constant at the end of planning horizon regardless of the assignment decisions. Therefore, scheduling of work package assignments is irrelevant for reaching company skill level targets φ s . In fact, scheduling of the work package assignments is only important for the objective function [START_REF] Alba | Software project management with GAs[END_REF]. In case of knowledge depreciation and company skill level targets the solutions do not only depend on the assignment decisions (i.e. who does what) alone but on the time line of assignments (i.e. scheduling, who does what in which period) as well. The decision of doing a skill rather in the beginning or in the end of the planning horizon in order to reach company skill level targets at minimum costs is obviously important. If an amount r s of work has to be done in one single period t by resource k subject to knowledge depreciation, it can be done in less time in period t = 1 (when f ks (•) = f ks (z ks0 )) than in period t = T (when f ks (•) f ks (z ks0 ) due to knowledge depreciation). In this context, it makes sense that in Figure 8 With costs of c i k = 500 the allocation of internal resources to cheap skills (c e 1 = 400) is avoided, even if internal resources are not completely utilized. For the parameters of the relevant test instances internal resources cannot compete with external resources in terms of costs for skill 1 and 2 despite learning effects. This can be seen in the results when comparing ρ s of odd instance numbers (c i k = 0) with the values of even instance numbers (c i k = 500). However, in order to reach company skill level targets φ s > 0 in instances 4, 8, 12 and 16 the assignment of relatively expensive internal resources for these skills is mandatory.

Conclusions and Outlook

In this paper we presented a constrained non-linear continuous optimization model to address the problem of cost-minimal assignment of project work to internal and external multi-skilled human resources. The effects of learning and depreciation of knowledge as well as the goal of reaching company skill level targets at the end of the planning horizon are taken into account.

Our results show that our heuristic implementation has a promising computational performance. The use of random starting points improved the quality of the solutions significantly for some instances.

We analysed the solutions w.r.t. the influence of problem parameters. For selected parameter combinations the intuition is supported that externally expensive skills should be done by internal resources while cheaper skills should be outsourced. However, several observations showed that an implementation of this policy as a greedy assignment heuristic will not provide optimal results. Furthermore, we found out that the shape of the learning curve leads to different qualification strategies which depend on the level of analysis: Faster learning makes human resources more specialized which paradoxically leads to a broader qualification of the company. Company skill level targets are an effective tool to broaden the qualification of the company but they can induce high costs. Due to the objective of cost-minimization, depreciation of knowledge together with company skill level targets cause assignment decisions to depend on the time line (i.e. scheduling): Skills which are less expensive when outsourced are used more intensively by internal resources in the beginning than in the end of the planning horizon in order to minimize costs.

Further research could compare our random starting point approach with a more sophisticated selection of starting points which ensures diversification. Furthermore, other global search techniques like e.g. genetic algorithms might also be well-suited to solve the problem. The fractional assignment decisions could be coded as a string and their fitness can then be evaluated easily (cf. Yoshimura et al. [START_REF] Yoshimura | Decision-making support system for human resource allocation in product development projects[END_REF]). Finally, the problem could be extended to include project scheduling decisions in order to facilitate the possibility of demand levelling.

Figure 1 :

 1 Figure 1: Production possibility frontier

b

  |R e s | = 1 ∀s and R kt = ∞ ∀t, k ∈ R e s . External resources will not be subject to learning nor forgetting and their unit production time is 1, i.e. a ks = 0 and ks = 1 ∀s, k ∈ R e

Figure 2 :# λ ks β kst φ s c i k 1

 21 Figure 2: Plots of the flat (λ = 0.012) and the steep (λ = 0.02) learning curve f (z)

x C 3 ,

 3 . . .} be sets that combine the data points of the different skills. The coefficient of variationCV(X C ) = VAR (X C ) E (X C )

Figure 3 (

 3 b) depicts the assignments on an individual level. Each subdiagram shows the assignments of work requiring one skill to the resources. The work done by internal resources is shown as stacked bars on the ordinate while the time line is on the abscissa. For example in Figure 3(b) the top right subdiagram depicts skill 2. The work requiring that skill is done by resource 1 and furthermore by resource 5 in a growing amount. The height of the stacked bars correspond to the strictly monotone increasing curve for skill 2 in Figure 3(a).

Figure 3 :

 3 Figure 3: Solution of instance 1 with λ ks = 0.012, φ s = 0, β kst = 0, and c i k = 0

13

 13

Figure 4 :

 4 Figure 4: Solution of instance 9 with λ ks = 0.02, φ s = 0, β kst = 0, and c i k = 0

  138), 5 and 7 (0.167 vs. 0.117), 9 and 11 (0.146 vs. 0.092), as well as 13 and 15 (0.151 vs. 0.120) in Table

Figure 5 :

 5 Figure 5: Solution of instance 3 with λ ks = 0.012, φ s = 2.5, β kst = 0, and c i k = 0

Figure 6 :

 6 Figure 6: Influence of company skill level targets φ s on costs

  (b) resource 3 is completely dedicated to skill 1 during the entire planning horizon while resources 1 and 2 perform skill 1 only in the first periods to reach the company skill level targets.

Figure 7 :

 7 Figure 7: Solution of instance 5 with λ ks = 0.012, φ s = 0, β kst = 10, and c i k = 0

Figure 8 : 0 Page

 80 Figure 8: Solution of instance 7 with λ ks = 0.012, φ s = 2.5, β kst = 10, and c i k = 0

Table 1 :

 1 Parameters and skill matrix of the base case

				Skill
			1	2	3	4
	c e s	400 500 600 700
		1 x	x
	Resource	2 x 3 x 4 5	x x	x x	x x
	s	6			x	x

Table 2 :

 2 Variations of the base case

	8

Table 3 :

 3 Computational resultsIn order to analyze the limits of the Ipopt-algorithm we additionally employed larger instances as given in Table4. The instances are based on instance #1 given in Table2, where the number of resources |R i | and skills |S| has been increased and two random skills were assigned to each resource. For each combination of |R i | and |S| we generated 5 test instances. The results are given in Table

		sec]	best avg dev % max dev %
	1	13.5 159,164	0.00	0.00
	2	21.8 516,994	0.06	0.18
	3	18.1 163,787	0.00	0.00
	4	24.9 523,407	0.00	0.02
	5	17.9 177,638	0.06	1.38
	6	32.7 531,690	0.00	0.00
	7	18.7 191,577	0.35	2.07
	8	23.6 551,577	0.11	0.55
	9	14.8 147,619	0.13	0.14
	10	18.1 506,487	0.18	0.26
	11	14.0 168,480	0.00	0.00
	12	17.4 528,480	0.00	0.00
	13	18.5 173,990	0.60	8.30
	14	30.6 527,117	0.72	1.85
	15	20.3 182,632	1.45	3.52
	16	27.7 542,631	0.23	0.91
	mean	20.8 349,579	0.24	1.20

  To gain managerial insights the results are analysed on the individual level of a resource and on the aggregate level of the company, in which the values of all internal resources are summed up. Let

	|S| avg. time [sec] max dev %
	25 10	45.5	0.00
	50 25	110.7	0.00
	100 50	494.5	0.00
	Table 4: Additional computational results
	4.2 Managerial insights
	4.2.1 Analysis methods
	x R ks =	x kst
		t
	be the amount of work done by resource k for skill s during the entire planning
	horizon and let
	x C s =
		k∈R i

Table 5 :

 5 Fraction ρ s of work done internally, specialization of the company CV(X C ) and specialization σ k of internal resources

  (a), expensive skills are done completely internally while cheap skills are outsourced. On the individual level, however, the assignments of work to resources shift. Of course, the fraction

		195
		190
	Total costs in thousands	165 170 175 180 185
		160
	155 F r o 2.6 2.65	2.7	2.75	2.8	2.85	2.9	2.95	3
		P
		e
		e r
				R
				e
				v i e
						w
							O n l
								y

Page 8 of 24 http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk International Journal of Production Research

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

Page 13 of 24 http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk International Journal of Production Research

Page 18 of 24 http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk International Journal of Production Research


M Bargès

S Loisel

X Venel

On finite-time ruin probabilities with reinsurance cycles influenced by large claims

Keywords: JEL classification codes: G22 C60 Finite-time ruin probability, reinsurance cycles, Erlangization, dependence in risk theory, phase-type distributions

Market cycles play a great role in reinsurance. Cycle transitions are not independent from the claim arrival process : a large claim or a high number of claims may accelerate cycle transitions. To take this into account, a semi-Markovian risk model is proposed and analyzed. A refined Erlangization method is developed to compute the finite-time ruin probability of a reinsurance company. Numerical applications and comparisons to results obtained from simulation methods are given. The impact of dependency between claim amounts and phase changes is studied.

Introduction

The classical Sparre-Anderson risk model describes the surplus process {R t } t 0 for an insurance or a reinsurance company as R 0 = u and

R t = u + ct -S t ,
where u represents the non-negative initial capital and c > 0 is the premium income rate. The compound renewal process {S t } t 0 represents the cumulated claim amount up to time t and is 1 expressed as

S t = 0, if N t = 0, Nt i=1 X i otherwise,
where the claim amounts X i , i 1, are non-negative and i.i.d. random variables. The number of claims up to time t is modeled by the renewal process {N t } t 0 with i.i.d. inter-arrival times denoted by W k , k 1. The X i 's, i 1 are independent from the process (N t ) t 0 . One quantity of interest (among others) to set risk tolerance and risk limits is the finite-time ruin probability ψ(u, T ) = P [∃t ∈ [0, T ], u + ct -S t < 0] , u 0, T 0.

Different methods to compute this finite-time (non-)ruin have been proposed in the literature. [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] derived an explicit formula for the finite-time non-ruin probability in a compound Poisson model where the claim amounts are integer-valued. This idea was then widely used and extended to more general models as in [START_REF] De | Explicit finite-time and infinite-time ruin probabilities in the continuous case[END_REF], [START_REF] Ignatov | A finite-time ruin probability formula for continuous claim severities[END_REF], [START_REF] Rullière | Another look at the Picard-Lefèvre formula for finite-time ruin probabilities[END_REF] or [START_REF] Lefèvre | Finite-time ruin probabilities for discrete, possibly dependent, claim severities[END_REF], among others. Asymptotic methods have also been used to evaluate the finite-time ruin probability as in [START_REF] Biard | Impact of correlation crises in risk theory: asymptotics of finite-time ruin probabilities for heavy-tailed claim amounts when some independence and stationarity assumptions are relaxed[END_REF] where the claim sizes are heavy-tailed and where stationarity and independence assumptions are relaxed in possible correlation crises. Previous general asymptotic results have been provided for sums of random variables as in [START_REF] Wüthrich | Asymptotic value-at-risk estimates for sums of dependent random variables[END_REF], Albrecher et al. (2006) or [START_REF] Kortschak | Asymptotic results for the sum of dependent nonidentically distributed random variables[END_REF].

The method we focus on is called Erlangization. It was first introduced in Avram and Usabel (2003) but principally developed by [START_REF] Asmussen | Erlangian approximations for finite-horizon ruin probabilities[END_REF]. The Erlangization method consists in replacing the deterministic time horizon T by an Erlang-distributed random variable H with mean T and independent from the claim process. [START_REF] Asmussen | Erlangian approximations for finite-horizon ruin probabilities[END_REF] have shown that the ruin probability before the independent Erlang time H coincides with the ruin probability in a certain Markovian fluid model with a matrix-exponential form. They considered the compound Poisson risk model with phase-type claim amounts. Stanford et al. (2005a) proposed an extension to the general Sparre-Anderson risk model. In [START_REF] Yu | Erlangian approximation to finite time ruin probabilities in perturbed risk models[END_REF], the risk process with phase-type distributed interclaim times and claim amounts is perturbed by a standard Brownian motion. [START_REF] Ramaswami | The Erlangization method for Markovian fluid flows[END_REF] also extended the Erlangization method to Markovian modulated fluid flow models which will prove to be essential our work. Erlangization has also been used in some other domains as in [START_REF] Stanford | Erlangized fluid queues with application to uncontrolled fire perimeter[END_REF] where the forestry problem of the evolution of uncontrolled fire perimeter over time is treated. For a complete review on fluid flow matrix analytic methods applied in ruin theory, see [START_REF] Badescu | Applications of fluid flow matrix analytic methods in ruin theory, a review[END_REF].

We propose here to consider the reinsurance market. In the reinsurance business, the presence of underwriting cycles cannot be ignored. As said in [START_REF] Cummins | The costs and benefits of reinsurance[END_REF], reinsurance cycles are composed of soft phases, during which reinsurance supply is extensive and reinsurance premium income rates are quite low and decreasing, and hard phases where the supply is too low and the premium income rate becomes higher and higher. Phase transitions are hard to explain, and seem to happen due to both endogenous and exogenous events. In [START_REF] Weiss | Underwriting cycles: a synthesis and further directions[END_REF], it is mentioned that these changes can come from institutional features as well as the effects of real phenomena. These phenomena can be asymmetric information in capital markets, capital surpluses and shortages in insurance or interest and/or loss shocks.

A simple generalization of the above-mentioned Erlangization model could be obtained with a Markovian environment process which would modulate the premium income rate, following the idea of Asmussen (1989) and of many others. Nevertheless, this first generalization leads to a model in which the cumulated claim amount process has no feedback control on the phases of the Markovian environment process {J t }. This is not the case in real world, as changes in cycle phases may be triggered by a large claim. In this paper, we propose a generalization of the above risk model that takes into account this possible trigger.

We consider a simplification of a multi-phases cycle as it is exposed in Subsection 3.1 into a 2-phases reinsurance cycle consisting of a soft phase with low premium rate and a hard phase with high premium rate. In financial insurance, it is usual to suppose only two environment states as in the model introduced by [START_REF] Hardy | A regime-switching model of long-term stock returns[END_REF] where the volatility of a stock is either high or low according to the state environment. Considering two environment phases permits of course not only model simplifications but also avoids many estimation problems which appear when more phases are assumed.

In the Solvency II framework, the capital required to cover risk is determined in such a way that a 1-year Value-at-Risk is controlled. Nevertheless, for reinsurance risk, and in particular for catastrophe risk, 1-year occurrence probabilities of extreme events are small, and risk must be controlled within a larger time period, 5 to 10 years, to determine not only solvency requirements but also long term economic capital in internal risk models. The average length of reinsurance cycles typically ranges from 4 to 10 years, see [START_REF] Venezian | Ratemaking methods and profit cycles in property and liability insurance[END_REF], [START_REF] Cummins | An international analysis of underwriting cycles in propertyliability insurance[END_REF] and [START_REF] Weiss | Underwriting cycles: a synthesis and further directions[END_REF]. It is thus necessary to use a risk model that incorporates this cyclic behavior.

Even if the numerical application provided in Section 6 is detailed when only the premium income rates changes according to the cycle phases, our model is developed for the inter-claim time and claim amount distributions to be dependent from the environment cycle process. The feedback of the claim amounts on the environment process thus involves some dependence between inter-claim times and claim amounts. Other types of dependence have already been explored in the study of infinite-or finite-time ruin probability. [START_REF] Boudreault | On a risk model with dependence between interclaim arrivals and claim sizes[END_REF] and [START_REF] Albrecher | A ruin model with dependence between claim sizes and claim intervals[END_REF] assumed that each claim amount depends on the previous inter-claim time or inversely that each inter-claim time depend on the previous claim amount. This has been extended to claim amounts depending on several preceding inter-claim times as can be seen in [START_REF] Biard | Asymptotic finite-time ruin probabilities for a class of path-dependent claim amounts using Poisson spacings[END_REF]. In [START_REF] Albrecher | Exponential behavior in the presence of dependence in risk theory[END_REF] the dependence between the inter-claim times and claim amounts is introduced by an arbitrary copula.

One of the contributions of the paper is to propose a more suitable model for long term economic capital management for reinsurance companies. Within this framework, the dependence of the premium income rate, the claim amount distribution and the inter-claim time distribution with the cycle environment is taken into account. The model also permits a feedback control from the claim amounts on the cycle phases. We adapt the Erlangization method to this context in order to approximate the finite-time ruin probability. We use intensive calculation techniques to this aim and compare our results to Monte-Carlo simulations. We evaluate the model error in several examples and quantify the impact of the claim amounts feedback over a 5-years time horizon. We show that for some reasonable values, the finite-time ruin probability observes a decrease of 6% when this correlation is not considered. This is of course just an illustrative example, and not a general fact. This paper begins with a recall on the Erlangization method and related results. In Section 3 reinsurance cycles are discussed. General notations are given and a first naive model is proposed where we show how the Erlangization method can be adapted to deterministic phase changes. The main model incorporating the impact of the claim amounts on the change in cycle phases is presented in Section 4. Section 5 is dedicated to the control of the precisions of the proposed main model. Finally, an extended and practical illustration is provided in Section 6 where fitting methods and finite-time ruin probability approximations are given.

Erlangization method

We recall in this section the principle of the Erlangization method and the central theorems from [START_REF] Asmussen | Erlangian approximations for finite-horizon ruin probabilities[END_REF] which lead to an approximation for the finite-time ruin probability. We follow the presentation of [START_REF] Asmussen | Erlangian approximations for finite-horizon ruin probabilities[END_REF] and [START_REF] Ramaswami | The Erlangization method for Markovian fluid flows[END_REF] and slightly adapt the notation for it to be consistent with our purpose. First, the surplus process R t is embedded in an initial continuous Markovian additive fluid process {( Jt , Ṽt )} where Jt is a fluid flow and Ṽt a Markovian environment process. The state space Ẽ of Ṽt consists of two subspaces Ẽand Ẽ+ such that when Jt is in state i ∈ Ẽ-, Ṽt is in a decreasing inter-claim time state (no jump) with fluid rate r i = -c i and when Jt is in state i ∈ Ẽ+ , Ṽt is in an increasing claim state (jump) with fluid rate r i = 1. The transition intensity matrix of Jt is noted as

Λ = Λ --Λ -+ Λ +-Λ ++ ,
where Λ i 1 i 2 correspond to the transition from a state of subspace Ẽi 1 to a state of subspace Ẽi 2 with i 1 and i 2 ∈ {-, +}.

The aim of the Erlangization method is to evaluate the finite-time ruin probability before an Erlang horizon H ∼ Erlang(a, L). The finite-time ruin probability can thus be seen as the probability that the fluid flow reaches the initial reserve u in an increasing state (claim state) conditioned on starting in a decreasing state (inter-claim time state) before the expiry of the Erlang horizon H. For this purpose, we construct an expanded absorbing Markovian fluid process {(J t , V t )}. The state space E of the Markovian environment process J t is now composed of three subsets E 0 , E -and E + where E 0 is the one state absorbing space, E -is the decreasing inter-claim time state space and E + is the increasing claim state space.

The transition intensity matrix for the Markovian environment process

J t is Θ =   0 0 0 h -H ⊕ Λ --I L ⊗ Λ -+ 0 I L ⊗ Λ +-I L ⊗ Λ ++   =   0 0 0 h -Θ --Θ -+ 0 Θ +-Θ ++   , ( 1
)
where H is the phase generator of the Erlang horizon H and h -is the vector of absorbtion probabilities for the fluid process in a decreasing state. Operators ⊗ and ⊕ respectively denote the Kronecker product and the Kronecker sum.

The incorporation of the premium income rates which correspond to the fluid rates for the fluid process V t is done by time scaling. These fluid rates are r i = 0 for i ∈ E 0 , r i = -c i for i ∈ E -and r i = 1 for i ∈ E + . The matrix Θ which takes the premium income rates into account associated with Θ is

Θ =   0 0 0 (I L ⊗ R -) -1 × h -(I L ⊗ R -) -1 × Θ --(I L ⊗ R -) -1 × Θ -+ 0 Θ +- Θ ++   =   0 0 0 (C -) -1 × h -(C -) -1 × Θ --(C -) -1 × Θ -+ 0 Θ +- Θ ++   =   0 0 0 h -Θ --Θ -+ 0 Θ +-Θ ++   ,
where R -is a diagonal matrix formed with the decreasing fluid rates.

Denote by η ij the probability that when J 0 = i ∈ E -, then the first upcrossing of level 0 of {V t } occurs at time t with J t = j ∈ E + . The matrix η corresponds to the matrix of the "upcrossing phase probabilities" at the completion of a downwards excursion for {V t }, conditioned on starting in a decreasing state. We now cite some results from [START_REF] Asmussen | Erlangian approximations for finite-horizon ruin probabilities[END_REF], which can also be seen in [START_REF] Ramaswami | The Erlangization method for Markovian fluid flows[END_REF] which aim to give an evaluation for the probability ψ(u, H) of reaching u conditioned on starting in a decreasing state before the expiry of the Erlang horizon. We first reccal a general result on the upcrossing probabilities of fluid models.

Theorem 1 Let {(J t , V t )} be a Markovian fluid process and Θ the transition intensity matrix of the Markovian environment process {J t }. Let also η denote the "upcrossing phase probabilities" at the completion of a downwards excursion for {V t }, conditioned on starting in a decreasing state i ∈ E -. Note by ψ(u) the vector of the infinite horizon phase-distribution at u (the initial surplus, starting at 0),conditioned on starting in a decreasing state. Then we have : (a) the matrix η satisfies the Riccati equation:

ηΘ ++ + ηΘ +-η + Θ --η + Θ -+ = 0 (2) and ψ(u) = ηe U u 1 where U = Θ ++ + Θ +-η.
(b) let now δ > 0 be a "convergence improving parameter", satisfying δ > max(-(Θ --) ii ) for all state i of H. Then, the solution η of the previous Riccati equation can be computed by the iteration scheme η = lim n→∞ η (n) where η (0) = 0,

η (n+1) = ((δI + Θ --)η (n) + Θ -+ )(δI -(Θ ++ + Θ +-η (n) )) -1 . (3)
Focusing on our case where we look at the fluid process before the expiry of the Erlang horizon time H with initial probability vector γ, we have the following corollary.

Corollary 2 Including the Erlang time H with the initial probability vector γ, it follows that

ψ(u, H) = γηe V u 1, ( 4
)
where

V = Θ ++ + Θ +-η.
Finally, [START_REF] Asmussen | Erlangian approximations for finite-horizon ruin probabilities[END_REF] obtained the following convergence result.

Theorem 3 For T > 0, H the Erlang distribution with L stages and mean T , then

ψ(u, H) → ψ(u, T ) as L → ∞, at rate O(L -1 ).
The proof of Theorem 1 is based on the fluid flow theory and on the study of the steadystate distribution of the general fluid process {(V t , J t )}. Indeed, [START_REF] Asmussen | Stationary distributions for fluid flow models with or without Brownian noise[END_REF] proved that the distribution of the steady-state variable (V, J) is described by V being phase-type given J = i, with phase generator U independent from i, the unit column vector e as initial vector for i ∈ E + and η i as the ith row of η, for i ∈ E -. The Corollary 2 is then deduced when including the Erlang time horizon. The details of the proof of Theorem 3 can be found in [START_REF] Asmussen | Erlangian approximations for finite-horizon ruin probabilities[END_REF].

We show in appendix how to deal with the matrix inversions that are involved in Equation (3).

A few words on reinsurance cycles modeling

Reinsurance cycles

In order to give an idea of a possible generalization of the model we propose in this paper, we first expose here a cycle modeling with several phases. Then we recall the reasons of our simplification choice to a 2-phases cycle.

Reinsurance cycles can be composed of several phases with soft phases where the premium rates are low and hard phases with high premium rates. Intermediate phases have premium rates varying according to their proximity with the soft or the hard phases. We suppose that a reinsurance cycle has, for example, 10 phases where i = 1 corresponds to the softest phase and i = 10 to the hardest phase. We note c i , i = 1, ..., 10, the premium income rate for the 10 phases with c 1 c 2 ... c 10 and assume that when the cycle is in phase i, it could only go to phase i -1, i + 1, i + 2 or stay at i depending on the claim amount. When a large claim occurs, the cycle have high probability to go to phase i + 1 and if the claim is very large, the cycle can also directly jump to phase i + 2. Figure 1 shows a possible cycle phase path using these ideas. We could also double the phase states such as 1 ↓ , 1 ↑ , 2 ↓ , 2 ↑ , 3 ↓ , 3 ↑ , ... (or 1, 2 ↓ , 2 ↑ , 3 ↓ , 3 ↑ , ..., 10) in order to keep the general trend of the cycle. For example, we could set probabilities at state 3 according to the cycle trend as in Table 1. But, in this paper, we focus our study on a reinsurance cycle that is composed of only a soft phase with low premium rate and a hard phase with high premium rate. Indeed, dealing with models with more than two phases in a cycle would require very large transition matrices and would highly increase computational issues. Moreover, it is conventional in finance insurance to consider environment processes with 2 states, as in [START_REF] Hardy | A regime-switching model of long-term stock returns[END_REF], which furthermore avoids many estimation problems.

We assume that large claims have high probabilities to cause a change from the soft phase to the hard phase. This impact of the claim amounts on the environment process represents quite well the general reinsurance loss market as most of insurance and reinsurance companies hold high expositions in the most risky regions (Western Europe, Japan, Eastern United-States, South-East Asia, China,...). In consequence, a big event (catastrophe,...) impacts all the actors with a very large claim.

General notation

We fix now some notation for all the following sections. We suppose that in each cycle phase the claim amounts are i.i.d. and phase-type distributed such as X i ∼ P H(G i , α i ) where G i is the phase generator matrix and α i the initial row vector in cycle phase i. The inter-claim times are also i.i.d. and phase-type distributed. The inter-claim time random variable for cycle phase i is noted as W i ∼ P H(K i , β i ) with K i the phase generator matrix and β i the initial row vector. We write by g i = -G i e and k i = -K i e the rates of absorption for respectively X i and W i where e is the unit column vector. Denote by T the fixed time horizon. The phase-type parametrization of the Erlang horizon time H with L stages, mean T and rate a = L/T is H ∼ P H(H, γ = (1, 0, ..., 0)) where

H =      -a a 0 0 0 . . . . . . 0 0 0 -a a 0 0 0 -a      .

Deterministic phase changes

We give first a naive model which assumes that only one change of phase happens and that this time of phase change is known in advance. We present here the case where the cycle environment process begins in the soft phase and switches to the hard phase at a deterministic time T 1 and stops at the horizon T . But this model could be adapted with the cycle beginning in the hard phase and going to the soft phase. The times spent in the soft phase and in the hard phase are respectively T 1 and T 2 = T -T 1 . Here, the Erlang horizon time H with mean T and L stages is the sum of two intermediate Erlang horizon times H 1 and H 2 with means respectively T 1 and T 2 and L 1 and L 2 stages. We can then write the transition intensity matrix Θ as

Θ =   0 0 0 h -Θ --Θ -+ 0 Θ +-Θ ++   =       0 0 0 0 0 0 H 1 ⊕ K 1 M 12 Id L 1 ⊗ k 1 α 1 0 h 2 ⊗ e k 2 0 H 2 ⊕ K 2 0 Id L 2 ⊗ k 2 α 2 0 Id L 1 ⊗ g 1 β 1 0 Id L 1 ⊗ G 1 0 0 0 Id L 2 ⊗ g 2 β 2 0 Id L 2 ⊗ G 2      
, where M 12 = ((h 1 × γ 2 ) × e k 1 ) × β 2 is the bloc matrix of transitions between the two phases and h i = -H i e the absorption rate vector for the phase-type parametrization of the Erlang horizon H i .

Remark 4 The way to determine the numbers of stages L 1 and L 2 for respectively H 1 and H 2 is shown in appendix.

Example 5 Consider that H 1 ∼ Erlang(1, a) and H 2 ∼ Erlang(2, a), then we have

Θ =           0 0 0 0 0 0 0 0 K 1 -aId n1 (a × e n1 ) × β 2 0 k 1 α 1 0 0 0 0 K 2 -aId n2 aId n2 0 k 2 α 2 0 ae n2 0 0 K 2 -aId n2 0 0 k 2 α 2 0 g 1 β 1 0 0 G 1 0 0 0 0 g 2 β 2 0 0 G 2 0 0 0 0 g 2 β 2 0 0 G 2           ,
where n 1 is the number of stages of W 1 and n 2 the number of stages of W 2 .

Using the same idea, it is also possible to construct an intensity matrix Θ for the situation where the environment process begins in the hard phase and goes to the soft phase at a deterministic time. But in fact, cycle phase changes are random and may be influenced by large claims.

Phase change instants dependent from claim amounts

We now focus on the main model of the article where we let the cycle switch from the soft phase to the hard phase when a claim amount is greater than a threshold x 0 0. In what follows, we suppose that the cycle begins in the soft phase but using the idea of the previous simple model, one can add an initial hard phase with a deterministic length. As noted in introduction, the going back to the soft phase should also depend from the claim amount, for example according to the severity of the claim that caused the passing to the hard phase. But for simplification reasons, we suppose here that it happens after a deterministic time t 0 0. The fact that the time spent in the previous soft phase is not known, we can not proceed as in the deterministic phase change model developed above.

We require here that the claim amounts are distributed as mixtures of Erlang distributions. Mixtures of Erlang distributions include the exponential distribution, the Erlang distribution and mixtures of Exponential distributions. The canonical random variable for the claim amounts in cycle phase i, i = 1, 2, is then written as

X i = n i j=1 α ij Y ij ,
where Y ij are n i independent and Erlang distributed random variables such that Y ij ∼ Erlang(m ij , λ ij ) with shape parameters m ij > 0 and rate parameters λ ij > 0. The scalars α ij 0 are the mixture proportions and must satisfy α i1 +...

+α in i = 1. We assume that E[Y i1 ] = m i1 /λ i1 > ... > E[Y in i ] = m in i /λ in i .
This section is decomposed into 2 subsections. The first one introduces the transition from the soft phase to the hard phase which depends on the claim amounts. The second subsection exposes the way the go back to the soft phase.

From the soft phase to the hard phase

We first introduce the transition from the soft phase to the hard phase which is dependent from the claim amounts. To modelize the fact that the environment process goes to the hard phase when a claim amount is greater than the threshold x 0 , we use the probabilities

p 1j = P (Y 1j > x 0 ) = 1 -F Y 1j (x 0 ) (5)
where Y 1j ∼ Erlang(m 1j , λ 1j ), for j = 1, ..., n 1 . As the random variables Y 1j are ordered by decreasing mean, we have p 11 > ... > p 1n 1 . Then we integrate these probabilities into the transition matrix Λ of the initial Markov process {I t } as follows. Construct line vectors p 1j = (0, ..., 0, p 1j ) of length m 1j and concatenate them in p 1 = ( p 11 , ..., p 1n 1 ). Then, set

P 1 = (p 1 ⊗ e k 1 ),
the matrix containing the probabilities for the claim amounts to be greater than x 0 when the environment process is in the soft phase as needed to incorporate them into Λ. The latter becomes

Λ =     K 1 0 k 1 α 1 0 0 K 2 0 k 2 α 2 P - 1 g 1 β 1 P 1 g 1 β 2 G 1 0 0 g 2 β 2 0 G 2     ,
where

P - 1 = ((e n 1 -p 1 ) ⊗ e k 1 ),
and the operator is the Hadamard product (entrywise product). Thus the transition intensity matrix Θ of {J t } can be obtained as exposed in Section 2. The concept of dependent changes is here probabilistically respected to approximate the ruin probability.

From the hard phase to the soft phase

We now include the fact that the cycle process can go back to the soft phase from the hard phase. As said previously, we do not suppose that this transition is dependent from the claim amount but that the expected time spent in the hard phase is fixed at time t 0 . This transition can happen at any state of the Erlang time H when the process is in the hard phase. The possibility to go from the hard phase to the soft phase is introduced by a transition probability q in the block Θ --of the Markovian process {J t }. It corresponds to the probability of changing from the hard phase to the soft phase at any state of the Erlang time H and is given as follows.

Note by Z the number of Bernoulli trials with success probability q needed to get one change from the hard phase to the soft phase. The random variable Z is indeed geometrically distributed with mean 1 q . We set the expected value of Z equal to the number of stages of the Erlang time H at which we want the cycle to jump to the soft phase. We write this as

E[Z] = a × t 0
which ensures that the time spent in the hard phase is t 0 on average. The rate parameter a of the Erlang time H corresponds to the number of states needed to have a one year long time. We thus deduce the probability q as

q = 1 a × t 0 = T L × t 0 .
For a better overview, suppose that H ∼ Erlang(L, a) with L = 2 and that the inter-claim times W 1 and W 2 are identically distributed. Their canonical random variable W is distributed as a mixture of exponential distributions with intensity matrix K such as

K = -λ 1 0 0 -λ 2 .
The block Θ --becomes

Θ --=             -λ 1 -a 0 0 0 a 0 0 0 0 -λ 2 -a 0 0 0 a 0 0 0 0 -λ 1 -a 0 aq 0 a(1 -q) 0 0 0 0 -λ 2 -a 0 aq 0 a(1 -q) 0 0 0 0 -λ 1 -a 0 0 0 0 0 0 0 0 -λ 2 -a 0 0 0 0 0 0 0 0 -λ 1 -a 0 0 0 0 0 0 0 0 -λ 2 -a             = Λ ---aI M * 0 Λ ---aI .
The environment process can go from the hard phase to the soft phase without affecting the run of the Erlang time H. This can of course be extended for L 2 leading to a block matrix

Θ --=       Λ ---aI M * 0 0 • • • 0 0 0 Λ ---aI M * 0 • • • 0 0 • • • • • • • • • • • • • • • • • • • • • 0 0 0 0 • • • Λ ---aI M * 0 0 0 0 • • • 0 Λ ---aI       ,
where the block M * contains the transition probabilities from the hard phase to the soft phase.

Computational aspects

Ramaswami et al. ( 2008) proposed a recursion procedure to solve Equation (2) in η using its particular structure which permits to reduce the dimension of the matrices involved in Theorem 1 and gain significative computational time. Indeed, the upcrossing phase probability matrix η has the following upper triangular block Toeplitz structure

η =       η 0 η 1 η 2 • • • η L-1 0 η 0 η 1 • • • η L-2 0 0 η 0 • • • η L-3 • • • • • • • • • • • • • • • 0 0 0 • • • η 0       .
Denote that in our model, the block Θ --of the intensity matrix Θ can be expressed as

Θ --=       Λ ---aI M * 0 0 • • • 0 0 0 Λ ---aI M * 0 • • • 0 0 • • • • • • • • • • • • • • • • • • • • • 0 0 0 0 • • • Λ ---aI M * 0 0 0 0 • • • 0 Λ ---aI       ,
where the block M * contains the transition probabilities between the two phases as explained above. Following the idea of [START_REF] Ramaswami | The Erlangization method for Markovian fluid flows[END_REF], writing out Equation (2) leads to the recursion equations for the blocks of η

η 0 Λ +-η 0 + R -1 --(Λ ---aI) η 0 + η 0 Λ ++ + R -1 --Λ -+ = 0 (6) η l Λ ++ + l j=0 η j Λ +-η l-j + R -1 --(Λ ---aI) η l + R -1 --M * η l-1 = 0; l = 1, ..., L -1. (7)
Define for a matrix A the vector vec(A) as the vector obtained by concatenating successive columns of A from a single matrix. From [START_REF] Graham | Kronecker products and matrix calculus: with applications[END_REF], we can write for a matrix product AXB that vec(AXB) = (B t ⊗ A)vec(X), where B t is the transpose of the matrix B. As in [START_REF] Ramaswami | The Erlangization method for Markovian fluid flows[END_REF], we use this notation and Equations ( 7) to derive a recursive expression for η l , l = 1, ..., L -1 and obtain

vec(η l ) = [Λ ++ + Λ +-η 0 ] t ⊗ I + I ⊗ R -1 --(Λ ---aI) + η 0 Λ +- -1 × vec   -   R -1 --M * η l-1 + l-1 j=1 η j Λ +-η l-j     . ( 8
)
Then, we only need to solve the nonlinear equation ( 6) using Theorem 1 to approximate η 0 and use the recursive expression (8) for η l , l = 1, ..., L -1 in order to derive an approximation of η.

Control of precisions

The control of precisions of our approximation method is based on the following points : the number L of stages of the Erlang horizon H, the cycle change probabilities and the fitting algorithms that we use to obtain claim amounts distributed as mixtures of Erlang distributions..

Model errors

As shown is Theorem 6 of [START_REF] Asmussen | Erlangian approximations for finite-horizon ruin probabilities[END_REF], the following convergence result for the finite-time ruin probability stands:

ψ(u, H L ) = ψ(u, T ) + ψ (2) (u, T ) L2! + O(L -2 ),
where T > 0 and H L denote the Erlang distribution for H with L stages and mean T . The model error can be estimated using the numerical computation for the second derivative of ψ(u, T ) as

ψ (2) (u, T ) = ψ(u, T + h) + 2ψ(u, T ) -ψ(u, T -h) h 2 + O(h 2 ).
The same authors also suggest to use the Richardson extrapolation leading to the improved estimate

ψ(u, T ) ≈ (L + 1)ψ(u, H L+1 ) -Lψ(u, H L ) (9)
with an refined convergence rate L -2 .

Using probabilities p 1j = P r(Y 1j > x 0 ), j = 1, ..., n 1 , as probabilistic modeling for the cycle environment to change from the soft phase to the hard phase can provide some errors. Indeed, the closer to 0 or to 1 all the p 1j 's, the most efficient the approximation. So, the threshold x 0 must be appropriately fixed. Giving such phase change probabilities leads to consider the Y 1j as if they are separate random variables. Medium values for the p 1j 's includes some perturbations and the modeling of the general claim behavior can fail.

The change of the environment process from the hard phase to the soft phase as explained in the Subsection 4.2 introduces a possible change of phase even when the process enters in the hard phase in the t 0 last years of the horizon time T . There, the process should not be able to change the phase but setting different probabilities for the phase change breaks the repeating block structure of η and then the recursion equation ( 8) could not be used anymore. We suggest then to use the possibility for the environment process to go back to the soft phase only for a large horizon time T and a large ratio T /t 0 .

An other part of the potential error comes from the fitting algorithm. To ensure a good fit, several statistical tools such as Chi-squared test, Kolmogorov-Smirnov test, qq-and pp-plots can be used.

Numerical errors

Some numerical issues can also include some deviation as the matrix inversion required in Equation (3). This matrix inversion is done using the matrix right division operator "/" in Matlab R2007a software. Furthermore, we use the software package Expokit from R. B. Sidje (see [START_REF] Sidje | Expokit. A software package for computing matrix exponentials[END_REF]) which uses Krylov subspace projection methods in order two compute the matrix exponential of Equation ( 4). The resort to this numerical tool can also include a slight deviation but permits to save substantial computation time.

Numerical applications

We illustrate in this section the use of our model when it is considered that the phase change instants depend on the claim amounts. As discussed in introduction, the European Solvency II regulation project will impose the insurance and reinsurance companies to have a 1-year horizon time ruin probability less than a certain threshold which is fixed at 0.5%. For reinsurance companies dealing with large claims and catastrophic events, it seems more realistic to consider a 5-or 10years horizon time with a ruin probability respectively around 2.5% or 5%. We propose for this application to deal with a 5-years period with a ruin probability around 2.5%. We also suppose that the environment process begins in the soft phase and does not go back to it. This last hypothesis of not going back to the soft phase could be seen as restrictive but is in fact quite realistic. Indeed, our studying time horizon is smaller than the average underwriting cycle length of 6 years as explained in [START_REF] Weiss | Underwriting cycles: a synthesis and further directions[END_REF], [START_REF] Cummins | An international analysis of underwriting cycles in propertyliability insurance[END_REF] or [START_REF] Venezian | Ratemaking methods and profit cycles in property and liability insurance[END_REF]. The phases differ here only by their premium income rates, a low one for the soft phase and a high one for the hard phase.

Data simulation

We first simulate 5000 realizations from a Generalized Pareto Distribution (GPD) with cumulative distribution function (c.d.f.), for x > 0,

F GP D (x) = 1 -1 + ξx σ -1 ξ ,
where σ > 0 is the scale parameter and ξ > 0 the shape parameter. This last parameter is also referred as the extreme value parameter in extreme value theory where α = 1/ξ is called the index of the regular variation of the distribution. This distribution is heavy-tailed and permits to modelize large claims and catastrophic events as fire losses, wind storm, earthquake, etc. The smaller the index α and, consequently, the greater the shape parameter ξ, the heavier the right queue of the distribution. For this illustration, we set α = 2.5 (ξ = 1/2.5) and σ = 1. A histogram for the simulated data is given in Figure 2. The empirical mean and variance are here respectively 1.63130 and 11.06983.

Figure 2: Simulation of claim amounts from a GPD

Phase-type fitting

Our model requires claim amounts to be phase-type distributed. There are different propositions that have be done in the literature to fit phase-type distributions to heavy-tailed distribution or data sets. One could object that this does not allow for heavy tails. However, as we are interested in ruin probabilities of order 2.5% or 5% (for time horizon 5 years or 10 years), we are not necessarily in the domain of asymptotic approximation and it is possible to have an appropriate fit of the claim size distribution up to the required level in our case. If we wanted to study ruin probabilities of magnitude 10 -4 , this method would of course be inappropriate, and one has to be very cautious when using light-tailed models in risk theory. For example, Feldmann and Whitt (1998) proposed a heuristic-based approach to fit heavy-tailed distribution functions. But this necessitates to first fit an heavy-tailed distribution to the data set which introduces additional errors. The popular expectation-maximization (EM) algorithm permits to fit directly data sets. Several variants of fitting methods using the EM algorithm to obtain a general or a particular phase-type distribution have been explored as in [START_REF] Asmussen | Phase-type distributions and risk processes with statedependent premiums[END_REF], [START_REF] Khayari | Fitting world-wide web request traces with the EM-algorithm[END_REF] or [START_REF] Riska | An EM-based technique for approximating long-tailed data sets with PH distributions[END_REF]. For this application, we use the EM algorithm developed in [START_REF] Lee | Modelling insurance losses and calculating risk measures via a mixture of Erlangs[END_REF] to fit mixtures of Erlang distributions with common scale parameter to heavy-tailed data sets. As written in this reference, this class of mixtures of Erlang distributions with common scale parameter is dense in the space of positive continuous distributions and then can approximate any positive absolutely continuous distribution to any accuracy. Moreover, this class is closed in mixing, convolution, compounding and aggregation of data which offers a good tractability for actuarial purposes. [START_REF] Wang | A general model for long-tailed network traffic approximation[END_REF] built an EM fitting algorithm for mixtures of Erlang distributions with common shape parameter and [START_REF] Thummler | A novel approach for phase-type fitting with the EM algorithm[END_REF] with no required common parameter. But these two classes does not fulfill the closure properties mentioned above.

The density of a mixture of k Erlang distributions with common scale parameter θ > 0 is

f (x) = k j=1 α j x r j -1e -x/θ θ r j (r j -1)! ,
where r j > 0 and α j > 0 are respectively the shape parameter and the weight of the jth Erlang distribution. The distribution weights must observe the following constraint k j=1 α j = 1.

Applying the Lee & Lin algorithm to our simulated data set, it results that a density using a mixture of 14 Erlangs fits the empirical density very satisfactorily. The common scale parameter is evaluated at θ = 0.81585 and the shape parameters and weights for every Erlang distributions of the mixture are presented in To attest to the goodness of fit, a graphical comparison between a relative histogram for the simulated data and the fitted distribution can be seen in Figure 3, a pp-plot and a qq-plot are drawn in Figure 4. All these graphical tools already show that the fit is very satisfying even in the right queue. A Kolmogorov-Smirnov test is used to confirm the high quality of the fit. Indeed, with a test statistic of 0.00578 and a cutoff value of 0.01917, the p-value is evaluated at 0.99612 and the hypothesis that the data sample comes from the fitted distribution is not rejected. Again, this method could be dangerous if it were applied to the study of more extreme events.

Ruin probability approximation

We now use our model to evaluate the ruin probability for a horizon T = 5 years where the claim amounts are distributed according to the fitted mixture of Erlang distributions. Let the inter-claim time random variable W be exponentially distributed with parameter λ = 2. So that, the average time between two claim arrivals is 0.5 year and two claims per year are expected. As for the claim amounts, we suppose that the inter-claim time r.v. are identically distributed in the two phases of Only the premium income rates differs between the two phases. Indeed, we define the premium rate c i , with i = 1 for the soft phase and i = 2 for the hard phase, such as

c i = (1 + ρ i ) × E[ S 1 ] = (1 + ρ i ) × E[W ]E[ X 1 ],
where ρ i is the security loading corresponding to the phase i, S t is the surplus process where the claim amount random variables X k are i.i.d. with the empirical simulated distribution. [START_REF] Froot | The market for catastrophe risk: a clinical examination[END_REF] shows that in the market of large claims and catastrophic events the reinsurance premiums can be more than seven times the expected loss during a period. Following this idea, we set ρ 1 = 30% and ρ 2 = 300%. We fix the initial surplus at u = 30 and the claim amount threshold x 0 > 0 that provokes the passage from the soft phase to the hard phase at x 0 = 25. This last value represents the 99.663% quantile of the empirical cumulative distribution. Over the 5 years horizon with 2 claims per year, this does not seem to have an impact the ruin probability. In Table 3 is displayed the probability p 1j of causing a change from the soft cycle to the hard cycle for the jth Erlang distribution of the mixture, j = 1, ..., 14, as they are defined in Equation ( 5 We can now construct the required matrices Λ, R --and M * and use Equations ( 6), ( 8) and Theorem 1 to approximate the finite-time ruin probability ψ(u, T ). For comparison issues, Monte-Carlo (MC) simulations for ψ(u, T ) are performed using both empirical and fitted distributions for the claim amounts. The results for the three approximations are given in Table 4 where we used for our method a number of Erlang stages L = 50.

Method

Approximation for ψ(u, T ) Model 3 2.37199% MC sim. with fitted distr.

2.35448% MC sim. with empirical distr.

2.34998% Table 4: Approximated finite time ruin probabilities

We see on this Table that our approximation is quite accurate with a relative error of 0.73% from the Monte-Carlo simulation with the fitted distribution and 0.92% from the Monte-Carlo simulation with the empirical distribution. The error inserted from the distribution fitting is 0.19%. The computational time required for our model is here around 100 times less than the one needed for the two simulation methods.

The previous example is the reference parametrization that is used in the next subsection.

But, we also give now some results in

Comparison between model 3 and independent changes of cycle phase

We illustrate now the impact of dependence between phase changes and claim amounts. We compare our model to a so-called independent model in which the change of phase occurs independently from claim amounts. For this independent model, we construct an intensity matrix Θ for the environment Markov process J t in which the change from the soft phase to the hard phase can appear at any state of the Erlang horizon time H with probability p 0 ∈ [0, 1]. This is done by modifying the block matrix Θ --such that at any state i, i = 1, ..., L -1 of the Erlang time in the soft cycle phase, the process can go to state i + 1 of the Erlang time in the hard cycle phase with probability p 0 . Therefore, the change of phase is independent from the claim amount.

For exponentially distributed inter-claim times with parameter λ and for an Erlang time horizon H ∼ Erlang(2, a), the block Θ --becomes

Θ --=     -λ -a 0 a(1 -p 0 ) ap 0 0 -λ -a 0 a 0 0 -λ -a 0 0 0 0 -λ -a    
and can also be written as

Θ --= Λ ---aI M * 0 Λ ---aI ,
where, as in Subsection 4.2, M * contains the direct transitions between the phases. This can, of course, also be extended to a number of Erlang stages L 2.

In order to compare this independent model to our dependent model, we give to these two models the same probability to change at least one time from the soft phase to the hard phase. For the independent model, this probability is P 0 = 1 -(1 -p 0 ) L . In our dependent model, the probability to go to the hard phase after each claim is n j=1 α j p 1j .

Then the probability to have no change from the soft phase to the hard phase during [0, T ] is

E[(1 - n j=1 α j p 1j ) N ]
where N ∼ P oisson(T λ). Finally, the probability to change from the soft phase to the hard phase at least one time during [0, T ] is

P = 1 -E[(1 - n i=1 α i p i ) N ] = P 0
and we obtain

p 0 = 1 -E[(1 - n i=1 α i p i ) N ] 1/L .
We also construct an intermediate model such that the probability to change from the soft phase to the hard phase at least one time is still the same where the cycle can change either independently from claim amounts or because of a large claim. This intermediate model permits to observe the slide between the dependent and independent models and to see the impact of the claim amount dependence on the finite-time ruin probability. The probability to change from the soft phase to the hard phase at least one time during [0, T ] for this intermediate model is

P (θ) = 1 -(1 -p(θ)) L E (1 - n i=1 θα i p i ) N = P 0
where θ ∈ [0, 1] and p(θ) represents the part of the phase change probability that is independent from the claim amounts. Then we obtain

p(θ) = 1 -(1 -p 0 ) × E (1 - n i=1 θα i p i ) N -1/L .
The bounds p(0) = p 0 and p(1) = 0 correspond respectively to the independent model and the dependent model.

For our application, as we use the Richardson extrapolation recalled in Equation ( 9) for a number of Erlang stages L = 50, two probabilities p 0 are evaluated. For L = 50, we obtain p 0 = 6.63365×10 -4 and for L = 51, p 0 = 6.50362×10 -4 . Thus, the ruin probability is approximated at 2.50960% with the independent model versus an evaluation at 2.37199% for our dependent model. We recall that the MC simulation with the empirical distribution provides a finite-time ruin probability at 2.34998%. The error from the independent model on our dependent model can be evaluated at 5.48% and the one on the MC-simulation at 6.36%. We can imagine that this difference would be quite larger when looking at a greater time horizon but also when including more phases in a cycle.

Conclusion

This paper provides an extension for the Erlangization method taking into account the reinsurance cycles in order to evaluate the finite-time ruin probability in the Sparre-Andersen model with phasetype distributed claim amounts and inter-claim times. We assume the presence of two phases in a reinsurance cycle, a soft phase and a hard phase. We first propose a naive model with deterministic changes between the two phases. Then, a more realistic model is developed where the change from the soft phase to the hard phase depends on the severity of the claim amounts. We illustrate our model by an application and compare our results to Monte-Carlo simulations and to the case where the claim amounts have no impact on the phase changes. We show that this feedback control may change the ruin probability by 6% in our particular numerical example.

where .

M 1 =       B (n) 1 B (n) 2 . . . B (n) L-1 0 B (n) 1 . . . B (
Thus we have

(B (n) ) -1 = M -1 1 -M -1 1 M 2 M -1 4 0 M -1 4 .
The remaining difficulty stays in the inversion of the block matrix M 1 . To this aim, we proceed algorithmically using the same procedure as for B (n) .

B Initialization for the deterministic phase change instants model

For the computation of the approximation of the ruin probability in the case of deterministic phase change instants, we have some technical but not stopping constrains on the numbers of states L 1 for H 1 and L 2 for H 2 . Indeed, we have

H 1 ∼ Erlang(L 1 , a), H 2 ∼ Erlang(L 2 , a).
The final horizon time T can be decomposed as T = T 1 + T 2 such as

T 1 = E[H 1 ] = L 1 a , T 2 = E[H 2 ] = L 2 a .
Both T 1 and T 2 are fixed and then we have to find appropriate values for a ∈ R * + , L 1 ∈ N et L 2 ∈ N such that

L 1 = aT 1 , L 2 = aT 2 , L = L 1 + L 2 = a(T 1 + T 2 ) = aT.
These equalities can be fulfilled by choosing reasonable values for L, L 1 and L 2 knowing that improving this values leads to better approximations.

Figure 1 :

 1 Figure 1: Example of a cycle path with 10 phases

Figure 3 :

 3 Figure 3: Comparison between a relative histogram for the simulated data and the fitted density

Figure 4 :

 4 Figure 4: pp-plot and qq-plot for the simulated data versus the fitted distribution

  Figure 5 are drawn the slide of the ruin probability with the intermediate model and the slide of the underlying error made comparing to our dependent model.

Figure 5 :

 5 Figure 5: Comparaison between the model 3 and the independent model

Table 1 :

 1 Example of phase change probabilities at phase 3.

	Phase 2	0.001	Phase 2	0.09
	Phase 3	0.9	Phase 3	0.9
	Phase 4	0.09	Phase 4	0.009
	Phase 5	0.009	Phase 5	0.001

Phase 3 ↓ to Probability Phase 3 ↑ to Probability

Table 2 .

 2 This table also displays the expected value associated for each Erlang r.v. Y j of the distribution mixture.

	j r j	α j	E[Y j ]	j r j	α j	E[Y j ]
	1 75 0.00063 61.18865	8 22 0.00122 17.94867
	2 59 0.00021 48.13507	9 16 0.00469 13.05358
	3 48 0.00012 47.31922 10 15 0.00283 12.23773
	4 40 0.00199 32.63395 11 14 0.00166 11.42188
	5 39 0.00024 31.81810 12 8 0.03157 6.52679
	6 25 0.00078 20.39622 13 4 0.14131 3.26339
	7 24 0.00122 19.58037 14 1 0.81155 0.81585

Table 2 :

 2 Parameters of the fitted mixture of Erlang distributions.

Table 3 :

 3 ). Phase change probabilities of the Erlang distributions of the fitted mixture

	j	p 1j	j	p 1j
	1	≈ 1	8 0.04326
	2	≈ 1	9 0.00138
	3 0.99999 10 0.00064
	4 0.94060 11 0.00028
	5 0.91826 12	≈ 0
	6 0.13159 13	≈ 0
	7 0.09433 14	≈ 0

Table 5 :

 5 Table 5 when modifying one parameter. We see that our approximations are sill very close to the Monte-Carlo simulations obtained with the empirical claim size distributions. Approximated finite time ruin probabilities with a modified parameter

	Parameter modified MC simulation Model 3 Relative error
	x 0 = 30	2.42243%	2.42480%	0.098%
	ρ 2 = 200%	2.38310%	2.40007%	0.707%

Acknowledgements

This work has been partly funded by the research chair Management de la Modélisation sponsored by BNP Paribas Assurance, the research chair Actuariat Durable sponsored by Milliman, and the French Research National Agency (ANR) under the reference ANR-08-BLAN-0314-01.

APPENDIX A Matrix inversion

The Erlangization method used to evaluate the ruin probability or equivalently the probability for the fluid model to reach 0 before the horizon H requires the inversion of several matrices. As seen in Theorem 1, the matrix of interest η can be approximated by the iteration scheme η = lim n→∞ η (n) where η (0) = 0 and

Thus, to determine η (n+1) , we need to inverse the matrice B (n) . The Erlangization method, according to the specific structure of the intensity matrix of the Erlang horizon H, provides that the block matrices Θ --, Θ -+ , Θ +-and Θ ++ have all an upper triangular block Toeplitz structure and that the matrix B (n) takes this form

Indeed, this structure is stable for addition, multiplication and inversion and B (n) -1 , A (n) and η (n) also have the form.

Let M be a mn × mn matrix such as

where the m × m block matrix M 1 and the n × n block matrix M 4 are invertible. Then, we have

.

To inverse the matrix B (n) described before, we can split it into four blocks as follows: