The Skewness of Bornhuetter-Ferguson

1. Introduction

After the famous chain-ladder method, the Bornhuetter-Ferguson method (hereinafter “BF method”—see Bornhuetter and Ferguson 1972) is one of the methods practicing actuaries most use for the projection of non-life paid or incurred triangles. The BF method relies on two sets of parameters:

• the claims development patterns derived from the incurred or paid triangles; and

• the a priori estimates for the ultimate claims amount—such estimates can come from pricing models or from any other sources.

Using those sets of parameters, the BF method estimates the best estimate of the claims liability.

Following the development of the prediction error estimate for the chain-ladder method (Mack 1993, 1999), Mack (2008) proposed the estimate of the prediction error for BF. In addition, different skewness estimates for the chain-ladder method were developed (Salzmann, Wuthrich, and Merz 2012; Dal Moro 2013). All such estimates are done in a distribution-free framework.

To complete our knowledge of the BF method in a distribution-free framework, the missing piece is the skewness of the BF method. With knowledge of the first three moments, we can then estimate the different quantiles of the reserving distributions (see Dal Moro and Krvavych 2017). Such quantiles will certainly be useful in the context of International Financial Reporting Standard 17 (IFRS 17), which directs (re)insurance companies to publish additional quantile disclosures.

In this paper, we propose an approach to estimating this skewness relying mainly on the stochastic BF framework proposed by Mack (2008). The paper is divided into the following sections:

• Section 2 describes the stochastic BF model of Mack (2008) and extends it to cover the assumptions necessary to estimate the third moment.

• Section 3 provides an estimate of the skewness of the BF method per accident year.

• Section 4 provides the estimate of the skewness over all accident years.

• A final section provides a few numerical examples.

Remark: An Excel sheet developed to estimate the presented formulae is available at https://drive.google.com/open?id=1iRPEnD8eVOECPyR4oZl_Ezd89agVHYAa.

2. The stochastic model underlying the BF method

Let \(C_{i,k}\) denote the cumulative claims amount (either paid or incurred) of accident year i after k years of development, \(1 \leq i,\ k \leq n,\) and \(\nu_{i}\) be the premium volume of accident year i where n denotes the most recent accident year. Then \(C_{i,n + 1 - i}\) denotes the currently known claims amount of accident year i. Let further \(S_{i,k} = C_{i,k} - C_{i,k - 1}\) denote the incremental claims amount (with \(C_{i,0} = 0\)) and U_i the (unknown) ultimate claims amount of accident year i. Then \(R_{i} = U_{i} - C_{i,n + 1 - i}\) is the (unknown true) claims reserve for accident year i. Finally, let \(S_{i,n + 1} = U_{i} - C_{i,n}\) be the incremental claims amount after development year n (tail development).

Bornhuetter and Ferguson (1972) introduced their method to estimate R_i as follows:

\[{\widehat{R}}_{i}^{\text{BF}} = {\widehat{U}}_{i}\left( 1 - {\widehat{z}}_{n + 1 - i} \right),\]

where \({\widehat{U}}_{i} = \nu_{\text{i}}{\widehat{q}}_{i}\) with a prior estimate \({\widehat{q}}_{i}\) for the ultimate claims ratio \(q_{i} = \frac{U_{i}}{\nu_{i}}\) of accident year \(i\), \({\widehat{z}}_{k} \in \lbrack 0;1\rbrack\) is the estimated percentage of the ultimate claims amount that is expected to be known after development year k.

The BF stochastic model developed in Mack (2008) relies on the following assumptions related to the increments \(S_{i,k}\):

• BF1: All increments \(S_{i,k}\) are independent.

• BF2: There are unknown parameters x_i, y_k such that

  • \(E\left( S_{i,k} \right) = x_{i}y_{k}\); and

  • \(y_{1} + \ldots + y_{n + 1} = 1\).

• BF3: There are unknown proportionality constants \(s_{k}^{2}\) with \(\text{Var}\left( S_{i,k} \right) = x_{i}s_{k}^{2}\).

On the basis of these three assumptions, one can estimate the prediction error of BF (see Mack 2008). The prediction error, usually denoted as the MSEP (mean squared error of prediction), consists of two components: the process error and the estimation error. Whereas one can basically always calculate the estimation error via the laws of error propagation, for the process error Mack (2008) developed a stochastic model of the claims process.

