Victoria Fernández
email: victoria.fernandez-abrevaya@inria.fr

Abrevaya Stefanie

Wuhrer Edmond Boyer

Spatiotemporal Modeling for Efficient Registration of Dynamic 3D Faces

We consider the registration of temporal sequences of 3D face scans. Face registration plays a central role in face analysis applications, for instance recognition or transfer tasks, among others. We propose an automatic approach that can register large sets of dynamic face scans without the need for landmarks or highly specialized acquisition setups. This allows for extended versatility among registered face shapes and deformations by enabling to leverage multiple datasets, a fundamental property when e.g. building statistical face models. Our approach is built upon a regression-based static registration method, which is improved by spatiotemporal modeling to exploit redundancies over both space and time. We experimentally demonstrate that accurate registrations can be obtained for varying data robustly and efficiently by applying our method to three standard dynamic face datasets.

Introduction

Facial motion provides significant non-verbal signals that give cues for communication and social interaction. Considering such dynamic face information expands the scope and performance of automatic facial analysis systems, for instance with expression or identitiy recognition [START_REF] Fang | Kakadiaris. 3d/4d facial expression analysis: An advanced annotated face model approach[END_REF][START_REF] Alashkar | A grassmann framework for 4d facial shape analysis[END_REF], pain detection [START_REF] Zhang | Three dimensional binary edge feature representation for pain expression analysis[END_REF] or realistic expression synthesis [START_REF] Yu | Perceptiondriven facial expression synthesis[END_REF]. To recover coherent face information from inconsistent visual observations, the data must generally be first registered to ensure that anatomically corresponding points are consistently identified. Such registration can be purely spatial, yielding a static face pose parametrization that is agnostic to time information. Yet, in the case of faces in motion, the registration can account for the temporal aspects through a spatiotemporal face pose parametrization. The interest is to better capture face deformations with a temporal tracking, where static registrations provide only coarse and noisy motion information. Spatiotemporal, or dynamic, registration is however more complex than its static counter-part since tracking robustly and reliably is still challenging in practice.

Face registration is also confronted with the issue of possibly large and heterogeneous datasets. This is particularly true with recently released large 4D face databases, e.g. [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF][START_REF] Zhang | Bp4d-spontaneous: a high-resolution spontaneous 3d dynamic facial expression database[END_REF]. While providing a rich and diverse source of information for the study of both spatial and temporal patterns in facial motion, the body of available datasets raises difficulties for global registration, such as dealing with data coming from different capturing systems and different acquisition protocols. Face registration methods in this context should be able to handle datasets in a fully automatic way, as manual intervention is not feasible for thousands of frames, in addition to being robust to different types of noise resulting from different acquisition scenarios.

This work addresses these objectives by proposing an automatic method to register dynamic 3D face data with the longer-term motivation of analyzing large datasets of facial motion. The method is purely geometry-based to allow leveraging any available temporal 3D face scan, even when RGB data is not provided e.g. for privacy reasons. We do not require pre-determined landmarks as input, which in the case of 3D data is more challenging to obtain, thus removing a possible source of error. The main innovation is to extend a recent efficient and robust static 3D face registration method [START_REF] Abrevaya | Multilinear autoencoder for 3d face model learning[END_REF] to dynamic data by exploiting the spatiotemporal coherence [START_REF] Akhter | Bilinear spatiotemporal basis models[END_REF] of the data. This enables registrations that both fix identities over temporal sequences and regularize observed motions to prevent high-frequency flickering. The approach presents the following advantages: it can register multiple datasets into a single representation; it does not require color information, as in e.g. [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF]10,14], and is robust to occlusions by construction; it runs an order of magnitude faster than recently proposed methods based on parametric face models [7, 20] while achieving comparable accuracy; and provides compact representations of the results.

We evaluate our method qualitatively and quantitatively on the three publicly available datasets D3DFACS [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF], BU-4DFE [31] and BP4D-Spontaneous [START_REF] Zhang | Bp4d-spontaneous: a high-resolution spontaneous 3d dynamic facial expression database[END_REF], and show that we can efficiently obtain accurate representations. We provide comparisons to three recent methods and show that we achieve similar or better results in terms of vertex-to-scan error and in the semantic parametrization, while remaining either more general in terms of requirements of the datasets, or more efficient in terms of computational times.

Related Work

Numerous works have studied the registration of static 3D face scans [6,[START_REF] Patel | 3d morphable face models revisited[END_REF]. While a static method can be applied independently to each frame of a motion sequence, this is known to lead to artifacts including high-frequency jitter. We consider here 4D face data in the form of sequences of 3D face scans and therefore focus our review on methods that take advantage of the temporal redundancy captured by 4D data. Another line of research that has recently received considerable attention is the reconstruction of 4D facial motion based on monocular 2D video. These works solve an underconstrained reconstruction instead of a 3D registration as addressed in this work; the interested reader is anyway referred to the survey of Zollhöfer et al. [38]. In the following we discuss strategies for the registration of 4D face data.

Registration of 4D face data

Sun et al. [START_REF] Sun | Tracking vertex flow and model adaptation for three-dimensional spatiotemporal face analysis[END_REF] propose a geometry-based coarse-to-fine approach to register planar meshes using conformal mapping, aided by automatically computed feature points. For expression recognition, Fang et al. [START_REF] Fang | Kakadiaris. 3d/4d facial expression analysis: An advanced annotated face model approach[END_REF] perform pairwise registration of consecutive frames using an Annotated Face Model (AFM) [START_REF] Kakadiaris | Threedimensional face recognition in the presence of facial expressions: An annotated deformable model approach[END_REF]. Temporal information is exploited by initialization using the result of the previous frame.

For real-time expression transfer, Weise et al.

[30] introduce a system based on a non-rigid Iterative Closest Point (ICP) method, from which a person-specific blendshape model is built and used to track sequences of the same actor. Follow-up work [START_REF] Bradley | High resolution passive facial performance capture[END_REF] improve on this by using color cues and a probabilistic animation prior to handle noise. The methods of [19,[START_REF] Bouaziz | Online modeling for realtime facial animation[END_REF] further remove the need for calibration by updating an initial blendshape model on-the-fly. Recent improvements on this line of work include robustness to occlusions and pose [START_REF] Hsieh | Unconstrained realtime facial performance capture[END_REF] and more detailed blendshape models by the use of displacement maps [START_REF] Blanz | A morphable model for the synthesis of 3d faces[END_REF].

More recent methods follow two main lines. The first performs registration in texture space by computing correspondences between sparse landmarks predicted using an Active Appearance Model (AAM), which are densified using thin-plate spline deformations [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF]10]. The method of Cosker et al. [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF] achieves inter-sequence correspondence by registering each frame towards a manually selected neutral expression, and intra-sequence correspondence by registering these neutral frames to a template. To better handle texture variations, Cheng et al. [10] extend the previous by using session-and-subject specific AAM, and non-rigid ICP [4] between manually selected neutral frames. Since these methods operate on color information, they require careful acquisition setups with controlled lighting conditions, as e.g. moving shadows can lead to inaccuracies.

The second line of work takes advantage of low dimensional parametric shape spaces learned from large databases of static 3D face scans and used as prior during registration, e.g. [START_REF] Thies | Real-time expression transfer for facial reenactment[END_REF][START_REF] Zollhöfer | Automatic reconstruction of personalized avatars from 3d face scans[END_REF]7,20]. Most related to our work, multilinear models of identity and expression [7] and a linear articulated model with expressions [20] have been used for this purpose. These works achieve registrations of relatively high accuracy and report running times of 30s to 2min per frame. Although the overall facial shape is recovered, finescale details such as wrinkles are not modeled. Our work shares this property, as well as the robustness and accuracy of these methods while allowing for a gain in efficiency.

Joint registration and reconstruction

In the context of multi-view reconstruction, several works use optical flow to recover a temporally consistent geometry from synchronized multi-view 2D videos, e.g. [START_REF] Zhang | Spacetime faces: high resolution capture for modeling and animation[END_REF]9,[START_REF] Beeler | High-quality passive facial performance capture using anchor frames[END_REF]. Bradley et al. [9] jointly solve for registration and reconstruction by sequentially tracking the initial frame. Beeler et al. [START_REF] Beeler | High-quality passive facial performance capture using anchor frames[END_REF] improve on this method by using a set of "anchor frames", which are automatically selected key-frames that define subsequences in which optical flow is performed; results are of very high quality and achieve pore-level details. Both these methods require a dense setup of synchronized video cameras. Valgaerts et al. [27] simplify these requirements by introducing a method that achieves results of similar quality from a single pair of stereo cameras, combining scene flow to compute the global registration with shading-based refinement to compute fine-scale details. More recently, Fyffe et al. [14] proposed a method to reconstruct a registered model that includes the full head, eye balls, and the interior cavities of mouth and nose. All of these methods achieve temporally coherent results which are of high quality and include fine-scale details. However, they are limited to specific acquisition setups as the input to the methods are synchronized and calibrated 2D videos. In this work, we consider the more general problem of registering the geometry of 4D face scans without the need for reliable color information.

Background

This section recalls the two components our method builds upon: a multilinear autoencoder for static registration and a spatiotemporal model to handle time aspects.

Multilinear autoencoder

