Log-normal approximations in the Solvency II Standard Formula

(1)

Log-normal approximations in the Solvency II Standard Formula

Peder T. Haagensli

Master’s Thesis, Spring 2017

(2)

This master’s thesis is submitted under the master’s programme Modelling and Data Analysis, with programme optionFinance, Insurance and Risk, at the Department of Mathematics, University of Oslo. The scope of the thesis is 30 credits.

The front page depicts a section of the root system of the exceptional Lie group E8, projected into the plane. Lie groups were invented by the Norwegian mathematician Sophus Lie (1842–1899) to express symmetries in differential equations and today they play a central role in various parts of mathematics.

(3)

Abstract

On the 1st of January 2016, the Solvency II Directive regulating the European insurance industry came into force. Along with a host of reporting requirements, quantitative algorithms for calculating key quantities are provided in the Solvency II documentation.

One of these key quantities is the Solvency Capital Requirement (SCR) providing an insurance company with a "soft" oor in terms of capital that it needs to hold. A vital element in calculating the SCR for an insurance company is the so-called Standard formula based on Gaussian risks which is used iteratively to aggregate the dierent risk elements in the business.

An alternative to the Standard formula was proposed by Bølviken and Guillen who introduced log-normal distributions to better capture the skewness often present in insurance risks. They used classical moment matching to approximate the distribution of a sum of log-normal risks. In this thesis we will consider another log-normal alternative making use of moment-generating functions. This method will then be shown to oer superior accuracy compared to the Standard formula and the method of Bølviken and Guillen.

Keywords: Cholesky decomposition, Gauss-Hermite quadrature, moment-generating function, principal components, Solvency II, spectral decomposition, Standard formula, sum of log-normals.

(4)

1 Introduction

Solvency II is a set of regulations from the European Insurance And Occupational Pensions Authority (EIOPA) placed upon insurance companies conducting business in the EU. The directive institutes a collection of algorithms for calculating a set of quantities which will dictate the amount capital an insurance company will be required to hold to be allowed to oer insurance. Solvency II also issues a set of protocols with respect to reporting to national authorities and improving public disclosure in the insurance business.

With its introduction EIOPA seeks to establish an EU wide license supporting the idea of a single market in the EU.

Risk aggregation in Solvency II takes a modular approach. Referring to Figure 1 which is taken from Bølviken (2017). Consider the overall risk of a modern insurance company which may be broken down into many dierent types of risks. Take for example the risk linked with market investments. This market risk may be viewed as a sum of sub-risks associated with the interest rate, equity, property and so on. On the other hand, the overall market risk may also be viewed as a part of a larger aggregated risk along with life risk, non-life risk etc.

The Solvency Capital Requirement (SCR) is a key quantity in Solvency II. The SCR represents the amount of capital an insurance company is required to hold in order to reliably be able to meet its liabilities. The 99.5 percent principle is a key feature in Solvency II which means that there should in principle be no more than a 0.5 percent chance of the insurance company going bankrupt within one year. In mathematical terms, we may think of the SCR as a percentile of an associated risk variable. An insurance company is required to calculate its top-level company SCR regularly. This is done by starting at the lowest possible level of Figure 1. Solvency II provides us with algorithms for approximating the SCRs connected to all of the input nodes at the lowest level.

Then after obtaining all the risks making up a particular node one level up we aggregate these risks into a single new SCR for that node. For example after obtaining SCRs for Health^li, Health^nli and Catastrophe these are used to aggregate up one level, giving us the Health SCR. After doing this for the other nodes Market, Life, Non-Life etc. these SCRs are used to aggregate up one level again. This continues until all of the risks making up the business have been aggregated to one single parent node. Details beyond what is presented in this thesis regarding this process can be found in EIOPA (2014) and Commission Delegated Regulations (2015).

This aggregation process uses a mathematical formula referred to as the Standard formula which is mathematically correct if the underlying risks are Gaussian. Conse- quently, the Standard Formula may seem tting when considering for example nancial risk which often exhibits approximate normality. Insurance risk however, will generally not t any Gaussian assumption. This breach of assumption may reveal that the 99.5 percent principle in Solvency is misleading.

Bølviken and Guillen (2017) proposed an extension based on log-normal risks. They used the method of matching the moments of the sum of the underlying log-normal risk variables with that of a new log-normal variable. There is however another alternative from engineering literature presented by Methe et al. (2007). This method also creates a new log-normal variable, but through matching the moment-generating function of the new variable with that of the sum of the underlying log-normal variables. The purpose of this thesis is to develop and apply this technique to the model in Bølviken and Guillen (2017). Furthermore, the performance of the resulting method will be compared to the

(7)

method in Bølviken and Guillen (2017) and the Standard formula. The procedures will be implemented using the programming language R for which the source code will be provided on page 25. A large part of the time dedicated to this thesis was spent producing the code for these methods from scratch.

Figure 1: Risk evaluations in Solvency II and its aggregation. There is also another level below for some nodes (not in gure).

2 The Solvency II approach

2.1 Introduction

We will now try to formalize the methodology described in section 1. Consider a random variable Z representing the total company net loss for the coming year. If it is negative then the business is protable that year. The expectation E[Z] is predictable and the Solvency II approach is to examine

Y =Z−E[Z].

The Solvency Capital Requirement is then the solution of the equation

P(Y >SCR) = (1)