To estimate the skewness of the BF method, we need a fourth assumption, which is derived by analogy with the skewness estimation of the chain-ladder model (see Dal Moro 2013):

• BF4: There are unknown proportionality constants \(t_{k}^{3}\) with \(K\left( S_{i,k} \right) = x_{i}^{\frac{3}{2}}t_{k}^{3}\) where \(K\left( S_{i,k} \right) = E\left\lbrack \left( S_{i,k} - E\left( S_{i,k} \right) \right)^{3} \right\rbrack\).

Following Mack (2008), we have the following with x₁,… x_n known:

\[\ {\widehat{y}}_{k} = \frac{\sum_{i = 1}^{n + 1 - k}S_{i,k}}{\sum_{i = 1}^{n + 1 - k}x_{i}}\tag{1}\]

is a best linear unbiased estimate of y_k, and

\[{\widehat{s}}_{k}^{2} = \frac{1}{n - k}\sum_{i = 1}^{n + 1 - k}\frac{\left( S_{i,k}\ - \ x_{i}{\widehat{y}}_{k} \right)^{2}}{x_{i}} \tag{2}\]

is an unbiased estimate of s²_k.

And, following BF4:

\[{\widehat{t}}_{k}^{3} = \frac{1}{n - k}\sum_{i = 1}^{n + 1 - k}\frac{\left( S_{i,k}\ - \ x_{i}{\widehat{y}}_{k} \right)^{3}}{x_{i}^{\frac{3}{2}}} \tag{3}\]

is an unbiased estimate of t³_k.

Having estimated the parameters \({\widehat{y}}_{k},\ {\widehat{s}}_{k}^{2},\ and\ {\widehat{t}}_{k}^{3}\) (see Mack (2008) for details on the estimation process), we can calculate the BF claims reserve by

\({\widehat{R}}_{i}^{\text{BF}} = {\widehat{U}}_{i}\left( {\widehat{y}}_{n + 2 - i} + ... + {\widehat{y}}_{n + 1} \right) = {\widehat{U}}_{i}\left( 1 - {\widehat{z}}_{n + 1 - i} \right)\) with \({\widehat{z}}_{k} = {\widehat{y}}_{1} + ... + {\widehat{y}}_{k}.\)

The properties of these estimators are given below:

• \({\widehat{y}}_{1},...,{\widehat{y}}_{n + 1}\) are pairwise (slightly) negatively correlated as they have to add up to unity.

• \({\widehat{y}}_{1},...,{\widehat{y}}_{n + 1},\) and \({\widehat{z}}_{1},...,{\widehat{z}}_{n + 1}\) are independent from \({\widehat{U}}_{1},...,\ {\widehat{U}}_{n}\).

• \({\widehat{R}}_{i}^{\text{BF}}\) and R_i are independent (due to BF1).

• \(E\left( {\widehat{U}}_{i} \right) = E\left( U_{i} \right) = x_{i}\text{ for }1 \leq i \leq n\).

• \(E\left( {\widehat{y}}_{k} \right) = y_{k}\text{ for }1 \leq k \leq n + 1\) and therefore
  \(E\left( {\widehat{z}}_{k} \right) = z_{k}\text{ for }1 \leq k \leq n + 1\).

• \(E\left( {\widehat{s}}_{k}^{2} \right) = s_{k}^{2}\text{ for }1 \leq k \leq n + 1\).

• \(E\left( {\widehat{t}}_{k}^{3} \right) = t_{k}^{3}\text{ for }1 \leq k \leq n\).

For the last four bullet points, we simply assume that the actuary’s selections are unbiased.

Having described the stochastic BF model and its assumptions, we devote the next section to estimate the skewness of the BF method for one accident year.

3. Skewness of the BF method per accident year

The mean skewness of any reserve estimate \({\widehat{R}}_{i}\) is defined as

\[\text{SKEW}\left( {\widehat{R}}_{i} \right) = E\left\lbrack \left( {\widehat{R}}_{i} - R_{i} \right)^{3}│S_{i,1},...,S_{i,n + 1 - i} \right\rbrack.\]

Due to BF1, \({\widehat{R}}_{i}^{\text{BF}}\) and R_i are taken to be commonly independent from \(S_{i,1},...,S_{i,n + 1 - i}\). Hence,