Multilinear models are suitable for 3D face modeling and have been shown to effectively separate identity and expression [28]. These models are based on tensor algebra, and obtain decoupled representations by parameterizing a face using a set of weights for each modeled factor of variation, e.g. identity and expression, whose relationship is encoded using a core tensor M. More formally, multilinear models allow to approximate a registered face x ∈ R 3n , stored as vector of coordinates associated with n vertices of a mesh, using identity and expression coefficient vectors

w id ∈ R d id and w exp ∈ R dexp by x ≈ x + M × 2 w id × 3 w exp , (1)
where M ∈ R 3n×d id ×dexp is the core tensor of the model, x is the mean face over the model's training data, and × i denotes mode-i multiplication (we refer to [18] for an overview of tensor operations). Typically, the dimensions d id and d exp are chosen to obtain a compact representation that yields a low approximation error of the training data. Abrevaya et al. [START_REF] Abrevaya | Multilinear autoencoder for 3d face model learning[END_REF] recently proposed an autoencoder architecture that allows to robustly and efficiently regress raw 3D face scans on identity and expression using a multilinear model. This method can be used to register an arbitrary 3D facial scan by representing it as a heightmap, which is processed by a convolutional neural network in order to get multilinear model parameters. We leverage this to obtain initial registrations for each frame of an input sequence.

Bilinear spatiotemporal model

Akhter et al. [START_REF] Akhter | Bilinear spatiotemporal basis models[END_REF] propose the use of a bilinear model to leverage both spatial and temporal redundancies of a single 3D motion sequence. The model requires that the sequence is registered, and we denote its F frames by x 1 , . . . , x F ∈ R 3n in the following. The temporal sequence is organized into a matrix S ∈ R F ×3n containing each frame in a row, S = [x 1 , . . . , x F ] T , which implies that the columns of S describe vertex trajectories. Let B ∈ R 3n×ds contain a set of shape basis vectors that encode each frame into a space of dimension d s , and similarly, let Θ ∈ R F ×dt contain a set of temporal basis vectors that encode the trajectory of each vertex into a space of dimension d t . Then the sequence matrix S can be decomposed as

S ≈ ΘCB T , (2)
where C ∈ R dt×ds is a matrix of spatiotemporal coefficients that compactly encode S. The dimensions d s and d t allow to trade off the compactness of the representation and the approximation error of the input sequence.

In the original formulation, both shape basis B and coefficients C need to be optimized for each individual sequence. However, they show that the temporal basis Θ can be fixed to the Discrete Cosine Transform (DCT), as this approaches the optimal PCA-learned basis when the data is generated from a stationary first-order Markov process. Akhter et al. [START_REF] Akhter | Bilinear spatiotemporal basis models[END_REF] empirically demonstrate this to hold for sparse facial data. We leverage this model to build a compact spatiotemporal space that allows to efficiently register 4D faces in an iterative way.

Method

Our goal is to register a large number of sequences of 3D face scans. Each of these sequences may contain many frames, and each frame a large number of vertices, making the problem high-dimensional and difficult to optimize. Furthermore, datasets captured with different acquisition setups present different levels of noise, missing data and occlusions, which plagues naive template fitting. To keep the method as general as possible, we do not assume availability of either landmark or color information.

To tackle this challenging problem, we proceed in two major steps. First, we perform an initialization independently on each frame that robustly and efficiently regresses each scan against identity and expression coefficients of the multilinear model. The resulting approximations correctly capture the general structure of the motion and shape but are still overly smooth and, because each frame is treated independently, exhibit high-frequency jitter. To remedy this, in a second step, we use the multilinear face model to build a bilinear spatiotemporal model, which iteratively improves the initial approximation. This regularizes the motion of the vertices, and turns the problem into a much lower-dimensional one compared to a frame-wise formulation. After a few iterations we obtain registered data which is regularized both in space and time and fits more closely to the scans. Figure 1 summarizes this process.

Frame-wise initialization

Given a sequence of observations [o 1 , . . . , o F ] consisting of F frames, we register each of the frames independently by regressing identity and expression coefficients using the multilinear autoencoder, as outlined in Section 3.1. We thus obtain a sequence of identity and expression weights that represent the face motion:

W id = [w 1 id , . . . , w F id ], and (3)
W exp = [w 1 exp , . . . , w F exp ]. (4)
Pre-processing Although the neural network we use [START_REF] Abrevaya | Multilinear autoencoder for 3d face model learning[END_REF] was trained with raw scans presenting different types of noise, the quality of the registration degrades when the input data differs in form and orientation from the original training data. To ensure results of high quality, we therefore pre-process the face as follows. We first detect the nose tip by training a neural network for this task on depth data with the same architecture and training data as [START_REF] Abrevaya | Multilinear autoencoder for 3d face model learning[END_REF]. We next crop the face a radius of 100mm around this point, and perform a coarse frontalization step so that it approximately looks towards the z-direction. For frontalization, we consider the direction of the normals of each vertex, which gives a distribution of orientations over the sphere that sample a semisphere for a cropped face. The directional mean of this distribution gives a coarse approximation of where the face is "looking at", and we align this to the directional mean of the model's mean face. This makes the weak assumption that the face is not upside-down and works well as long as the cropped face does not contain too many holes or extra parts. Note that only a coarse alignment is required here since the autoencoder was trained to have some degree of robustness to pose variation in the depth image. As a result, this preprocessing step can be replaced with any other method that will produce a cropped face and a rough frontalization.

The resulting multilinear model representations [x 1 , . . . , x F ], obtained by reconstructing the faces using Equation 1 with the coefficients from Equations 3 and 4, are in the coordinate system in which the multilinear model was learned. For further refinement of these approximations they need to be compared to the original scans. To align the observations [o 1 , . . . , o F ] to the model coordinate system, we take advantage of the height-map images generated during regression to find initial correspondences. In particular, we consider the depth image of the cropped and frontalized scan and the depth image of the registered mesh, and establish preliminary correspondences by assigning pixels at the same location. This correspondence is used to rigidly transform o i to x i (i = 1, . . . , F ). Subsequently, a few iterations of regular iterative closest point alignment are performed. Once each frame is aligned, we discard the cropped version and go back to the original raw scan; this allows to remove the quality of the crop as possible source of error in the subsequent steps.

Iterative spatiotemporal registration

The previous section computes, for each frame i, spatially aligned observations o i along with registered faces x i which approximate the geometry of o i . We leverage spatiotemporal cues to refine x i jointly over all frames, thereby obtaining temporally smooth results that can be represented compactly. By iterating between refining the geometry of x i to match o i , and building a spatiotemporal sequence representation, we improve the geometric approximation of the registrations.

Geometric refinement

To improve the quality of the approximation x i of o i , we non-rigidly deform the registrations to the scans. The following discussion omits the frame index i to simplify notation. The registration x is warped to o by optimizing for displacements δ x of the vertices of x along their normal directions with Laplacian regularization. In particular, we solve for

min {δ x 1 ,...,δ x n } α n j=1 wj||v x j + δ x j n x j -p x j || 2 + β n j=1 ||L(δ x j )|| 2 , (5)
where v x j ∈ R 3 is a vertex in frame x, p x j ∈ R 3 the closest point to v x j in o, n x j ∈ R 3 the normal vector of v x j , L the cotangent discretization of the Laplace-Beltrami operator [START_REF] Meyer | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF], and w j , α, β are scalar weights. We discard closest points whose Euclidean distance is greater than 5mm and whose deviation in normal vector is greater than 45 • by setting w j to zero. This formulation can be efficiently minimized by solving a linear system of equations.

Spatiotemporal sequence projection

We build upon the model explained in Section 3.2 to compute a spatiotemporal representation of the input sequence for given [x 1 , . . . , x F ]. Unlike the original formulation of [START_REF] Akhter | Bilinear spatiotemporal basis models[END_REF], in which both the shape basis B and the model coefficients C need to optimized, we leverage the multilinear model to obtain a person-specific spatiotemporal representation. Specifically, we summarize the regressed identity coefficients W id into a unique coefficient for the entire sequence, given that the identity of the subject is fixed for any given motion. Let w id denote this unique coefficient vector. We compute w id as a mean over the regressed results, w id = W id . Given w id we can create a shape basis B by multiplying the core tensor with w id as B T = M × 2 w id . For the temporal basis, we follow Akhter et al. and fix Θ to the DCT basis.

For each iteration, we gather the approximations x i obtained in the previous step into a sequence matrix S, and use Equation 2 to compute C by solving

BC T = S T Θ. (6)
This can be performed efficiently, as B is fixed and can be factorized once for all the iterations.

Sequence representation

To obtain more detailed results we complete the iteration with a final geometric refinement step. This allows the results to leave the bounds of the multilinear model, thereby providing more accurate approximations of o i . To prevent artifacts, a stronger regularization weight β is used in the last geometric refinement. This allows closer fits towards the data, but in turn looses both the compactness of the representation and the motion regularization. To rectify this, we project the trajectory of the displacements δ x i into a second DCT basis (of possibly different dimensionality), obtaining a displacement coefficient vector d i ∈ R d t for each of the n vertices, with d t << F . We thus retain compactness while allowing for more detailed registrations, as well as preventing flickering in the final trajectory of the vertices.