with = 0.5%. In order to reasonably assume the company to be able to meet its obligations Solvency II require the insurance company to hold assets valued at leastE[Z]+

SCR at the start of the year.

Consider the upper layer in Figure 1 with risk variables Z =Zôper+Z^basic+Zâdjust. Subtracting the means then yields the variables of interest Y =Yôper+Y^basic+Yâdjust.

(8)

Equation (1) may be solved through the use of Monte Carlo, but that involves much modelling and is not the Solvency II way. Instead, SCRs on the lowest levels are calculated and aggregated iteratively in order to end up with the SCR for the entire company.

Consider for example the top-level SCR^company which is immediately comprised of and calculated from SCR^oper, SCR^basic and some adjustment term Adj. Now, SCR^basic is computed from the SCRs from market, lif e and so on, and those are in turn computed from the SCRs of their own underlying variables. Those may have other risk variables under them and the same technique of merging risks from their percentiles are used there.

This scheme deals with many types of risk and the true distribution of these dierent risks may vary. Financial risk will often emit Gaussian traits while risk connected with say natural catastrophes in insurance will be much more skewed. With these considerations notwithstanding, the risks are handled in the same way in regards to risk aggregation.

How is this carried out?

The general case is K risk variables Y_i, ...Y_K with mean equal to zero for which their Solvency Capital Requirements SCR1, ...SCRK are to be aggregated to the corresponding SCR for the sum

Y_S =Y₁+...+Y_K.

Of course, for the bottom nodes where the scheme starts the SCRs must be found by other methods. How this is done is not relevant to the scope of this thesis, but one may consult EIOPA (2014) to see how this is done for each node. Our objective is the aggregation of percentiles from one level in Figure 1 to another.

2.2 The Standard Formula

We still assume that we have K nodes representing the dependent risk variablesY₁, ..., Y_K making up a parent node and risk variableY_S. These random variablesY_khave expectation equal to zero and 99.5 percentiles SCR1, ...,SCRK. The dependence structure is handled through the use of correlation coecients

ρ_ij =cor(Y_i, Y_j)

which are supplied by Solvency II documentation. See EIOPA (2014) or some of the samples on page 23. We then seek the corresponding SCRS or 99.5 percentile of the sum Y_S =P

Y_k that we obtain through the Standard formula SCRS =

K

X

i=1 K

X

j=1

ρ_ij×SCRi×SCRj

!1/2

. (2)

If the risk variables Yk are Gaussian with mean zero, then (2) is mathematically exact.

To see this, observe that if we let φ be the 99.5th percentile of the standard normal distribution and denote the standard deviation of risk k = 1, ..., K as σ_k then

SCRk=φ×σ_k, k = 1, ..., K.

Now take the variance formula for the sum Y_S, which is

σ_S =

K

X

i=1 K

X

j=1

ρ_ijσ_iσ_j

!1/2

,

(9)

and multiply both sides by φ and you arrive at (2).

The Standard formula may therefore seem tting when considering nancial risk which often exhibits approximate normality. Insurance risk, however, will generally not t any Gaussian assumption. This breach of assumption may mean that the 99.5 percent principle in Solvency II is misleading. However, in spite of this, there is no denying the simplicity of (2) in terms of implementation and its ability to quickly aggregate many sources of risks quickly and easily. There is also the question of whether or not describing the dependence through the use of correlation coecients is appropriate or not for all the dierent types of risk. We will not address that question in this paper. In the next section we will explore alternative methods to the Standard formula which will abandon any Gaussian assumption of the underlying risks Y_k making up the sum. We will then compare the accuracy of these methods and eventually discuss any trade-os with these methods in regards to practical implementation.

3 A log-normal standard formula

3.1 Introduction

As noted in the previous section the mathematics behind the Standard Formula which the risk aggregation in Solvency II is built upon assumes Gaussian risks. An alternative approach is to assume the underlying risks to be log-normal. Our model (3) is taken from Bølviken and Guillen (2017) and should better accommodate non-Gaussian risks skewed to the right which is common in insurance. We will assume risks of the form

Y_k =σ_ke^−τ^k²^/2+τ^k^k−1 pe^τ^k² −1

, k = 1, ..., K (3)

where the elements in = (_k) are Gaussian with mean zero, standard deviation 1 and correlation matrix Σ. Furthermore, σk is the standard deviation of Yk and τk is a shape parameter. Given the shape parameter τ_k the skewness γ_k is given by the formula

γ_k = (e^τ^k² + 2)p

e^τ^k² −1. (4)

One may show by application of L'Hôpital's rule that Y_k would approach a normal distribution if we were to let τ_k tend to zero. Choosing the log-normal distribution over a Gaussian represents a step towards more skewed distribution which should be closer to reality when it comes to insurance risk. The log-normal distribution in particular also has some properties that makes it easy to work with. Note that Y_k has expectation zero.

Our goal is then to approximate the sum PK

k=1Y_k with a new log-normal variable Y_S =σ_Se^−τ^S²^/2+τ^S^S −1

pe^τ^S² −1

for some new parameters (σ_S, τ_S). This problem has no exact solution. However, it has been studied a great deal and several methods have been developed.

3.2 The Fenton-Wilkinson approximation

Bølviken and Guillen (2017) have already studied this particular problem through the use of the so-called Fenton-Wilkinson scheme. They considered the skewness and standard