\[\text{SKEW}\left( {\widehat{R}}_{i}^{\text{BF}} \right) = E\left\lbrack \left( {\widehat{R}}_{i}^{\text{BF}} - R_{i} \right)^{3} \right\rbrack.\]

In Appendix A, it is proved that

\[\text{SKEW}\left( {\widehat{R}}_{i}^{\text{BF}} \right) = K\left( {\widehat{R}}_{i}^{\text{BF}} \right) - K\left( R_{i} \right)\tag{4}.\]

Following BF4, we have

\[K\left( R_{i} \right) = K\left( S_{i,n + 2 - i} \right) + ... + K\left( S_{i,n + 1} \right) = x_{i}^{\frac{3}{2}}\left( t_{n + 2 - i}^{3} + ... + t_{n + 1}^{3} \right),\]

which is estimated by

\[\widehat{K}\left( R_{i} \right) = {\widehat{U}}_{i}^{\frac{3}{2}}\left( {\widehat{t}}_{n + 2 - i}^{3} + ... + {\widehat{t}}_{n + 1}^{3} \right).\]

As for \(K\left( {\widehat{R}}_{i}^{\text{BF}} \right) = K\left\lbrack {\widehat{U}}_{i}\left( 1 - {\widehat{z}}_{n + 1 - i} \right) \right\rbrack\), we use the formula below, which holds when X and Y are independent:

\(\begin{align} K\left(XY\right) &= K(X)K(Y) + K(X)E(Y) \lbrack{3\text{Var}(Y) + E(Y)^{2}}\rbrack \\&\ \ \ \ + K(Y)E(X)\lbrack {3\text{Var}(X) + E(X)^{2}} \rbrack \\&\ \ \ \ + 6E(X)E(Y)\text{Var}(Y)\text{Var}(X). \end{align}\)

Using the above formula, we have

\(\begin{align} K\left( {\widehat{R}}_{i}^{\text{BF}} \right) =& - K\left\lbrack {\widehat{U}}_{i} \right\rbrack K\left( {\widehat{z}}_{n + 1 - i} \right) \\&+ K\left\lbrack {\widehat{U}}_{i} \right\rbrack E\left( 1 - {\widehat{z}}_{n + 1 - i} \right)\left\lbrack 3Var\left( {\widehat{z}}_{n + 1 - i} \right) + E\left( 1 - {\widehat{z}}_{n + 1 - i} \right)^{2} \right\rbrack\ \\& - K\left\lbrack {\widehat{z}}_{n + 1 - i} \right\rbrack E\left( {\widehat{U}}_{i} \right)\left\lbrack 3Var\left( {\widehat{U}}_{i} \right) + E\left( {\widehat{U}}_{i} \right)^{2} \right\rbrack \\& + 6\ E\left( {\widehat{U}}_{i} \right)E\left( 1 - {\widehat{z}}_{n + 1 - i} \right)\text{Var}\left( {\widehat{U}}_{i} \right)\text{Var}\left( {\widehat{z}}_{n + 1 - i} \right). \end{align}\)

Estimators of \(E\left( {\widehat{U}}_{i} \right),\ Var\left( {\widehat{U}}_{i} \right), E\left( 1 - {\widehat{z}}_{n + 1 - i} \right),\) and \(Var\left( {\widehat{z}}_{n + 1 - i} \right)\) are provided in Mack (2008) (see Appendix B for details on these estimators). There remains to estimate \(K\left\lbrack {\widehat{U}}_{i} \right\rbrack\text{ and K}\left( {\widehat{z}}_{n + 1 - i} \right)\).

As for \(K\left\lbrack {\widehat{U}}_{i} \right\rbrack\), we have to make an assumption on the underlying distribution of \({\widehat{U}}_{i}\). As a generally accepted distribution used in reserving contexts, we will therefore assume a lognormal distribution with parameters μ and σ derived from \(E\left( {\widehat{U}}_{i} \right)\text{ and Var}\left( {\widehat{U}}_{i} \right)\). The parameters and resulting skewness are given by

\(\sigma = \sqrt{\ln\left( 1 + \frac{\text{Var}\left( {\widehat{U}}_{i} \right)}{{E\left( {\widehat{U}}_{i} \right)}^{2}} \right)}.\)