After this process, we can store the results compactly given the multilinear model, which is stored once for all sequences. For each registered sequence the following information suffices to reconstruct x i : (1) the identity coefficient w id ∈ R d id that compactly encodes the shape matrix B; (2) the spatiotemporal coefficients C ∈ R dt×ds ; (3) the dimensions of the temporal basis d t and d t ; and (4) the displacement coefficients {d 1 , . . . , d n }, d i ∈ R d t . This significantly reduces the storage requirements of large datasets (e.g. from 9.1GB to less than 1MB in the example from Figure 4), while still retaining a reasonable level of detail.

Evaluation

We validate our method by registering D3DFACS [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF], BU-4DFE [31] and BP4D-Spontaneous [START_REF] Zhang | Bp4d-spontaneous: a high-resolution spontaneous 3d dynamic facial expression database[END_REF] datasets. D3DFACS contains 519 sequences of 10 subjects performing different types of facial action units, while BU-4DFE presents 101 subjects with 6 sequences each performing the six prototypical emotions; in both cases the average sequence length is around 100 frames and the meshes contain around 30K and 35K vertices respectively. BP4D-Spontaneous is made of 328 sequences with 41 subjects performing 8 tasks each, which were designed to elicit spontaneous emotions. The average sequence length is around 1100 frames and the average number of vertices is around 37K. In the following we provide both qualitative and quantitative evaluations over these.

Implementation details

The code was implemented in C++ using the Eigen3 library [START_REF] Guennebaud | [END_REF]. We use the publicly available autoencoder of [START_REF] Abrevaya | Multilinear autoencoder for 3d face model learning[END_REF] that was trained on BU-3DFE [START_REF] Yin | A 3d facial expression database for facial behavior research[END_REF] and Bosphorus [START_REF] Savran | Bosphorus database for 3d face analysis[END_REF] datasets. The dimensions of identity and expression spaces are set to 65 and 20. For the dimension of the temporal basis, we observe that the impact depends on the length of the sequence. We thus set it to a factor of the sequence length; in particular F/5. We set d t = 5, and unless otherwise specified, we fix the number of iterations to 5. For Equation 5we set β = 1, α = 0.9 during iterations and α = 0.8 for the final step. Since BU-4DFE presents noisier scans, we set α = 0.5 during iterations and α = 0.2 for the final step to avoid overfitting. The mesh used to register all face scans contains 5996 vertices.

Qualitative results

Figure 3 shows an example of the results obtained on each of the steps of the method: regression, spatiotemporal registration, and final refinement. Figure 2 shows a few more examples of registrations obtained on D3DFACS, BU-4DFE and BP4D-Spontaneous. They illustrate that accurate cross-dataset registrations can be obtained, while still being robust to different types of noise in the data.

Quantitative results

We evaluate the quality of the registrations with two commonly used metrics: median per-vertex error towards the input scan and landmark distances. The median pervertex error is taken across all registered frames in the dataset, and shows how close the registrations are to the real scans. We also evaluate semantic accuracy by manually placing landmarks on 5 key-frames over 10 randomly selected sequences of D3DFACS, and measuring the Euclidean distances between these and the landmarks defined over the template. We use in particular 11 landmarks: nose tip, inner and outer eye corners, left and right mouth corners, and upper/lower outer and inner center of lips. The chosen key-frames sample the sequence by taking the first and last frame, the peak frame, and two intermediate ones.

We also evaluate the stability of the motion by using a compactness measure [START_REF] Davies | Statistical models of shape: Optimisation and evaluation[END_REF] as follows. For each sequence, we align the frames using generalized Procrustes analysis, perform PCA, and measure the amount of variabilty captured by each principal component. If the registrations exhibit high-frequency jitter, we expect to see less variability retained by the first principal components, as the variations coming from flickering vertices would have to be encoded by higher-order principal components. To be able to summarize this over the entire dataset, we determine the mean variability obtained as a function of the percentage of principal components considered.

Figure 5 shows a cumulative plot of the median pervertex error on D3DFACS obtained after the initial regression, after 1, 3, 5 and 10 iterations of spatiotemporal registration, and for the final result. We can see that the iterative process improves the initial regressions in terms of surface fit. Furthermore, Table 1 shows the mean landmark error over the 11 landmarks for the final results obtained after 1, 3, 5 and 10 iterations. Despite being a landmark-free registration method, the method obtains a good semantic accuracy that is improved with each iteration.

We evaluate the benefits of the temporal regularization by comparing the full model with a static version of our framework. For this, instead of projecting onto the spatiotemporal model we independently project each frame onto the shape basis B, and measure the results in terms of vertex error and compactness. Figure 6 shows cumulative plots obtained for these registrations. Note that while using a spatiotemporal model achieves similar accuracies in terms of vertex error, the compactness of each sequence improves with the spatiotemporal model, implying less highfrequency jitter with the latter. This results can also be qualitatively assessed in the accompanying video.

Finally, we show the ability of the method to track long videos by registering the sequences from BP4D-Sponteanous, many of which consist of more than 1000 frames. We obtain a mean vertex error of 0.33mm over all registered sequences and frames. Figure 4 further shows the median per-vertex error for each frame on one example. This error stays between 0.1 and 0.4mm and does not increase with time, suggesting that no drift is occuring. This is expected, as the regression step is performed independently on each frame. A textured visualization of this sequence can be found in the accompanying video.

With respect to the running times, we report a mean perframe processing time of 578ms on the D3DFACS dataset, 637ms on BU-4DFE and 399ms for BP4D, for 5 iterations in all cases. Computation times were measured on an Intel Xeon 3.30GHz with NVidia GeForce GTX 1080 GPU.

Comparisons

We compare our method to the previous works of Bolkart and Wuhrer [7] , Li et al. [20] and Cosker et al. [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF] using registrations provided by the authors.

Bolkart and Wuhrer This work [7] also registers motion sequences in a fully-automatic manner by using a multilinear model and geometric information only. We compare to this method on 497 sequences from BU-4DFE, which are the sequences that were correctly registered by [7]. Figure 7a shows cumulative plots of the median per-vertex error on all the registered sequences in BU-4DFE, comparing [7] to our registration without and with the final refinement step, since [7] has no refinement step. Figure 8a further shows a qualitative comparison over a challenging example. Results reveal similar accuracy for both methods without the refinement step, whereas [7] requires around 30 seconds per-frame to process.

Li et al.

The method of [20] was used to register D3DFACS, and thus we compare our results over this dataset. For a fair comparison we crop their full-head model so that it contains only the face, to be similar to our registrations. We obtain a mean vertex error of 0.13 mm for our method, and 0.33mm for [20]. In Figure 7b we show the cumulative plots for the median per-vertex distance for both methods. They demonstrate that our approach achieves higher accuracy while reducing both the running time ([20] reported 155 seconds per-frame) and the requirements on the dataset. We also compare with respect to landmark errors; results can be found in Table 1. Our approach achieves better results with a single iteration, even though it requires no pre-determined landmarks to guide the process. This also confirms that our registrations faithfully preserve the anatomic semantics. Figure 8b shows a qualitative comparison.

Cosker et al.

Finally, we compare to the work of Cosker et al. [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF]. Comparisons are done over 3 sequences of D3DFACS that were provided by the authors. We obtain a mean error of 0.13mm for our method and 0.18 mm for [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF]. Figure 7c shows the results in terms of median pervertex error. They demonstrate similar accuracy although our method is more general as it does not require a controlled capture setup. Figure 8c shows a qualitative comparison.

Comparisons on efficiency

The efficiency of our method comes from both the regression step and the spatiotemporal model optimization. The regression is essential to get a good starting point that is already close to a local minimum, and it can be performed efficiently on the GPU thanks to the heightmap representation. Furthermore, thanks to the spatiotemporal model we need to optimize for much less parameters, while still remaining in a global sequence formulation. In particular, on each iteration we need to solve for the matrix C which is of size d s × d t , with d s = d exp . In our implementation d t = F/5 and thus we optimize for d exp F/5 parameters, reducing by a factor of 5 compared to a frame-by-frame formulation. The method of [7] optimizes for the parameters of each frame, which amounts to d id + F d exp variables to be solved. The main data term on the method of [20] optimizes for shape and expression parameters plus per-joint pose parameters of an articulated model on a frame-by-frame basis, increasing the complexity. Moreover, each frame is initialized from the previous one and thus it cannot be parallelized. As for UV-based methods such as [START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF], the computational complexity depends on the number of pixels of the image; while these are usually more efficient than 3D-based ones, we have shown that we can achieve similar accuracy, while remaining more general with respect to the acquisition setup.

Limitations

The use of regression for initialization allows the method to be robust to noise in the input data, but it comes with drawbacks. In particular, the use of a heightmap implies that the method is not rotation-invariant, and thus a proper pre-processing is needed to ensure that the face is "looking front". Although this does not require accurate pose detec- tion (the network was trained with data showing ±30 • of pose variation), the output will be more accurate the closer the input is to a frontal pose, and this in turn affects the final result. Furthermore, our choice of initialization can sometimes be a source of failure, particularly with the nose tip detection; see e.g. Figure 9. When this step fails all subsequent steps fail too, as regressions are inaccurate and initial correspondences cannot be found. In our experiments, this resulted in failure of some of the frames in BU-4DFE and BP4D datasets. On the other hand, dealing with motion data means that isolated failed frames can be ignored without resulting in failure of the entire sequence. In our implementation we automatically detect failed frames after ICP diverges, and this is fixed by interpolating pose and shape parameters using correct neighbouring frames. With this simple approach all sequences from BU-4DFE and 95% from BP4D were correctly registered (no errors were found during registration of D3DFACS).