(10)

deviation of the sum of the log-normals P

Y_k and selected the new parameters (σ_S, τ_S) in order to match the skewness of P

Y_k and Y_S. They selected σ_S through the ordinary variance-formula for a sum of correlated variables.

σ_S =

n

X

i=1 n

X

j=1

ρ_ijσ_iσ_j

!1/2

(5) The skewness coecientsγ_k is readily obtained through (4) and they showed the skewness of the sum P

Yk to be given by γ_S =X

i

α³_i(e^τⁱ² + 2)(e^τ^j² −1)²+ 3X

i6=j

α²_iα_jh

e^τⁱ²(β_ij² −1)−2(β_ij −1)i + 6 X

i<j<k

α_iα_jα_k[β_ijβ_ikβ_jk−β_ij −β_ik−β_jk+ 2] (6) where

α_i = σ_i σS

pe^τⁱ² −1 and β_ij = 1 +ρ_ij r

e^τⁱ² −1

e^τ^j² −1 . For their proof see the appendix of Bølviken and Guillen (2017)

They then used these results to match the skewness of the new random variableY_Swith PYk. Finally, with γS in hand they then obtain τS by solving equation (4) with respect to τ_k. This then fully species the new random variable Y_S such that skew(P

Y_k) = skew(Y_S).

This method may be generalized to distributions other than the log-normal. For eligible distributions, an exact expression for γ_S might not exist in which case numerical methods may be introduced in order to solve the equations. Bølviken and Guillen (2017) also demonstrate how copulas may be incorporated instead of correlation coecients to better capture dierent dependence structures. All in all this approach is a simple and exible solution to the problem and their results suggest promising accuracy for problems with K = 3 log-normal risks for dierent degrees of skewness.

3.3 Matching moment-generating functions

We will now develop another way of approximating the sum of log-normals with a new log-normals. The method will have similarities with the methods of moment-matching in Bølviken and Guillen (2017), but will rather obtain the parameters (σ_S, τ_S)by matching the moment-generating functions at some point t.¹ This concept is introduced in Metha et al. (2007), but we will also suggest a modication in an attempt to reduce the computational burden of problems with a large number of dimensions K.

We will consider the same model as Bølviken and Guillen. The properties of the log-normal distribution are well-studied. The pdf is given by

f_Y_k(y) = 1

˜

σf_V_ky

˜ σ + 1

, −˜σ_k < y <∞

1We will more specically use the characteristic function which is a special case of the moment- generating function. However, we will not make any distinction between these two and treat them as one.

(11)

where we have dened V_k := exp(−τ_k²/2 +τ_k_k) and σ˜_k := σ_k/p

e^τ^k² −1 to simplify the expression. The pdf of V_k is readily known as

f_V_k(v) = 1 vτk

√2πexp

−(lnv−(−τ_k²/2))² 2τ_k²

. The moment-generating function (MGF) is written as

M_Y_k(t) =E e^−tY

= Z ∞

−˜σk

exp(−ty)f_Y_k(y)dy. (7) We may here interpret the exponential function as a weight function with a tuning parameter t. If we then set the MGF of Y_S equal to the MGF of the sum P

Y_k we may solve for the parameters (τ_S, σ_S)at some point t.

We may also relate this to the skewness matching method in Bølviken and Guillen (2017) by noting that matching variance and skewness is equivalent to solving the equations

Z ∞

−∞

w_i(y)f_Y_S(y)dy= Z ∞

−∞

w_i(y)f^P_Y_k(y)dy (8) for the corresponding weight functions w₁(y) = (y− E[Y_k])² = y² and w₂(y) = (y − E[Y_k])³ =y³. Dierent weight functions will translate to dierent penalization schemes.

Using a weight function like(y−E[Y_k])³ will penalize the tail portion of the density more than for example ln(y) which will penalize the head portion more. By choosing exp(−ty) as our weight function we will then be able to tune which parts of the density we want to emphasize. In particular, our focus will be the 99.5th percentile in Solvency II. This makes the method more exible than similar methods which use a xed weight function.

On Figure 2 equation (8) is solved using the weight function exp(−ty)for dierent values oft with the Cholesky approach described in section 3.3.1. Simulation of the sum directly using Monte Carlo is also included. Here we see that by taking t= 5, or weight function exp(−5y), errors in the tail of the distribution is penalized less. While with t = 0.5 we t the tail much better. For values t = 0.1 and t = 0.02 it may seem that we need to go even further along the tail to a reach an area where the scheme benets from such small values. On the other hand, t= 5 ts the head portion of the distribution quite well in comparison to the smaller values of t. The next step in the process is now to derive approximations for the MGF of the sum P

Y_k and the tting variable Y_S. This will be done through the use of Gauss-Hermite expansions of the integrals.

3.3.1 Cholesky decomposition Now, the MGF of the sum P

Y_k may be written on the following form M^P_Y_k(t) =E

e^−t^P^Y^k

=EhY

exp(−t˜σ_k(exp(X_k)−1))i

= Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp(−tσ˜_k(exp(x_k)−1)) 1

(2π)^K/2|C|^1/2 exp

−1

2[x−(−τ²/2)]^TC⁻¹[x−(−τ²/2)]

dx₁...dx_K (9)

(12)

Figure 2: Main window: CCDFs resulting from solving the problem for dierent values of t. Top right window: Corresponding CDFs. Details on page 34 - Program 1