\(K\left\lbrack {\widehat{U}}_{i} \right\rbrack = \left( 2 +\text{exp}\left( \sigma^{2} \right) \right)\frac{{\text{Var}\left( {\widehat{U}}_{i} \right)}^{2}}{E\left( {\widehat{U}}_{i} \right)}\).

As for \(K\left( {\widehat{z}}_{n + 1 - i} \right)\), we note that due to the slightly negative correlations between \({\widehat{y}}_{1},...,{\widehat{y}}_{n + 1}\), we have to make the following approximations:

\(\ K\left( {\widehat{z}}_{k} \right) \approx min\left\lbrack K\left( {\widehat{y}}_{1} \right) + ... + K\left( {\widehat{y}}_{k} \right);K\left( {\widehat{y}}_{k + 1} \right) + ... + K\left( {\widehat{y}}_{n + 1} \right) \right\rbrack.\)

As we have

\[K\left( {\widehat{y}}_{k} \right)\ = K\left( \frac{\sum_{i = 1}^{n + 1 - k}S_{i,k}}{\sum_{i = 1}^{n + 1 - k}{\widehat{U}}_{i}} \right) = \frac{\sum_{i = 1}^{n + 1 - k}{K\left( S_{i,k} \right)}}{\left( \sum_{i = 1}^{n + 1 - k}{\widehat{U}}_{i} \right)^{3}}\]

and

\[\widehat{K}\left( S_{i,k} \right) = {\widehat{U}}_{i}^{\frac{3}{2}}{\widehat{t}}_{k}^{3},\]

we finally get

\[K\left( {\widehat{y}}_{k} \right)\ = \frac{{\widehat{t}}_{k}^{3}\sum_{i = 1}^{n + 1 - k}{\widehat{U}}_{i}^{\frac{3}{2}}}{\left( \sum_{i = 1}^{n + 1 - k}{\widehat{U}}_{i} \right)^{3}}\]

This completes the estimation of the different elements of \(\text{SKEW}\left( {\widehat{R}}_{i}^{\text{BF}} \right)\).

4. Skewness of the BF method over all accident years

Having estimated the skewness of the BF method by accident year, we now aggregate these elements over all accident years. To do so, we use the Fleishman polynomials (see Fleishman 1978). First, we assume that the centralized and normalized copy of the risk value X_i of the ith class, \(\widehat{X_{i}} = \frac{X_{i} - E\left( X_{i} \right)}{E\left( X_{i} \right)\ \text{CoV}_{X_{i}}}\) (where CoV denotes the coefficient of variation), is estimated by the Fleishman polynomial structure of a standard normal random variable. In particular, we consider the following two different cases:

• \({\widehat{X}}_{i} = P_{2}\left( Z_{i} \right) = a_{i}Z_{i} + b_{i}\left( Z_{i}^{2} - 1 \right)\), where \(Z_{i}\) denotes the standard normal distribution. Such a case is suitable for estimating the skewness of a risk portfolio profile when the confidence level is approximated using skewness only.

• \({\widehat{X}}_{i} = P_{3}\left( Z_{i} \right) = a_{i}Z_{i} + b_{i}\left( Z_{i}^{2} - 1 \right) + c_{i}Z_{i}^{3}\). This is suitable for estimating skewness and kurtosis of the risk portfolio risk profile when the confidence level is approximated using both skewness and kurtosis.

The coefficients of the polynomials P₂ and P₃ are calibrated using the method of moments by matching the second and third moments of P₂(Z_i) and the second, third, and fourth moments of P₃(Z_i) to 1 (standard deviation of \(\widehat{X_{i}}\)), γ_i (skewness of X_i), and ι_i + 3 (noncentralized or absolute kurtosis of X_i), respectively.

The coefficients of P₂ can be analytically expressed by solving the following system of equations:

\(\begin{align}\left\{ \begin{matrix} 1 = a_{i}^{2} + 2b_{i}^{2}. \\ \gamma_{i} = 6a_{i}^{2}b_{i} + 8b_{i}^{3}. \\ \end{matrix} \right.\end{align} \tag{5}\)

The system (5) is reduced to

\(\left\{ \begin{matrix} a_{i} = \sqrt{1 - 2b_{i}^{2}.} \\ \gamma_{i} = 6b_{i} - 4b_{i}^{3}. \\ \end{matrix} \right.\tag{6}\)

Such a system is easily solved. The roots of the cubic equation (6) can be found using Cardano’s formula (see, e.g., Abramowitz and Stegun 1972). If we denote

\[\varphi = arccos\left( - \frac{\gamma_{i}^{3}}{\sqrt{8}} \right).\]

Then, the only real root of equation (6) is

\[b_{i} = \sqrt{2}\text{cos}\left( \frac{\varphi}{3} + 4\frac{\pi}{3} \right).\]

Having estimated the above parameters of the Fleishman polynomial, we define the total reserve value across the portfolio of m risks as

\[X_{\Sigma} = \sum_{i = 1}^{m}X_{i}\]

where each ith risk value is approximated by the Fleishman polynomial of a standard normal random variable:

\[X_{i} \approx \text{CE}_{i}\ \left( 1 + \text{CoV}_{i}\ P_{3}\left( Z_{i} \right) \right),\]

where CE denotes the central best estimate of X.

It is clear that \(\text{CE}_{\Sigma} = \sum_{i = 1}^{m}\text{CE}_{i}\). It should also be noted that setting \(c_i = 0\) reduces the problem to simply approximating:

\[X_{i} \approx \text{CE}_{i}\ \left( 1 + \text{CoV}_{i}\ P_{2}\left( Z_{i} \right) \right).\]

As discussed earlier, all the stand-alone risks interact between each other according to a Gaussian dependence structure, of which the linear correlations ρ_ij (coefficients of a Gaussian copula) are given as in Mack (2008):

\[\rho_{\text{ij}} = \frac{{\widehat{z}}_{n + 1 - j}\left( 1 - {\widehat{z}}_{n + 1 - i} \right)}{{\widehat{z}}_{n + 1 - i}\left( 1 - {\widehat{z}}_{n + 1 - j} \right)}.\]

Skewness

We compute the third central moment of X_Σ:

\[E\left\lbrack \left( X_{\Sigma} - \text{CE}_{\Sigma} \right)^{3} \right\rbrack = E\left\lbrack \left( \sum_{i = 1}^{m}{\sigma_{i}\ P_{3}\left( Z_{i} \right)} \right)^{3} \right\rbrack\]

\[\begin{align}E&\left\lbrack \left( X_{\Sigma} - \text{CE}_{\Sigma} \right)^{3} \right\rbrack \\&= \sum_{i = 1}^{m}{\sigma_{i}^{3}\ \gamma_{i}} + 3\sum_{\text{ij}}^{}{\sigma_{i}^{2}\sigma_{j}E\left\lbrack {P_{3}\left( Z_{i} \right)}^{2}P_{3}\left( Z_{j} \right) \right\rbrack} \\&\ \ \ \ + 6\sum_{\text{ijk}}^{}{\sigma_{i}\sigma_{j}\sigma_{\text{k}}E\left\lbrack P_{3}\left( Z_{i} \right)P_{3}\left( Z_{j} \right)P_{3}\left( Z_{k} \right) \right\rbrack}\end{align}\]

where \(E\left\lbrack {P_{3}\left( Z_{i} \right)}^{3} \right\rbrack = \gamma_{i}\), as the Fleishman polynomial coefficients are calibrated so that the polynomial has skewness γ_i for the ith stand-alone risk profile. In formula (7), the summation term with multiple 3 has \(\begin{pmatrix} m \\ 2 \\ \end{pmatrix}\) different subterms, and the summation term with multiple 6 is relevant if \(m \geq 3\) and has \(\begin{pmatrix} m \\ 3 \\ \end{pmatrix}\) different subterms.

The following components of formula (7) above are provided here only for the partial case where P₂ is used—i.e., c_i = 0:

\[E\left\lbrack {P_{3}\left( Z_{i} \right)}^{2}P_{3}\left( Z_{j} \right) \right\rbrack = 2\rho_{\text{ij}}\left( 2a_{i}a_{j}b_{i} + \left( a_{i}^{2} + 4b_{i}^{2} \right)b_{j}\rho_{\text{ij}} \right).\tag{8}\]

\[\begin{align}E&\left\lbrack P_{3}\left( Z_{i} \right)P_{3}\left( Z_{j} \right)P_{3}\left( Z_{k} \right) \right\rbrack \\&= 2\left( a_{j}a_{k}b_{i}\rho_{\text{ij}}\rho_{\text{ik}} + a_{j}a_{i}b_{k}\rho_{\text{jk}}\rho_{\text{ik}} + a_{i}a_{k}b_{j}\rho_{\text{ij}}\rho_{\text{jk}} \right) \\&\ \ \ \ + 8b_{i}b_{k}b_{j}\rho_{\text{ij}}\rho_{\text{ik}}\rho_{\text{jk}}.\end{align}\tag{9}\]

The skewness γ_Σ is then calculated as follows:

\[\gamma_{\Sigma} = \frac{E\left\lbrack \left( X_{\Sigma} - \text{CE}_{\Sigma} \right)^{3} \right\rbrack}{\left( \text{CE}_{\Sigma}\ \text{CoV}_{\Sigma} \right)^{3}}.\]

In the numerical examples that follow, the parameters a_i, b_i are provided, the correlation matrices are provided in Appendix C, and the values for σ_i correspond to the msep(R_i), which are estimated using Mack (2008) formulae. Details of calculations can be found on the Excel sheet available at:

https://drive.google.com/open?id=1iRPEnD8eVOECPyR4oZl_Ezd89agVHYAa.

5. Numerical examples

The formulae above were tested on three triangles (see Appendix C). In all the cases, the coefficient of variation
of \({\widehat{U}}_{i}\) is kept constant at 10%. The correlation matrices used for applying the Fleishman
polynomials are also provided in Appendix C.

For the three triangles, it should be noted that one must add a tail factor as the cumulative patterns do not reach 100% (triangle 1 tail is 33.81%, triangle 2 tail is 32.45%, triangle 3 tail is 36.29%). As for these examples, the tail factors are the simple difference between the cumulated patterns coming from the application of the Mack (2008) method and 100%. Other methods could have been chosen to review the patterns such as increasing proportionally all incremental patterns to reach 100%.

As expected, the overall resulting skewness is positive in all cases. The skewness amounts are also reasonable as the corresponding distributions would not belong to extreme distributions (e.g., Pareto) but to smoother distributions (e.g., lognormal or gamma). Therefore, the proposed model seems to provide robust results on these three triangles.

6. Conclusion

This paper is a first attempt to estimate the skewness of the BF method. It relies on a few assumptions that have to be derived from external knowledge of the modeled reserving risk:

• we assume the value of \({\widehat{U}}_{i}\);

• we assume the coefficient of variation of \({\widehat{U}}_{i}\);

• a lognormal distribution is assumed to estimate \(K\left\lbrack {\widehat{U}}_{i} \right\rbrack\); and

• the slightly negative correlations between \({\widehat{y}}_{1},...,{\widehat{y}}_{n + 1}\) result also in some assumptions being made.

As for the first two points, the application of the proposed formulae could be refined if these parameters can be estimated externally, e.g., with the knowledge of the pricing distributions. However, such parameters are not always easily available. For the value of \({\widehat{U}}_{i}\), however, the application of the Cape Cod method can provide a good starting point.

Overall, refining the above assumptions may provide slightly better results, but we should bear in mind that in the context of IFRS 17 and of the quarterly disclosure requirements, simple and easily implementable models should be favored. Therefore, the proposed formulae will certainly fit with IFRS 17 requirements.

Further to refinements, it should be noted that the proposed skewness estimator is sensitive to the choice of the \({\widehat{U}}_{i}\) and to the tail factor retained for the pattern estimation. The sensitivity to the choice of \({\widehat{U}}_{i}\) is the result of the formulae given in section 3. As for the tail factor, it comes from the application of the parametrization of the BF model by Mack (2008). In general, the cumulative payment or incurred pattern does not reach 100% and, as a result, it is necessary to put a tail factor. This tail factor (which is automatically calculated in this paper) may be different from the one that an actuary would have applied based on his or her judgment. Such differences in the tail factor and in the overall pattern could affect the overall skewness. When applying the proposed formulae, it is therefore important to test the sensitivity of the results to both the choice of \({\widehat{U}}_{i}\) and to changes in patterns.

As a next step, we must recognize also that the BF method is usually used together with the chain-ladder method. Therefore, as for the hybrid chain-ladder method (see Arbenz and Salzmann 2014), a skewness formula for the mixed BF/chain-ladder method should be developed in a further paper.