Another intrinsic limitation comes from the limited scope of our trained model. Particularly for expressions that are far from this scope, our framework will provide only a coarse approximation and even the final refinement step can fail to compensate. Finally, the approach does not yet perform explicit tracking of subtle movements such as eyelids or the mouth during speech, which we leave for future work.

Conclusions

We presented a method that automatically registers large sets of dynamic 3D face scans. The main novelty is to combine a robust and efficient static registration approach with a spatiotemporal model, allowing to exploit redundancies over space and time and enabling compact representations for temporal face sequences. Remarkably, we successfully registered over three standard datasets using our method, without loosing accuracy and with significantly better time performances than related works, demonstrating therefore the efficiency of our model. In the future, we plan to use this method to build versatile statistical face models that cover large ranges of identities and expressions.

Figure 1 :

 1 Figure 1: Registration pipeline.

Figure 2 :

 2 Figure 2: Registration examples on (from top to bottom): D3DFACS, BU-4DFE and BP4D datasets.

Figure 3 :

 3 Figure 3: Results for the successive approach steps (left to right): raw scan, regression, 1 iteration, 5 iterations, final.

It. 1

 1 It. 3 It. 5 It. 10 Li et al. [20] Mean error 2.92 2.77 2.69 2.65 3.13

Table 1 :

 1 Mean landmark error in mm., for 1, 3, 5 and 10 iterations, and comparison with Li et al. [20], over 10 selected sequences of the D3DFACS dataset.

Figure 4 :

 4 Figure 4: Median per-vertex error for each frame of a long sequence in BP4D.

Figure 5 :

 5 Figure 5: Cumulative plot of median per-vertex error over D3DFACS (46028 frames) for: regression results, spatiotemporal registration (1, 3, 5 and 10 iterations) and final refinement.

  (a) Median per-vertex error (b) Mean compactness

Figure 6 :

 6 Figure 6: Comparison between our method and a static version, in terms of vertex error and sequence compactness.

  (a) [7] on 49757 frames (b) [20] on 46028 frames (c) [11] on 292 frames

Figure 7 :

 7 Figure 7: Comparisons to Bolkart and Wuhrer [7] on BU-4DFE, Li et al. [20] on D3DFACS, and Cosker et al. [11] on a subset of D3DFACS. Cumulative plots for median pervertex error.

  (a) Bolkart and Wuhrer [7] (b) Li et al. [20](c) Cosker et al.[START_REF] Cosker | A facs valid 3d dynamic action unit database with applications to 3d dynamic morphable facial modeling[END_REF]

Figure 8 :

 8 Figure 8: Qualitative comparisons. From left to right: original scan, compared method, our result.

Figure 9 :

 9 Figure 9: Failure example. (a) Heightmap obtained after bad nose tip detection (left) and the following frame (right); (b) Regression results; (c) Recovery by interpolation.

Acknowledgements We would like to thank the authors of [20],

[7] and [11] for providing comparison data.


Queensley C Chukwudum
email: queensleyv@yahoo.com

Extreme Value Theory and Copulas: Reinsurance in the Presence of Dependent Risks

Keywords: Tail dependent risks, Reinsurance treaties, Copulas, Economic capital, Stochastic simulation, Extreme value

An insurer's ability to accurately estimate the accumulation of risk, particularly in the right hand tail, is vital in ensuring that his risk appetites matches his risk exposures. This paper, therefore, focuses on the modeling of the extremal dependence structure between insurance risks using the Generalized Pareto distribution and the copula technique. The results obtained after comparing the dependence between large losses from two lines of business (motor and fire) of the Nigerian insurance industry and two specific non-life insurance companies, indicates that the correlation coefficients vary and is generally weak. With the aid of the archimedean copula, the analysis makes use of the data pair exhibiting the highest correlation to draw particular attention to the importance of taking into account the extremal dependence structure when quantifying the risk capital, allocating risk and when estimating the net reinsurance premium under different reinsurance strategies.

Introduction

When we consider more than one extremal event of interest, it is most often not a true reflection of reality to assume that the events are independent. Risk dependency can be seen between individual risks within each insurance portfolio and across lines of business (LoB) (Gildas, 2016). According to [START_REF] Li | Modeling Dependency between different Lines of Business with Copulas[END_REF], there are two broad sources of dependency between LoB -inherent dependency and other dependency. The inherent dependency is driven by factors such as inflation, interest rate, exchange rate, economic cycles and weather patterns. The other dependency arises as a result of the practices implemented by the insurer across its various LoB. Similar reserving methods and claims management procedures applied to different LoB are a few examples of such practices. [START_REF] Isaacs | Capital Adequacy and Dependence[END_REF] emphasizes the need for the insurer to properly account for dependence between risks stating that it can lead to inappropraite risk strategies as the insurer may be lulled to believe that he is less risky than he actually is.

In this article, the tail dependence structure of large losses from the motor and fire insurance portfolios are evaluated. Using the copula model, the aggregation of risks which is an important factor in the determination of the capital requirement for both the insurer and the reinsurer, is addressed. The reason for choosing the motor and fire portfolios stem from the fact that in Nigeria, cases exist where truck tanks carrying flammable products, have collided with other cars leading to multiple motor accidents and fire explosions. A recent example occured in Lagos at the Odetola bridge where over 50 vehicles were burnt including lives lost due to a fuel tanker that had rammed into another car, skidded off the road, emptied its content on the ground and exploded [START_REF] Ibrahim | Lagos Tanker Explosion: Smell of Burnt Human Flesh filled the Air -Reporter's Account[END_REF]. Another event that could lead to motor and fire risks being dependent is the act of terrorism especially in crisis-ridden areas. There are situations were vehicles are shot at, some are battered and some are set ablaze. All these lead to car crashes, people injuries and burns. Thus, there is the need to evaluate and model the tail dependence structure of these losses in order to obtain more accurate inferences.

Using copulas, [START_REF] Frees | Understanding Relationships using Copulas[END_REF] studied the dependence between liability insurance losses and allocated loss adjustment expenses (ALAE). Varying the insurer's retention at different values, they computed the reinsurance premium. [START_REF] Depuis | Multivariate Extreme Value Theory and its Usefulness in Understanding Risk[END_REF] used the same method and data to model the dependence of the joint upper tail of the distribution. Simulated hurricane losses in two U.S. states were investigated by [START_REF] Venter | Tails of Copulas[END_REF]. [START_REF] Belguise | Tempêtes: Étude des dépendances entre les branches Automobile et Incendie a`l'aide de la théorie des copules[END_REF] also made use of the copula technique to model the tail dependence between automobile and fire insurance losses. Other authors that have modelled the dependence between risk factors include [START_REF] Constantinescu | Archimedean copulas in Finite and Infinite Dimensionswith Application to Ruin Problems[END_REF][START_REF] Denuit | Actuarial Theory for Dependent Risks: Measures, Orders and Models[END_REF] and [START_REF] Tang | Economic Capital and the Aggregation of Risks using Copulas[END_REF].

To the best of the author's knowledge, there has been no study making use of African data that evaluates the extremal dependence structure of multiple lines of insurance claims and assesses its impact on reinsurance premium. This is the first research work done in this regard.

The structure of the paper is as follows -important copula concepts are presented along with the basic theory underlining the generalized Pareto distribution (GPD) in section 2. In section 3, the strength of the dependence for the three data pairs of losses (motor and fire) are analyzed. Illustrating using the motor and fire loss pair exhibiting the highest correlation, section 4 studies the impact of tail dependence on risk capital quantification and risk allocation. With the aid of stochastic simulation, section 5 focuses on estimating the net reinsurance premium under different reinsurance treaties in the presence of dependent risks and in each case, compares the results with those obtained when the risks are assumed to be independent. A summary of the analysis is presented in section 6 together with pathways to related future ideas.

Extremes and Extremal Dependence

This section reviews some basic concepts underlying the GPD and copulas. Tail dependence and some of the methods that can be used to estimate it are discussed. Then, the section is concluded with a description of the selected techniques for the application.

The Generalized Pareto distribution

This is one of the more common techniques used in the EVT field for correctly choosing the data that lies in the extremes. [START_REF] Pickands | Statistical Inference using Extreme Order Statistics[END_REF] developed the method. Suppose X 1 , X 2 , ... is a sequence of independent random variables with common distribution function F and

M n = max{X 1 , X 2 , ..., X n } Then, IP[M n ≤ z] ≤ z → G(z)
as n → ∞

given that X = {X i } n i=1 satisfies the generalized extreme value conditions.

G is a non-degenerate distribution function given as

G(z) = exp -1 + ξ z -µ σ -1 ξ +
for some µ, σ > 0 and ξ. For a large enough u the distribution function of X conditional on X > u can be approximated by

IP[X ≤ y|X > u] → H(y) as u → ∞
where

H(y) = 1 -1 + ξ y -u σ u -1 ξ + y > u (2.1)
H(y) is the GPD with the modified scale parameter σ u = σ + ξ(u -µ) corresponding to the excess of the threshold u. z + = max(z, 0). For X > u we can set X = X -u and thus rewrite the equation as

IP[X ≤ y|X > u] → H (u) (y) as u → ∞
where