whereτ = (τ_k),X= (X_k)is the Gaussian vector withX_k =−τ_k²/2 +τ_k_k and the matrix C is the covariance matrix of X.

There is no general closed-form expression for (7) and (9) but we may apply Gauss- Hermite approximations to the integrals for real values of t. For the 1-dimensional case of (7) we perform a standard change of variable to get an integrand compatible with Gauss-Hermite integration. We then end up with the following approximation.

M_Y_k(t)≈Ψ_Y_k(t;σ_k, τ_k) =

N

X

n=1

wn

√π exph

−t˜σ_k

exp√

2τ_ka_n−τ_k²/2

−1i

(10) a_n and w_n are the abscissas and weights of the Gauss-Hermite expansion of order N, respectively. The derivation of (10) may be found in the appendix section B on page 18.

The MGF of the sum is more complicated because of the covariance matrix C which we need to get rid of before Gauss-Hermite expansion may be applied. One way of doing this is by applying the following Cholesky transformation as done in Mehta et al. (2007)

xk=√ 2

K

X

i=1

lkizi−τ_k²/2, k= 1, ..., K

where L = (l_ki)^K_k,i=1 is the Cholesky decomposition of C, i.e. C = LL^T. This removes the covariance matrix from the last exponent in (9) and makes the integration variables independent in a sense. After iteratively applying Gauss-Hermite expansion with respect to the variables of integration this approach yield the following approximation to the MGF of the sum P

Y_k. M^P_Y_k(t)≈Ψ^P_Y_k(t;σ,τ, C)

=

N

X

nK=1

· · ·

N

X

n1=1

w_n₁· · ·w_n_K π^K/2

K

Y

k=1

exp (

−t˜σ_k

"

exp √ 2

K

X

j=1

l_kja_n_j−τ_k²/2

!

−1

#)

(11)

(13)

Details in the appendix section C on page 19. Then given parameters τ, σ, σ_S and covariance matrix C we may take some value of t and solve

Ψ^P_Y_k(t;σ,τ, C) = Ψ_Y_S(t;σ_S, τ_S) (12) with respect to τ_S. σ_S is already given by

σ²_S =VarX Y_k

=

K

X

i,j=1

cov(Y_i, Y_j) =

K

X

i,j=1

ρ_ijσ_iσ_j.

3.3.2 Spectral decomposition

We will now introduce an alternative to the expression in (11) based on spectral decomposition rather than Cholesky decomposition. This will open up the possibility of easing the computational burden by reducing the number of sums in an expression like (11) from K to some smaller number L. The expression in (11) sums N^K products. In Solvency II we will encounter K = 12 dimensional problems and these computations will quickly become non-trivial even for a modern computer.

Consider the correlation matrix of ,Σ, and its spectral decomposition Σ = ΓΛΓ^T

where thei-th column vector ofΓis thei-th eigenvector ofΣwith corresponding eigenvalue λ_isorted by descending eigenvalue. Λis the matrix with eigenvalues along the diagonal in descending order. Instead of the Cholesky decomposition used above we instead consider the transformation

δ = Γ^T or = Γδ.

We may then write the model using this transformation Y_k = ˜σ_k exp −τ_k²/2 +τ_k

K

X

j=1

γ_kjδ_j

!

−1

!

with γ_kj being the element (k, j) in Γ. This is useful because since δ ∼ N(0,Λ) the elements δ_j are independent. We may also write the model as

Y_k = ˜σ_k exp −τ_k²/2 +τ_k

K

X

j=1

γ_kjp λ_jη_j

!

−1

!

where η∼N(0, I). Using this expression the MGF of the sum P

Y_k can be written as M^P_Y_k(t) =

Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp (

−tσ˜_k exp

"

−τ_k²/2 +τ_k

K

X

j=1

γ_kjp λ_jη_j

#

−1

!)

× 1

(2π)^K/2 exp

−1 2η^Tη

dη₁...dη_K. (13) Then utilizing a small change of variable and applying Gauss-Hermite expansion iteratively, as done in appendix D on page 20, yields

M^P_Y_k(t)≈Ψˆ^P_Y_k(t;σ,τ,Σ)

=

N

X

nK=1

· · ·

N

X

n1=1

w_n₁· · ·w_n_K π^K/2

K

Y

k=1

exp (

−t˜σ_k exp

"

−τ_k²/2 +τ_k

K

X

j=1

γ_kjp 2λ_ja_n_j

#

−1

!) . (14)

(14)

3.3.3 Dimension reduction

The spectral decomposition gives us the option of dropping the smallest principal components and only including the L < K largest ones. This gives us the following approximation.

_k=

K

X

j=1

γ_kjδ_j ≈

L

X

j=1

γ_kjδ_j

The idea is that for L < k ≤ K, Var(δ_k) = λ_k will be close to zero and δ_k will be approximately N(0,0). Moving forward with this yields the approximate MGF of the sum

M^PY_k(t)≈ Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp (

−t˜σk exp

"

−τ_k²/2 +τk L

X

j=1

γkj

pλjηj

#

−1

!)

× 1

(2π)^L/2exp

−1 2η^Tη

dη₁...dη_L (15) and the Gauss-Hermite expansion

M^P_Y_k(t)≈Ψˆ^P_Y_k(t;σ,τ,Σ)

=

N

X

nL=1

· · ·

N

X

n1=1

w_n₁· · ·w_n_L π^L/2

K

Y

k=1

exp (

−t˜σ_k exp

"

−τ_k²/2 +τ_k

L

X

j=1

γ_kjp 2λ_ja_n_j

#

−1

!) . (16) The eciency gain in using (16) instead of (14) will be examined below.

4 Numerical experiments

4.1 Accuracy

In Figure 3² the absolute error of the dierent methods are compared in a problem with K = 6 risks and K = 12 risks. The x-axis denotes the number of large risks with a large scale (σ_k = 4) and heavy skewness (τ_k = 1.2). The rest of the K risks where smaller and with much lighter skewness (σ_k = 1 and τ_k = 0.25). For each set of parameters the resulting 99.5th percentiles where compared with that of a Monte Carlo simulation of the sum with a size of 10⁷ simulations. For theK = 6 problem at-value of0.1where chosen, while in the K = 12 problem t was set to 0.08 for both the Cholesky and full Spectral (PC) approach. For the PC model with L=K−1 t was set equal to 0.05. These values are believed to provide good results based on analysis in section 4.2. The correlation matrices of (_k) where set equal to the 6- and 12-dimensional correlation matrices from Solvency II documentation, see EIOPA (2014) or page 23. In Solvency II these matrices are meant to dene correlation on the risk-level, i.e. correlation of (Y_k). However, the resulting correlation matrix of (_k), given the correlation of (Y_k), will not always end up being positive denite for any choice of σk and τk. We avoid this by directly dening Σ = cor(_k) equal to the correlation matrix in Solvency II. By starting with a valid correlation structure of (_k) the resulting correlation matrix of (Y_k) will be necessarily positive denite. In this experiment cor(Yk)will change for each set of(τk)and(σk)along

2This section refers to gures placed towards the end of the document on page 14.

(15)

the x-axis. However, the resulting matrices are easily calculated through the results in the appendix on page 23, specically using result (E.8), or via the R-code on page 35.

Firstly note the poor performance of the classical Solvency II Standard formula. The Standard formula especially has problems with the heavily skewed risks. Since the risks stray further and further away from normality asτ_kincreases this is expected. The method introduced in Bølviken and Guillen (2017) provide very good results in the case of only lightly skewed risks. However, when larger, more heavy-tailed risks are introduced the Bølviken-Guillen (BG) method deteriorates abruptly in performance. This method along with the Standard formula seem to do worse and worse as more and more heavy-tailed risks are introduced. The BG method also seems to have trouble when introduced to a large number of risks as in the K = 12 problem.

In the K = 6 problem in Figure 3a the MGF-methods arguably perform better than BG as they handle the large risks better. Figure 3a might also suggest that the values of σ_k and τ_k plays a part in determining the value of t that minimizes the absolute error.

The K = 12 problem provides similar results. The BG method performs very well when all risks are small and lightly skewed. However, when introduced to large risks the method again immediately suers. The MGF methods also show a slightly similar tendency of performing the best when only having small risks with small skewness. Else the absolute accuracy seems relatively constant as we introduce more and more large- scaled risks. This shows considerable improvement on BG and especially the Standard formula. Using the PC-decomposition and discarding the last component also yield good results, albeit slightly worse than using the full model L = K. One can note for zero large risks the optimal t in the L = K −1 case is relatively far from the selected value t = 0.05in 3b which explains the low accuracy in that point. The sensitivity in the design parameter t will be discussed below in section 4.2.

One may dierentiate the large and small risks even more by building a new correlation structure. Consider for example letting the large risks be more strongly correlated with each other and much less correlated with the smaller ones. However, this does not seem to have a great impact on the performance of the dierent methods. This was tested, but not plotted in this document.

The MGF methods also seem to scale down to small-dimensional problems well with K as low as 3 being tested. The resulting plots showed a similar story to that of Figure 3 (Not plotted).

4.2 Design of the MGF method

When applying the MGF method to a problem with given parameters we need to make a choice of N, the Gauss-Hermite degree, and the point t for which to solve (12). We start by considering t which will dictate what part of the density we will approximate the best as discussed on page 6.

Consider theK = 12problem on Figure 4 where dierentt-values were tested and the corresponding 99.5th percentile is plotted for the MGF methods. An even mix of large and small risks were used as in section 4.1. Increasingt makes the resulting density of the new variable Y_S t the lower end of the distribution well, but the 99.5th percentile will be underestimated. On the other hand, a t-value that is too small will provide a good t further out than our target percentile. In which case the lower part of the distribution will be aected by this which will mean that our target 99.5th percentile will be overestimated.

For this problem we seem to hit a sweet spot somewhere between 0.05and0.1. Choice of

(16)

N seem to interact with the resulting plots when comparing theN = 4 case in Figure 4b and the N = 3 case in Figure 4a.

Discarding principal components has a noticeable eect, especially on low t-values.

However, decent accuracy is still obtainable in the K = 12 problem. Another interesting aspect is how much the Cholesky approach degrades for very high t-values compared to the PC approach. This is also observable in the K = 6 problem in Figure 3a.

For theK = 6 problem the similar balanced large-small risk distributions were tested along with letting every risk be small. For the balanced scenario, the results were similar to above with the exception of the reduced model L = 6 which performed very close to the full models. When using only light risks both the Cholesky and PC approach (with L=K) are very exible in terms oft-values that yield near-perfect results. On Figure 5b we observe that any t-value up to 0.5 yield great results, with signicant error building up only with values higher than this. This suggests that the nature of the risks has a signicant impact on how the methods react to dierent values of t.

4.3 Computing time experiments

Along with K, which is already given by the problem one will be facing, N will dictate the computational complexity of the MGF methods. By considering the rst K sums in the expression in (11) we already end up with N^K terms. This is before including the product and sum in the middle and right end of the expression. With this in mind, N and K will both dramatically increase the number of operations necessary to evaluate the resulting expressions. Similar observations are valid for the Spectral approach with the expressions (14) and (16). Minimizing N as much as possible while still achieving a reasonable accuracy should therefore be a high priority.

The time needed for the dierent procedures in the K = 12 and K = 7 problems is plotted on Figure 6 both for the Cholesky approach and the Spectral approach with dierent values of L. This is shown as a function of N. For this problem an even mix of large, heavy-tailed risks and smaller, lightly-tailed risks were chosen. The details may be found in the R-code on page 42 and 44. Here,twas set equal to1. Note that changing the risk parametersσ_k,τ_k,Σorthave negligible eect on the time spent computing. However, it is worth mentioning that not all values of t provide the problem with a solution. Low values of Lin particular does not yield a solution to (12) for low values of t.

First we will consider theK = 12problem on Figure 6a. The computation time seems to suddenly spike for the Cholesky approach and Spectral approach for large values of L. This makes sense considering the number of terms N^K (or N^L) noted above. Note that the Cholesky approach seems to be more ecient than using a Spectral decomposition with L = K. However, by letting L = 11 we may potentially reduce the time required even further. For the K = 12 problem going from N = 3 toN = 4 represents a big jump in computational time for the Cholesky and Spectral approaches. Yet, we notice from section 4.2 that N = 3can achieve good results with less than 30 seconds of computation time.

Now for the smaller K = 7problem in Figure 6b we obtain reasonably small computation times even for very large N. Further experiments indicate that here is little to no improvement in increasing N further than 4. For the K = 7 this problem will then be solved in 1 second. We still note that the Cholesky approach edges out the full Spectral approach in terms of time needed to solve the problem.

(17)

5 Conclusion

The MGF methods consistently provide the best results in terms of accuracy. In the previous section the MGF methods have been shown to provide good accuracy in across a number of dierent scenarios. We have looked at scenarios with only lightly skewed risks, with only heavily skewed risks and with a mixture of lightly and heavily skewed risks. The method of Bølviken and Guillen (2017) handles the rst scenario well, but its accuracy suers when introduced to more heavily skewed risks. The MGF methods presented here does not seem to have this problem but instead produce consistently good results across all of the scenarios.

The Standard formula handled the log-normal risks quite poorly, but on the other hand it is, along with the BG method, quite simple to implement. The MGF methods, on the other hand, are relatively demanding to implement and require a design parameter t to be set. The sensitivity of the choice of t is dependent on the risks and should be studied further if one were to consider implementing the MGF methods into a system like Solvency II which handles dierent types of risk. EIOPA would either have to set the value of t pre-determinately or provide some form of automatic procedure which would provide a t-value. If an insurance company could freely choose its t-values then it could knowingly underestimate the SCRs by intentionally setting them too high as shown on Figure 2.

For a problem with a high number of risks computation time for the MGF methods may be a concern for some. The time required was shown to be non-trivial in the 12- dimensional problem. A dimension reduction method was developed with such problems in mind and tested along with the full models. Combining this with reducing the Gauss- Hermite parameter N to 3 was shown to signicantly reduce computation time while still oering good accuracy mostly on par with the full-model approaches.

(18)

References

Bølviken, E., (2014) Computation and modelling in insurance and nance. Cambridge University Press, page 150.

Bølviken, E., (2017) "Solvency II quickly". In preparation.

Bølviken, E., & Guillen, M. (2017). Risk aggregation in Solvency II through recursive log-normals. Insurance: Mathematics and Economics, Vol. 73, pages 20-26.

Commission Delegated Regulation (2015). In the Ocial Journal of the European Union. Available as http://eur-lex.europa.eu/legal-content/EN/TXT/PDF/?uri=

OJ:L:2015:012:FULL&from=EN

EIOPA (2014). Technical Specication For The Preparatory Phase (Part I).

Available as https://eiopa.europa.eu/Publications/Standards/A_-_Technical_

Specification_for_the_Preparatory_Phase__Part_I_disclaimer.pdf

Mehta, N. B., Wu, J., Molisch, A. F., & Zhang, J. (2007). Approximating a sum of random variables with a lognormal. IEEE Transactions on Wireless Communications, Vol.6, No. 7, pages 2690-2699.

Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. (2007). Numerical recipes: The Art of Scientic Computing. 3rd ed. Cambridge University Press, Cam- bridge.

(19)

Appendices

A Figures - Numerical experiments

(a) The K= 6 problem. N = 4.

(b) The K= 12 problem. N = 3

Figure 3: Distance between the 99.5th percentiles of the Monte Carlo simulations and for the methods in the K = 6 and K = 12 problems with increasing number of large, heavy-tailed risks.

(20)

(a) N = 3.

(b) N = 4.

Figure 4: Resulting 99.5th percentiles for dierent values of t in the K = 12 problem.

Note non-linear x-axis.

(21)

(a) 3 heavy risks with (σ_k, τ_k) = (4,1.2) and 4 light risks with (σ_k, τ_k) = (1,0.25).

(b) Only light risks

Figure 5: Resulting 99.5th percentiles for dierent values of t in the K = 7 problem.

Note non-linear x-axis.

(22)

(a) The K = 12 problem. Using 6 risks with (σ_k, τ_k) = (4,1.2) and 6 with (σ_k, τ_k) = (1,0.25).

(b) The K = 7 problem. Using 3 risks with (σ_k, τ_k) = (1.4,4) and 4 with (σk, τk) = (0.25,1).

Figure 6: Time used to solve Ψ^P_Y_k = Ψ_Y_S for dierent values N. Comparing the Cholesky approach to the spectral decomposition approach for dierent values of L. Using correlation matrix from Solvency II documentation.

(23)

B PDF and MGF

B.1 PDF of our log-normal model

Our model variable Y is given as a function of a regular log-normal variable V through the transformation g

Y =σe^−τ²^/2+τ−1

√

e^τ² −1 = ˜σ(V −1) =g(V) (B.1) whereσ˜= σ

√

e^τ² −1 andV =e^−τ²^/2+τ ∼lnN(−τ²/2, τ). Then by the change of variable formula the density of Y is readily obtained through

f_Y(y) =

d

dyg⁻¹(y))