H (u) (y) = 1 -1 + ξ y σ u -1 ξ + y > u
For all cases, the GPD is given as

IP(X < y|X > u) =      1 -1 + ξ y-u σu -1 ξ + ξ = 0 1 -exp -y-u σu + ξ = 0 (2.2)

The Concept of Tail Dependence

Generally, dependence between two random variables X and Y , imply that there is some kind of relationship between the variables. This implies that when we have information on X we can say something about the Y variable. Specifically, tail dependence focuses on the likelihood of one extreme event or a series of extreme events giving rise to another extreme event. Hence, according to [START_REF] Shaw | Measurement and Modelling of Dependencies in Economic Capital[END_REF], we are able to model the dependency between risks conditional on the underlying values of the risks themselves. Given that dependencies create riskier worlds, it is important to understand how the changes in an economy affect the dependency structure of risks and by extension, the economic capital of a firm. This will help in the evaluation of the amount of loss we may experience under very bad economic conditions versus under normal conditions.

Tail Dependence Coefficient (TDC)

This is a copula-based parameter developed by Ledford and Tawn (1996) which is defined as the probability that one risk variable X i takes an extreme value assuming that another risk variable X j also takes on an extreme value at the same probability level. Given the two dimensional random vector (X i , X j ) with marginal distribution functions F i and F j , then the coefficient of upper tail dependence λ u and lower tail dependence λ l of (X i , X j ) respectively are:

λ u = lim u↓1 IP(X i > F -1 i (u)|X j > F -1 j (u)) λ l = lim u↓0 IP(X i < F -1 i (u)|X j < F -1 j (u)) (3)
provided that the limit exists in both cases. F -1 i (u) and F -1 j (u) are the quantiles of risks X i and X j at level u. If C is the copula of (X i , X j ) and C is the survival copula, we can represent λ u and λ l in terms of the copula function.

λ u = lim u↓1 C(u, u) 1 -u = lim u↓1 1 -2u + C(u, u) 1 -u λ l = lim u↓1 C(u, u) u (4)
Therefore, in general, we would have to use copulas with a higher tail dependence parameter if the variables are a lot more dependent than in a normal situation.

Measures of Tail Dependence

Traditionally, the linear (Pearson) correlation coefficient is used to describe the dependence between risks. [START_REF] Shaw | Measurement and Modelling of Dependencies in Economic Capital[END_REF] argue that this measure falls short when we want to measure tail dependence because it does not represent all the possible relationships between risks. Moreover there is the need to understand the nature of the dependence structure. Other drawbacks of the Pearson correlation coefficient include the fact that it measures linear dependence and the assumption that must be satisfied by the random variables i.e. normally distributed. The Spearman rho and Kendall tau measures are better than the Pearson. More accurate measures have been developed to model tail dependence and compute the tail dependence coefficient.

Copulas

A copula C is a joint distribution function of standard uniform random variables

C(u 1 , ..., u n ) = IP{U 1 ≤ u 1 , ..., U d ≤ u n } (5)
where U i = U (0, 1) for i = 1, ..., n. Copulas allow one to characterize the dependence structure of a set of random variables separately from the marginal distributions.

Sklars theorem: This theorem forms the foundation for the use and application of copulas. Let F be the joint distribution function of the random variables X 1 , ..., X n with continuous marginal distribution functions F 1 (x 1 ), .., F n (x n ). Then, there exists a unique copula function C given by

C(u 1 , ..., u n ) = F (x 1 , ..., x n ) = F (F -1 1 (u 1 ), ..., F -1 n (u n ))
The quantile function

F -1 i is defined by F -1 i (u) = inf {x : F i (x) ≥ u}.
Furthermore, the survival copula denoted as C is defined as,

IP{X 1 > x 1 , ..., X n > x n } = C(IP{X 1 > x 1 }, ..., IP{X n > x n })
.

In this study we will focus on the bivariate case. Thus for two random variables X 1 and X 2 , the copula is a parameterization of the joint cumulative distribution function

F (X 1 , X 2 ) = IP(X 1 ≤ x 1 , X 2 ≤ x 2 ) = C(F 1 (x 1 ), F 2 (x 2 ))
and the survival function can be expressed as

C(u 1 , u 2 ) = u 1 + u 2 -1 + C(1 -u 1 , 1 -u 2 )
It should be noted that C is also a copula. It can also be used to model the dependence structure of two or more random variables. In this case it is termed the 'associated' copula. In general, though copulas are used to model the dependence structure of the whole distribution, it is also used to model tail dependence. The Archimedean and EV copulas are useful in this regard.

Families of copulas i. Elliptical copulas

This family consists of the bivariate normal and bivariate t copulas. The formula for the bivariate normal copula is

C(u 1 , u 2 ) = Φ ρ (Φ -1 (u 1 ), Φ -1 (u 2 )) (6)
Where Φ ρ is the standard normal distribution function for the bivariate case, having linear correlation ρ. Φ -1 is the inverse of the standard normal distribution function.

ii. The Archimedean copulas This family of copulas is very useful when modeling tail dependence on one side, for instance, the upper tail dependence. The general representation is given by

C(u 1 , u 2 ) = ϕ -1 (ϕ(u 1 ) + ϕ(u 2 )) (7)
The generator, ϕ : (0, 1] → [0, ∞) is a strictly decreasing convex function with ϕ(1) = 0. We assume that ϕ -1 (t) = 0 for all t ≥ lim u↓0 ϕ(u). Some well-known copulas in this family are the Clayton, Gumbel and Frank copulas.

iii. Extreme value copulas

For an independent and identically disributed (iid) sequence of random variables, the bivariate extreme value (BEV) copula come into play when we focus on the limit distribution of linearly normalized component wise maxima. It is represented by the formula:

C(u 1 , u 2 ) = exp -V -1 log u 1 , -1 log u 2 (8)
where

V (x, y) = 1 0 max w x , 1 -w y 2dH(w)
H is a distribution function with support in [0, 1] and mean 1 2 . The main limitation with the BEV copulas is that it cannot capture dependence structures that exhibit dependence at extreme values but independence in the limit [START_REF] Depuis | Multivariate Extreme Value Theory and its Usefulness in Understanding Risk[END_REF].

For an in-depth review on copulas see [START_REF] Frees | Understanding Relationships using Copulas[END_REF] and Embrechts et.al. (2003) Archimedean copulas for insurance claims

The most common Archimedean copulas are employed to compute the extremal dependence between two classes of insurance claims-auto and fire. That is, it is assumed that the portfolio is made up of these two classes of risks. The Archimedean copulas are used because we are particularly interested in understanding the dependence structure of two random variables in one region i.e. the joint upper tail. We will compare the results obtained with that from theindependence assumption to highlight the importance of incorporating tail dependence when modeling dependent insurance losses.

Clayton copula:

C(u 1 , u 2 ) = u -θ 1 + u -θ 2 -1 -1 θ θ > 0 and the generator is ϕ(t) = t -θ -1
Gumbel copula:

C(u 1 , u 2 ) = exp -(-lnu 1 ) θ + (-lnu 2 ) θ 1 θ
θ ≥ 1 and the generator is ϕ(t) = (-logt) θ

Frank copula:

C(u 1 , u 2 ) = 1 θ log 1 + (e θu 1 -1)(e θu 2 -1) e θ -1 -∞ < θ < ∞ and the generator is ϕ(t) = log e θt -1 e θ -1
In this paper, only the Clayton and Gumbel copulas are used in the example.

Aggregation of risks

Two random variables X i ∼ F i i = 1, 2 are considered. F i = F , the marginal distribution function, is known. In this case, it is the GPD. X i represents seperate classes of business that are correlated. If we set X 1 as X and X 2 as Y the the aggregate homogenous portfolio is given as

S = X + Y
The joint distribution of (X, Y ) is also specified. Thus we can uniquely determine the risk measures of the aggregate portfolio.

Reinsurance premiums

The XL reinsurance, stop loss (SL) reinsurance and a mix of both reinsurances are considered in this study. The XL reinsurance is represented as

XL = IE n i=1 [(X i + Y i ) -r]
while stop loss is given as

SL = IE [ n i=1 (X i + Y i ) -r]
where r > 0 is the retention level. X i and Y i represent motor and fire losses respectively.

Risk measures

To ensure that a firm continues operating under a worst-case scenario, the economic capital is the amount of capital that the firm must hold. A common way to determine this is by estimating the VaR at a remote probability. Other risk metrics that can be used to determine the economic capital include tail VaR (TVaR), excess TVaR (XTVaR) and expected policyholder (EPD).

T V aR = IE[X|X ≥ VaR] = ∞ p xf (x)dx 1 -p
Closely-related to TVaR is XTVaR given as TVaR minus the mean. In this case we assume that the mean is already financed by some other funding, thus we only need capital for the losses that exceed the mean.

Assuming the capital is set at VaR, EPD is then defined as the expected value of defaulted losses.

That is,

EPD = (TVaR -VaR)(1 -p)
EPD estimate is of importance to the policyholder because the policyholder is expected to bear the amount of loss in the situation where the insurer defaults on payment.

Bivariate Analysis

Data is obtained from Nigerian Insurance Commission (NAICOM 2013(NAICOM , 2014(NAICOM , 2015(NAICOM and 2016) ) available at www.naicom.gov.ng. The large claims (3.5 million naira and above) are considered for years 2013 to 2016. First, the correlation between the individual losses for the industry-wide data and data from the Custodian and American International Insurance Company Limited (AIICO) are evaluated. Then, the pair with the highest correlation is used for illustrative purposes.

Nonlinear correlations

The dependence between motor and fire losses for the Nigerian insurance industry is positive but lower than that of AIICO, while that of the Custodian insurer is negative. The dependence structure varies from negative to positive correlations although the dependence is generally weak.

Illustrative example

Amongst the selected insurers, the losses from AIICO exhibits the highest correlation (table 1). Thus it is used to illustrate the impact of dependency on some quantities that are of interest to the insurer.

The four year data is made up of 132 losses for each risk.

The rank transform plot (figure 1) displays a distinct lack of data points in the upper-left (especially)

and lower-right corners of the plot, although a general weak dependence can be observed. This indicates that the probability of a small fire claim accompanying a large motor claim and vice versa is low. In other words , the likelihood is higher for large losses or quantiles of motor to follow large losses or quantiles of fire and vice versa.

From the mep (figure 2), it can be seen that the same threshold of 4e6 is an appropraite threshold for Two bivariate distributions are created using the GPD parameters (table 3) and each of the selected copulas, then data is simulated from the distributions. For both the MLE and IFM rows in each case, the first line involves the Clayton copula while the second line is when the Gumbel is applied. The results indicate that the choice of the dependency structure plays an important part in the estimation of the bivariate distribution parameters. Though the values are similar under the MLE, there is a noticeable difference in the case of the IFM method.

Quantifying Risk Capital and Allocating Risk

The first six (of a thousand) entries of the simulated data used in each case (dependent and independent) to quantify the risk capital, allocate risk and compute the reinsurance premium is shown in tables 5 and 6. Aggregate losses in each line follows a GPD with the parameters given by the MLE estimated parameters. For the dependent case, the Clayton and Gumbel copulas are included in the model while the independence assumption is replicated by simply removing the copula dependency structure from the model. The density plots for the both the independent and dependent (Clayton) cases look very much alike except at the tails. A noticeable difference can be seen at the points where the arrow is ponting at. The density function at this point is smoother with a slight bulge for the dependent case than that of the independent case which exhibits a slightly roughened edge. This is clearly depicted in the histograms below. 145,111,384 474,897,549 595,986,491 The VaR, TVaR, XTVaR and EPD risk measures are evaluated (table 7). It is observed that when the tail dependence structure is quantified using the Gumbel copula, the risk measure estimates are all higher than in the case of the independence assumption . For the case of the Clayton copula, it is also generally higher for all other risk measures with the exception of the EPD risk is lower for for extreme quantiles when compared to the case of independence. In literature it has been shown that dependence among LoB increases the risk of the aggregate portfolio (see [START_REF] Cerchiara | Aggregate Loss Distribution and Dependence: Composite Models, Copula Functions and Fast Fourier Transform for the Danish Fire Insurance Data[END_REF][START_REF] Arachi | Solvency Capital for Insurance Company: Modeling Dependence using Copula[END_REF]. This is confirmed by the Clayton (V aR 0.97 being an exception) and Gumbel copulas. Overall, less capital will be reserved when dependence is disregarded and this can lead to unanticipated solvency problems for the insurer. Another observation is that the VaR of the aggregate portfolio differs (is lower) from that of the sum of the individual losses and the difference widens at higher quantiles (table 8).

As an example, TVaR is allocated among LoB, given that it is a more useful risk measure than VaR (see Embrechts et al. 2014). The proportion of risk allocated to motor is computed using the simple formula T V aR α (motor) T V aR α (motor) + T V aR α (f ire)

Fire is calculated in a similar way or we simply subtract allocation for motor loss from one. Allocating risk to different LoB is very useful to the insurer as it helps him to be able to set capacity controls, for example, premium targets and limits for each business unit. It also enables him to allocate capital that is in proportion to the specified risk measure. If independence between the loss variables is assumed when losses are actually correlated, the insurer will be led by this wrong assumption to allocate much more capital to motor claims instead of fire claims since the result indicates that the motor claims carry more weight when independence is assumed (table 9). In addition, an increasing trend can be noticed under the independence hypothesis while a decreasing trend is observed when tail dependence is accounted for. The choice of copula also affects the value of the estimate of the allocated risk.

Premium Valuation

In this section, 4 possible reinsurance treaties are considered under both the dependence and independence assumptions. They are as follows:

Option 1 No reinsurance treaty.

Option 2 Individual XL for each loss exceeding the specified layers.

Option 3 Mixed reinsurance treaty. Two cases are considered.

Case 1 -motor has proportional treaty covering 70% of the specified layers while fire has an aggregate XL treaty covering same layers.

Case 2 -fire has proportional treaty covering 70% of the specified layers while motor has an aggregate XL treaty covering same layers.

Option 4 Stop loss reinsurance covering the whole portfolio at the specified layers.

Only limited covers are considered for the three different layers in options 2 to 4. The layers represent millions of naira. The stability of the result is quantified using standard error (SE) which can be seen in the braces. 43,285,420 (187,224,002) The impact of the choice of copula is made visible from the results displayed on table 11. In general, when the independent case is compared to dependent case, the difference between the reinsurance premium estimates in each case is not much. Additionally, if the dependence structure is specified using the Gumbel copula the reinsurance net premium is higher than in the clayton example. The SE also follows the same trend. For option 3, case 2 reports a lower net premium and SE than case 1.

The reinsurance net premiums across the options indicate that option 3 case 2 gives the insurer the opportunity to pay the lowest premium along with the lowest SE while the highest premium to be paid occurs when option 2 is chosen. In a situation where the insurer decides not to take up reinsurance, the uncertainty associated with the mean aggregate claim shoots up (table 10).

A suprising find is observed under the mixed treaties and stop loss reinsurance (both independent and dependent assumptions). Here, the layer (30,50] is generally lower than the (50,100] layer with a much higher SE. This does not fall in line with the known fact which is, the higher the retention layer, the lower the reinsurance premium. This fact is only broadly confirmed under the individual XL reinsurance structure. [START_REF] Eling | Risk and Return of Reinsurance Contracts under copula Models[END_REF] performed a similar analysis where they integrated several copula models and reinsurance contracts. Their simulation study revealed large differences in risk assessment under different copulas and reinsurance treaties. A similar situation is noted in this paper. Although the introduction of nonlinear dependencies led to a higher number of large losses in their analysis for all copulas when compared to the independent case, the situation differs a bit in this study given that both higher and lower losses are experienced when different copulas are employed. A point of note is that they focus on the expected gain of the insurer under the different reinsurance contracts while the focus here is on the expected losses above specified retention levels. Dhane and Goovearts (1996) showed that stop loss premium is higher when dependence is accounted for than under the independence assumption. In this study, the situation depends on the selected dependence structure.

Gildas (2016) studied loss ALAE claims and generally noted that the portfolio is less risky when independence is assumed. This find is confirmed true for a portfolio incorporating option 1 in this research (table 10). An interesting point raised in his study was that as the retention level is increased, the reinsurance premium for the XL treaty under different copula models tend to the estimates obtained under the independent assumption and that of the stop loss reinsurance deviates more from the independence results as deductible increases. In contrast, the result here shows that the reinsurance premiums under the copula models are much closer to the independent estimates when the lowest retention level is considered (table 11).

Conclusion

A homogenous portfolio consisting of the aggregate bivariate losses from the Nigerian insurance sector was analyzed. The techniques employed were GPD and copulas. This paper showcased the impact of tail dependence on risk measures, risk allocation and reinsurance net premium. When losses are correlated, the value of the risk measure is increased and risk allocation is also affected. The choice of copulas also plays a major role as we can see reinsurance premium generally being lower than (in the case of the Clayton copula) and higher (in the case of the Gumbel copula) than the reinsurance premium obtained under the independent hypothesis.

A related problem would be to look at the inhomogeneous case where the F i 's are not necessarily the same. Another subject of investigation is dependence uncertainty which is a situation whereby the marginal distributions are fixed and known but the joint model is unspecified.

Figure 1 :

 1 Figure 1: Scatterplot of AIICO log claims (left) and rank transformed data(right)

Figure 2 :

 2 Figure 2: Left: Mean excess plot for motor (up) and fire (down). Center and right respectively: gof tests for motor and fire.

Figure 3 :

 3 Figure 3: Density function of the aggregate losses: independent (left) and dependent (right) cases

Figure 4 :

 4 Figure 4: Histograms of the aggregate losses: independent (left) and dependent (right) cases

Table 1

 1

	: Kendall tau and Spearman rho
	Motor and fire claims kendall tau Spearman rho
	Industry-wide data	0.044	0.0662
	Custodian	-0.054	-0.0838
	AIICO	0.124	0.122

Table 2 :

 2 Descriptive statistics for AIICO's individual losses

	Industry data	motor	fire
	Min.	3247500	42601
	1st Qu	4050000	5050086
	Median	4631500	7957982
	Mean	5599483	18099085
	3rd Qu.	5913000	15465777
	Max. 51810400 345156235
	AIICO
	Min.	3247500	3509775
	1st Qu.	4041112	5801860
	Median	4569500	9674585
	Mean	5206237	23048856
	3rd Qu.	5636156	27801345
	Max. 14073345 152991576

Table 3

 3

	: GPD parameters
	scale shape
	motor 1.158340e6 0.318
	fire 8.307612e6 0.759

Table 5

 5

		: Simulated data for the independent case
	sim number motor loss	fire loss aggregate portfolio loss
	1	5771942 22318941	28090883
	2	4286236	9309488	13595724
	3	5198114 36983903	42182017
	4	4147002	5399101	9546103
	5	4071796	5461946	9533742
	6	10088735 12223174	22311909

Table 6 :

 6 Simulated data for the dependent case (Clayton)

	sim number motor loss	fire loss aggregate portfolio loss
	1	4414760	6749374	11164134
	2	6326208 17286991	23613199
	3	4663138	5707404	10370542
	4	7565872 47799368	55365241
	5	9294422 46886811	56181233
	6	4054413	7178607	11233020

Table 7 :

 7 Risk capital for total risk

	Indep	97%	99%	99.5%
	VaR 156,678,722	291,569,402	575,659,669
	TVaR 377,052,095	706,533,031	959,021,998
	XTVaR 340,661,947	670,142,883	922,631,851
	EPD	6,611,201	4,149,636	1,916,812
	Dep (Clayton)
	VaR 153,178,320	483,822,375	605,701,183
	TVaR 448,168,109	772,008,670	984,129,772
	XTVaR 410,753,821	734,594,381	946,715,483
	EPD	8,849,694	2,881,863	1,892,143
	Dep (Gumbel)
	VaR 167,138,545	349,112,469	633,809,791
	TVaR 588,770,049 1,272,882,384 2,122,309,339
	XTVaR 545,484,629 1,229,596,964 2,079,023,919
	EPD	12,648,945	9,237,699	7,442,498
	Indep means independence, Dep means Dependence.
	Table 8: Risk capital for individual risk
		97%	99%	99.5%
	Indep
	VaR(motor)	11,071,264	16,566,065	21,289,669
	VaR(fire ) 152,545,832 286,563,060 571,014,146
	Dep (Clayton)
	VaR(motor)	11,362,522	15,698,533	17,976,197
	VaR(fire)

Table 9 :

 9

	Risk allocation for motor loss
	97%	99% 99.5%
	Indep 0.466 0.547 0.560
	Dep (Clayton) 0.220 0.230 0.174
	Dep (Gumbel) 0.062 0.035 0.023

Table 10

 10

	: No reinsurance (option 1)
		Option 1
	Indep	36,390,148 (87,405,034)
	Dep (Clayton) 37,414,289 (92,627,942)
	Dep (Gumbel)

Author: Lecturer at the University of Uyo, Department of Insurance, Nigeria. Currently a doctoral student in financial mathematics at the Pan African University of Basic Sciences, Technology and Innovation Kenya.


PhD Queensley C Chukwudum
email: queensleyv@yahoo.com

Credibility Premium Estimation of Insurance Claims in Nigeria

Keywords: Mathematics Subject Classification 2010: 62C12 EBCT model 1 (Buhlmann Model), Credibility Premium Insurance Portfolios, Unbiased Estimators, Aggregate Net Claims, Chisquare goodness of fit test

In this paper one of the foundational theories in actuarial science -the Empirical Bayes credibility theory model 1 (Buhlmann model), is applied to estimate the pure risk premium based on the aggregate net claim data of 5 years (2009)(2010)(2011)(2012)(2013). This involves estimating the credibility premium in the next year (2014). Five different classes of risk are considered from the insurance industry in Nigeria -fire, accident, motor, marine and, oil and gas insurance portfolios. Verification of the unbiasedness of the structural parameters of the credibility factor is presented. A chi-square goodness of fit test is carried out to compare the estimated values with the actual claim data. The result shows that there is no significant difference between the estimated values and the actual values for the year 2014.

INTRODUCTION

One of the main branches of insurance mathematics is credibility theory. Intuitively, credibility theory generally deals with the measure of credence or 'trust' that should be assigned to a particular body of experience for the purposes of rate making. For example, if the loss experience in a given class of risk is termed as 'not credible enough', it means that we cannot rely on such data to forecast future events or costs because it will most probably produce incorrect results. Therefore being able to measure how reliable a given experience is, is very vital to the survival of any firm. It is an extremely useful tool in risk management practices. Furthermore, this study will aid in enlightening policyholders on issues regarding insurance claims.

Section 2 delves into the foundational principles of the EBCT model whereas section 3 deals with the data and methodology involved. In section 4, the data is analyzed and results are discussed. Conclusions are made in section 5.

THEORETICAL BACKGROUND

The advent of credibility dates back to Whitney, who in 1918 addressed the problem of assessing the risk premium m, defined as the expected claims expenses per unit of risk exposed, for an individual risk selected from a portfolio of similar risks. Advocating the combined use of individual risk experience and class risk experience, he proposed that the premium rate be a weighted average of the form

= z + (1 -z) μ
where is the observed mean claim amount per unit of risk exposed for the individual contract and μ is the corresponding overall mean in the insurance portfolio, the weight z was soon to be named credibility (factor) since it measures the amount of credence attached to the individual experience, and was called the credibility premium [START_REF] Norberg | Credibility Theory[END_REF]. Thus credibility theory is all about weighted averages. Different estimates of a quantity are to be weighted together. The more credible estimates get more weight. It is usually presented as a mathematically dense body of formulas [START_REF] Venter | Credibility Theory for Dummies 3[END_REF]. In a simpler way, [START_REF] Mahler | Chapter 8: Credibility[END_REF] describes it with a clear example. Rather than relying solely on recent observations, better estimates may be obtained by combining this data with other information. For example, suppose that recent experience indicates that carpenters should be charged a rate of $5 (per $100 of payroll) for workers compensation insurance. Assume that the current rate is $10. What should the new rate be? Should it be $5, $10, or somewhere in between? Credibility is used to weight together these two estimates. The basic formula for calculating credibility weighted estimates as:

Estimate = Z [Observation] + (1 -Z) [Other Information], 0 ≤ Z ≤ 1
where: Z is the credibility assigned to the observation. 1 -Z is commonly termed as the credibility complement.

If the body of observed data is large and not likely to vary much from one period to another, then Z will be closer to one. On the other hand, if the observation consists of limited data, then Z will be closer to zero and more weight will be given to other information. Therefore, solving the problem stated above: The current rate of $10 in the above example is the "Other Information." It represents an estimate or prior hypothesis of a rate to charge in the absence of the recent experience. As recent experience becomes available, then an updated estimate combining the recent experience and the prior hypothesis can be calculated. Thus, the use of credibility involves a linear estimate of the true expectation derived as a result of a compromise between observation and prior hypothesis. The carpenters' rate for workers compensation insurance is therefore Z * $5 + (1 -Z) *$10 under this model. Many different models having varying levels of complexity have been developed to solve credibility problems. Hans Buhlmann formulated the Empirical Bayes credibility theory (EBCT) model in 1967. The simplest of the EBCT models is the EBCT model 1. It provides a basic understanding of how EBCT models work and how the results obtained can be analyzed. One major advantage of EBCT models is that the collateral data speaks for itself-it is allowed to provide the 'between risks' rather than making distributional assumptions regarding the collateral data and the prior knowledge.

In summary EBCT= Bayesian credibility theory -distributional assumptions + collateral data, (Gray R.J. et al 2012). It can be regarded as a generalization of the normal/normal model -The Actuarial Education (2010). The literature covering the application of the EBCT model is quite rich especially in the developed countries. In developing countries, Africa in particular, the studies are quite scarce. [START_REF] Mawuli | Application of Buhlmann's Credibility Theory to an Automobile Insurance Claims[END_REF] applied Buhlmann credibility theory to an automobile insurance claims data from Ghana. Garnier (2016) has applied the model to life insurance. In order to estimate net premiums for catastrophic losses by world regions (continents), [START_REF] Jindrová | Empirical Bayes Credibility Models for Economic Catastrophic Losses by Regions[END_REF] applied both models 1 and 2.

DATA AND METHOD

In this section, we provide a brief description of the data and show how the method is applied.

Data

The data set was obtained from National Insurance Commission (NAICOM). Five different risk classes are chosen. The data shows the net claims paid over a 5 year period from the following insurance portfolios within the Nigerian insurance industry: Fire, accident, motor, marine, oil and gas. The data (in millions of naira) is given in Table 1.1.

Objective

The aim of this research is to use the EBCT model 1 to analyze the above data in order to obtain the credibility premium. The experience of this collection of five risks is the collateral information. The result obtained will be used to project (compute) the expected claim payment (pure premium) for each class in the next year (2014). Furthermore, a goodness of fit test is carried out using the Chi-square goodness of fit test.

Method

Let

… denote the aggregate claims, or the number of claims, in successive periods for each risk. Thus we have different risks which we denote by N. In this case N=5. The risk data is collected over a couple of years, i.e. Y1,Y2,…Yn. The number of years is denoted by n. Similarly n=5 in this case. We need to estimate the expected value of Yn+1. To do this, the credibility factor Z, will first be computed. The formula is given as:

Z = = =
Then the credibility premium frequency is computed. The formula is given as

+

The following observations about Z are noted: (i) Z is always between zero and one.

(ii) Z is an increasing function of n. This implies that if more data from the risk itself is available, then we will rely more on it when the credibility estimate of the pure premium or number of claims is computed.

(iii) Z is a decreasing function of . Hence, if is higher in value relative to , then the data from the risk itself will vary more when compared to the data from the other risks in the collective. Hence, data from the risk itself will be less reliable. (iv) Z is an increasing function of . Hence, if the value of is higher relative to , more variability will exist between the different risks in the collective and hence it will be less likely that the other risks in the collective will resemble the risk that is of interest. Thus, less reliance should be placed on the data from these other risks.

Definition of the above quantities

i)
is the mean of each of the risks given.. , =rows (risks), = column (years), represents the aggregate claims in the th year from the th risk.

ii)

=

. This is the overall mean of all the risks combined.

iii)