f_V(g⁻¹(y))

= 1

˜ σf_V y

˜ σ + 1

(B.2) where

f_V(v) = 1 vτ√

2π exp

−(lnv−(−τ²/2))² 2τ²

. The support of Y is (−˜σ,∞).

B.2 The moment-generating function of Y and Gauss-Hermite

The moment-generating function of a random variable is dened as M_Y(t) =E

e^−tY

, t∈R

wherever this expectation exists. Using equation (B.2) this becomes M_Y(t) =E

e^−tY

= Z ∞

−˜σ

exp(−ty)f_Y(y)dy

= Z ∞

−˜σ

exp(−ty)1

˜ σ

1 (y

˜

σ + 1)τ√ 2π

exp



−

(ln(y

˜

σ + 1) +τ²/2)² 2τ²



dy. (B.3) In order to end up with the expression in (10) we need to use substitution to rst arrive at an integrand compatible with Gauss-Hermite integration. Set

z = ln(y/˜σ+ 1) +τ²/2

√2τ giving y= ˜σh exp(√

2τ z−τ²/2)−1i . We then get

dz

dy = 1

√2τ(˜σ+y) ordy =√

2τ(˜σ+y)dz.

Finally, consider the integration limits. If y → −˜σ then ln(y

˜

σ + 1) → −∞and z → −∞. Now if y → ∞then ln(y

˜

σ + 1)→ ∞ and z → ∞. Therefore the integral (B.3) becomes M_y(t) =

Z ∞

−∞

√1

πexph

−tσ˜ exp(√

2τ z−τ²/2)−1i

exp(−z²)dz

≈

N

X

n=1

wn

√πexp h

−tσ˜

exp(√

2τ x_n−τ²/2)−1

i (B.4)

(24)

which is the expression in (10).

C The Cholesky decomposition

C.1 Change of variable

We have the full expression of the MGF of P

Yk as on page 6:

M^P_Y_k(t) = Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp(−t˜σ_k(exp(x_k)−1)) 1

(2π)^K/2|C|^1/2 exp

−1

2[x−(−τ²/2)]^TC⁻¹[x−(−τ²/2)]

dx₁...dx_K (C.1) We will now do two things with one change of variable. Firstly, we want to obtain inde- pendence in the integrand, i.e. we want expression like (C.1) but where x independent.

Secondly, we want to get an integrand compatible with Gauss-Hermite integration. In order to achieve the rst goal, we consider the Cholesky decomposition in order to decor- relate x. X = (X_k) is the Gaussian vector X_k = −τ_k²/2 +τ_k_k and the matrix C is the covariance matrix ofX. LetC =LL^T be its Cholesky decomposition. We then know that if we consider an arbitrary vector z with mean zero and covariance equal to the identity matrix that the transformation Lz will be a mean-zero vector with covariance matrix equal to C. If we furthermore multiply by √

2 and add−τ²/2 to the transformation our integrand will be Gauss-Hermite compatible. Our change of variable then becomes

x=√

2Lz−τ²/2 or

x_k=√ 2

K

X

i=1

l_kiz_i−τ_k²/2, k = 1, ..., K.

We then get dx_k dz_i =√

2l_ki and the Jacobian becomes J(z) =√

2Land |J(z)|= 2^K/2|L|= 2^K/2|C|^1/2. The MGF then becomes

M₍^P_Y_k₎(t) =E[e^−t^P^Y^k] = Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp (

−t˜σ_k exp

"

√2

K

X

i=1

l_kiz_i−τ_k²/2

#

−1

!)

× 1

(2π)^K/2|C|^1/2 exp

−1 2(√

2Lz)^TC⁻¹(√ 2Lz)

2^K/2|C|^1/2dz₁...dz_K

= Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp (

−t˜σ_k exp

"

√ 2

K

X

i=1

l_kiz_i−τ_k²/2

#

−1

!)

exp(−z^Tz)dz₁...dz_K. (C.2)

C.2 Applying Gauss-Hermite