. This represents the average variance of claim within the risk. The term in the parentheses on the right hand side of the equation represents the variance of each risk. iv)

=

-. This represents the variance of mean claims between (across) risks.

Assumptions for EBCT Model 1

1. The distribution of each depends on a parameter θ. The value of θ is fixed (and the same for all the ) but is unknown. 2. Given θ, the 's are independent and identically distributed (iid).

3. For different risk i ≠ j, the pairs ( ) and ( ) are iid.

The parameter θ is the risk parameter. It could be a real number or it could be a more general quantity such as a set of real numbers.

The structural parameters and their respective estimators

Table 1.2: A layout of the parameters and their estimators Structural parameter Estimator 1 2 3 -

Verification of the unbiasedness of the above estimators

One of the properties of a good estimator is that it should be unbiased. An estimator is unbiased if . This implies that the difference between the estimator's expected value and the true value of the parameter being estimated should be zero i.e.

. This is a very useful property because it tells us which estimator gives a more precise estimate.

Theorem 1 : If …, are independent with a common mean μ= and a common variance , 1≤i≤n and i is an integer .Then the sample mean and the sample variance respectively are: and . and are unbiased estimators for μ and respectively i.e. and . The former is easy to prove. The proof of the latter will help us understand 2 and 3 in table 1.2. A detailed proof is given in the appendix.

DATA ANALYSES AND DISCUSSION OF RESULTS

Each row represents observations for a different risk. The mean and variance for each risk is computed. Oil and gas displays the highest variance (Table 1.3). Further analysis shows that the variation within the risk is low compared to the variation between average claims across risks. This result is responsible for the high value of the credibility factor (Table 1.4). Given that overall variance is larger, it means that little weight can be put on the information that is provided by the data set as a whole but more weight can be put on the data from the individual risks.

The credibility premium for each risk is calculated. Table 1.5 displays the projected result for the expected claim payment (in millions of naira) under the EBCT model 1 for the year 2014. It should be emphasized that the risks are heterogeneous meaning that they have different risk parameters. Also we can observe that the credibility premium for each risk varies. This is a reflection of the claims experience of the risks. A quick inspection of the relationship between the mean claims in the available history of a risk and the credibility premium reveals that the higher the value of the mean claim of an individual risk set, the higher will be the credibility premium for that risk (Table 1.6) The actual claim payment data for the 2014 is obtained from NAICOM 2014 annual general business data as shown (Table 1.7). The largest difference between the expected and actual data exists in the motor insurance class. The plot (fig 1 .2) reveals that both the estimated and the actual (observed) data have similar trend with the estimated values being higher than the actual values. A Chi-square goodness of fit test is carried out to see how well the observed frequency matches the expected frequency where observed frequency is the actual data and expected frequency is the expected claim payment (credibility premium.. The difference between the observed and the expected frequencies has been computed in billions of naira to make the chi-square test computation easier.

Hypothesis testing

H0: O=E (i.e. there is no statistically significant difference between the observed (O) and the expected(E) frequencies)

H1: O≠E (i.e. there is a statistically significant difference between the observed (O) and the expected(E) frequencies)

Level of significance (α) =0.05

Decision rule: if p-value< α we reject H0 (i.e. the null hypothesis), if p-value>α we fail to reject H0. Diff. means difference, sq. means square.

Result:

The p-value = 0.28 > α =0.05, hence we fail to reject H0.

Report: There is no significant difference between the observed and expected frequencies at α =0.05.

This implies that the estimated value is acceptable.

Discussion of results

The estimated values for the expected claim payment are higher than the observed claim payments for the year 2014 in each class of risk. Notwithstanding this result, it can be said that it is a much safer estimate given to the fact that the result gives the insurer an idea of the amount to set aside as reserve for the payment of claims. It is better to set aside more reserves than less, so as to cover all claims that will be made in the coming year. Despite its shortcomings which includes the fact that the adjustments for the volume of business in each year is not taken into account and the assumptions are quite restrictive, the EBCT model 1 still gives an adequate estimate of the credibility premium for each class of risk.

CONCLUSION

This research showed how the EBCT model 1 can be used in estimating insurance premiums based on the data obtained from the Nigerian insurance sector. To improve on the work done, further research can be undertaken using the EBCT Model 2 which allows for changing levels of business. This reflects what happens in real life. Data for the volume of business in the insurance industry (on the given classes of risks), can be incorporated into the estimation process. This will give rise to better estimates and more approximate credibility premium results.

ACKNOWLEDGEMENT

The author is a lecturer at the University of Uyo, department of Insurance, Akwa Ibom State, Nigeria. This paper was presented at the 2018 African Finance Journal conference held in Kenya.

APPENDIX

Proof of theorem 1.

Double expectation: The concept of double expectation makes it possible to write as when the mean is conditioned on .

i) is an unbiased estimator for E(m(θ))

Proof:

Since we first find the mean of each risk before we find the mean across risks. We can write:

Fig

  Fig 1.1: Plot showing the trends for the 5 classes of risk.

Fig

  Fig 1.2: A graphical representation of the expected and observed claim payment.

  expectation of a sum is the same as the sum of the expectation, we therefore have: term and the second terms on the right hand side are unbiased estimators, it implies that is unbiased for .

Table 1 .1 General Insurance Net Claims (2009-2013) in millions of naira(N000,000)

 1

		FIRE	ACCIDENT	MOTOR	MARINE	OIL & GAS
	2009	6696.89	4945.07	10777.66	3951.33	8161.73
	2010	4176.38	4841.85	11889.22	2217.64	-872.75
	2011	4029.035	5643.9096	11721.66	1534.436	1456.384
	2012	10168.35	8906.28	14853.84	5204.59	11899.45
	2013	8231.68	7290.17	22342.93	4046.65	10634.25
	Source: NAICOM 2013 Industry Statistical Analysis

Table 1 .3: Results for the mean and variance of each risk.

 1

	RISK/ YEARS 2009	2010	2011	2012	2013		Variance of
								each risk
	FIRE	6696.89	4176.38	4029.03	10168.35 8231.68	6660.47	6967600.773
	ACCIDENT	4945.07	4841.85	5643.91	8906.28	7290.17	6325.46	3040595.831
	MOTOR	10777.66 11889.22	11721.66 14853.84 22342.93 14317.06 22465140.14
	MARINE	3951.33	2217.64	1534.44	5204.59	4046.65	3390.93	2214140.035
	OIL & GAS	8161.73	-872.75	1456.38	11899.45 10634.25 6255.81	32126199.77
	SUM						36949.73 66813676.54

Table 1 .4: The credibility factor

 1

	OVERALL MEAN (Ӯ)	7389.945534
	AVE. VAR OF CLAIMS WITHIN RISKS	13362735.31
	VAR OF MEAN CLAIMS ACROSS RISKS	19404702.23
	CREDIBILITY FACTOR	0.87894565

Table 1 .5: Credibility premium computation (in N000,000)

 1

	RISK	EXPECTED CLAIM PAYMENT IN 2014
	FIRE	6748.773535
	ACCIDENT	6454.317026
	MOTOR	13478.50517
	MARINE	3875.027463
	OIL & GAS	6393.10448

Table 1 .6: Relationship between the mean claim of each risk and their respective credibility premiums (in N000,000) Risk Mean claim of each risk set Credibility premium

 1

	FIRE	6660.47	6748.773535
	ACCIDENT	6325.46	6454.317026
	MOTOR	14317.06	13478.50517
	MARINE	3390.93	3875.027463
	OIL & GAS	6255.81	6393.10448

2: A graphical representation of the relationship in Table 1.6

Table 1 .7: Expected versus actual claim payment (for 2014) data EXPECTED CLAIM PAYMENT ACTUAL CLAIM PAYMENT DIFFERENCE

 1

Table 1 .8: Hypothesis test result. Diff. in billions of naira Expected claim in billions of naira Diff. sq./expected claim

 1

	1.12	6.75	0.19
	1.85	6.45	0.53
	5.93	13.48	2.61
	2.18	3.88	1.22
	1.79	6.39	0.50
		chi-sq. test statistic	5.05
		p-value of the statistic	0.28

We now make use of this standard result which generally states that the sample variance of a random sample of observations is an unbiased estimator for the population variance. By implication, if is an unbiased estimator for , then we can also say that is unbiased for meaning that, = . Furthermore, we note that, is simply the mean of . Hence is unbiased for .

iii) -is unbiased for .

Proof:

From the results obtained in (ii) above, we can also say that is unbiased for ………………………… (1*)

We then apply the law of total variance which states that = + as shown below.

Recall that . Similarly, .

Therefore applying the conditions of we can write:

………………………(1**) Given it implies that . Thus, equation 1** becomes