We factor out z1 and apply Gauss-Hermite expansion with respect to z1:

(25)

M₍^P_Y_k₎(t) = Z ∞

−∞

· · · Z ∞

−∞

1

π^(K−1)/2 exp −

K

X

i=2

z_i²

! _N X

n1=1

w_n₁

π^1/2f₍₁₎(z⁽⁻¹⁾, an1;t)dz2...dzK+R⁽¹⁾_N where we have some error term in the end and

f₍₁₎(z⁽⁻¹⁾, a_n₁;t) =

K

Y

k=1

exp −t˜σ_k

"

exp(√ 2

K

X

i=2

l_kiz_i+√

2l_i1a_n₁ −τ_k²/2)−1

#!

. Then again with z₂

M₍^P_Y_k₎(t)

= Z ∞

−∞

· · · Z ∞

−∞

1

π^(K−2)/2 exp −

K

X

i=3

z²_i

! _N X

n2=1 N

X

n1=1

w_n₁w_n₂

π^2/2 f₍₂₎(z⁽⁻²⁾, a_n₁, a_n₂;t)dz₃...dz_K +R⁽²⁾_N where

f₍₂₎(z⁽⁻²⁾, a_n₁, a_n₂;t) =

K

Y

k=1

exp −t˜σ_k

"

exp(√ 2

K

X

i=3

l_kiz_i+√ 2

2

X

m=1

l_ima_n_m−τ_k²/2)−1

#!

. By continuing this we end up the expression

M₍^P_Y_k₎(t)

=

N

X

nK=1

· · ·

N

X

n1=1

w_n₁· · ·w_n_K π^K/2

K

Y

k=1

exp (

−t˜σ_k

"

exp √ 2

K

X

j=1

l_kja_n_j−τ_k²/2

!

−1

#)

+R_N^(K). (C.3) This then justies the expression in (11) where the error term is discarded.

D The Spectral decomposition

D.1 Change of variables

The strategy and technique is the same as with the Cholesky decomposition. We start with the expression in (13) and apply a small change of variables before we iteratively apply Gauss-Hermite expansion with respect to the variables of integration. We start with

M^P_Y_k(t) = Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp (

−t˜σ_k exp

"

−τ_k²/2 +τ_k

K

X

j=1

γ_kjp λ_jη_j

#

−1

!)

× 1

(2π)^K/2 exp

−1 2η^Tη

dη1...dηK. Settingη=√

2zwill remove the fraction in the last exponent. This gives us the Jacobian J(z) = √

2I and |J(z)|= 2^K/2 and the expression above becomes

(26)

Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp (

−tσ˜_k exp

"

−τ_k²/2 +τ_k

K

X

j=1

Γ_kjp 2λ_jz_j

#

−1

!)

× 1

π^K/2 exp −

K

X

k=1

z_k²

!

dz₁...dz_K.

D.2 Applying Gauss-Hermite

First with respect to z₁: M^P_Y_k(t) =

Z ∞

−∞

· · · Z ∞

−∞

1 π^(K−1)/2

N

X

n1=1

w_n₁

π^1/2f₍₁₎(z⁽⁻¹⁾, a_n₁) exp −

K

X

k=2

z_k²

!

dz₂...dz_K +R⁽¹⁾_N where we have some error term in the end and

f₍₁₎(z⁽⁻¹⁾, a_n₁) =

K

Y

k=1

exp (

−t˜σ_k exp

"

−τ_k²/2 +τ_k

K

X

j=2

Γ_kjp

2λ_jz_j +τ_k

1

X

j=1

Γ_kjp 2λ_ja_n_j

#

−1

!) . Then again with respect z₂:

M^P_Y_k(t) = Z ∞

−∞

· · · Z ∞

−∞

1 π^(K−2)/2

N

X

n2=1 N

X

n1=1

w_n₁w_n₂

π^2/2 f₍₂₎(z⁽⁻²⁾, a_n₁, a_n₂) exp −

K

X

k=3

z_k²

!

dz₂...dz_K+R⁽²⁾_N where

f₍₂₎(z⁽⁻²⁾, a_n₁, a_n₂) =

K

Y

k=1

exp (

−t˜σ_k exp

"

−τ_k²/2 +τ_k

K

X

j=3

Γ_kjp

2λ_jz_j +τ_k

2

X

j=1

Γ_kjp 2λ_ja_n_j

#

−1

!) . Continuing this yields

M^P^K_Y

k(t) =

N

X

nK=1

· · ·

N

X

n1=1

wn1· · ·wn_K

π^K/2

K

Y

k=1

exp (

−tσ˜_k exp

"

−τ_k²/2 +τ_k

K

X

j=1

Γ_kjp 2λ_ja_n_j

#

−1

!)

+R^(K)_N (D.4) which then justies the expression in (14) where the error term is discarded.

D.3 Dimension reduction

The same strategy applies in this case, but still written out here for completeness. We start with expression (15).

M^P_Y_k(t)≈ Z ∞

−∞

· · · Z ∞

−∞

K

Y

k=1

exp (

−t˜σ_k exp

"

−τ_k²/2 +τ_k

L

X

j=1

γ_kjp λ_jη_j

#

−1

!)

× 1

(2π)^L/2 exp

−1 2η^Tη

dη₁...dη_L (D.5)

Log-normal approximations in the Solvency II Standard Formula