{'fraction_non_alphanumeric': 0.03537406710757842, 'fraction_numerical': 0.012418340311065267, 'mean_word_length': 4.459202826543005, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model error and its estimation, with particular application to loss reserving

This paper is concerned with forecast error, particularly in relation to loss reserving. This is generally regarded as consisting of three components, namely parameter, process and model errors. The first two of these components, and their estimation, are well understood, but less so model error. Model error itself is considered in two parts: one part that is capable of estimation from past data (internal model error), and another part that is not (external model error). Attention is focused here on internal model error. Estimation of this error component is approached by means of Bayesian model averaging, using the Bayesian interpretation of the LASSO. This is used to generate a set of admissible models, each with its prior probability and the likelihood of observed data. A posterior on the model set, conditional on the data, results, and an estimate of model error (contained in a loss reserve) is obtained as the variance of the loss reserve according to this posterior. The population of models entering materially into the support of the posterior may turn out to be thinner than desired, and bootstrapping of the LASSO is used to gain bulk. This provides the bonus of an estimate of parameter error also. It turns out that the estimates of parameter and model errors are entangled, and dissociation of them is at least difficult, and possibly not even meaningful. These matters are discussed. The majority of the discussion applies to forecasting generally, but numerical illustration of the concepts is given in relation to insurance data and the problem of insurance loss reserving.


Introduction
This paper is concerned with forecast error, particularly in relation to loss reserving. The majority of the concepts are fairly general, and might be applied in many forecasting environments, but the illustrations given here relate to loss reserving. To an extent, the paper is a sequel to McGuire, Taylor and Miller (2021), which dealt with estimation of a loss reserve by means of the LASSO. Both papers make applications of the LASSO, the first paper to point estimation, and the sequel to estimation of forecast error, particularly model error.
When a forecast is made on the basis of a model of observations, it will inevitably contain an error relative to the true value yet to be observed. Estimation of the properties of this error will provide some understanding of the reliability of the forecast. This has been an issue in the loss reserving literature since raised by Reid (1978), De Jong and Zehnwirth (1983) and Taylor and Ashe (1983).
In subsequent years, forecast error has been decomposed into a number of components, most notably parameter, process and model errors (see e.g. Taylor (1988), O'Dowd, Smith and Hardy (2005), Taylor and McGuire (2016), Taylor (2021)). Estimation of the first two of these three has become well understood, but there has been little development of the estimation of model error.
A notable exception was O'Dowd, Smith and Hardy (2005) and Risk Margins Task Force (2008), who provided a framework for the estimation of each component of forecast error using scorecards to score subjectively a range of factors identified as likely to influence the quantum of forecast error.
Where objective estimates are concerned, Taylor (2021) investigated model distribution error, a component of model error and Bignozzi and Tsanakas (2015) consider model risk in the context of VaR estimation, which is of course relevant to loss reserving, but loss reserving models as such are beyond their scope. Blanchet, Lam, Tang and Yuan (2019) estimate the effect of model error on a performance statistic in terms of the extrema of that measure over the set of admissible models. However, to the authors' knowledge, there has been no other progress in the actuarial loss reserving literature.
In the meantime, the subject has been addressed elsewhere in the economics and finance literature. Useful general overviews are given by Glasserman and Xu (2014) and Schneider and Schweizer (2015) in the context of financial risk management. The approach of Blanchet, Lam, Tang and Yuan (2019) is similar to the latter. Huang, Lam and Tang (2021) estimate the total of parameter error and model error (they call these data variability and procedural variability respectively) in relation to deep neural networks. Their approach equates more or less to bootstrapping, but where the replications are obtained by random variation of the network initialization rather than data re-sampling.
An example of this approach appears in the econometric literature in Loaiza-Maya, Martin and Frazier (2021), whose focus Bayesian prediction is similar to the Bayesian approach followed in the present paper, but with conditional likelihood replaced by a scoring rule. Martin, Loaiza-Maya, Maneesoonthorn, Frazier and Ramírez-Hassan (2021) give a non-Bayesian presentation of the same ideas. This literature has been highly valuable at the fundamental conceptual level. However, the concepts are not easy to operationalize, and the much of the literature does not address loss reserving specifically.
The present paper endeavours to fill some of the literature gaps identified above. It uses the LASSO (Hastie, Tibshirani and Friedman, 2009) to populate abstract concepts, such as model set and its prior distribution, that occur in the more theoretical literature. Since the LASSO may also be used as the source of a loss reserving model, this creates a direct nexus between that model and the estimation of its model error.
The result is that, using the procedures described here, one may perform the following entire sequence of operations: • model a data set of claim observations and extract a point estimate of loss reserve; • move on to estimating the distributions of several components of forecast error that constitute a major part of the total; • supplement these with the distributions of the missing components, derived from other sources; • apply these distributions to the calculation of loss reserve risk margins, or any other quantities of interest that depend on the distribution of forecast error.
The paper is arranged as follows. After the establishment of the necessary mathematical framework and notation in Section 2, the structure of the problem to be considered is established with a review of the components of forecast error, and thereafter the paper focuses on one particular component, internal model structure error ("IMSE") (Section 3). Section 4 discusses the estimation of IMSE in the abstract, and the ingredients required for it. Then Section 5 puts these concepts to work in the specific context of the LASSO. This derives an estimate of the distribution of IMSE, but this estimate is then strengthened by bootstrapping the LASSO in Section 6. The whole procedure is then applied to several synthetic data sets, of varying complexity, in the derivation of numerical results in Section 7. Finally, Section 8 summarizes and considers the successes and limitations of the paper in the attainment of its objectives.

Reserving framework and notation
As far as possible, the notation here will follow that of McGuire, Taylor and Miller (2021). Accordingly, the analysis below will be concerned with the conventional claim triangle. Some random variable of interest is labelled by accident period = 1,2, … , and development period = 1,2, … , − + 1. In this setup, a cell of the triangle refers to the combination ( , ), and the observation in this cell denoted . The payment period to which cell ( , ) relates will be denoted by = + − 1.
Let Δ denote the collection of cells (i.e. ordered pairs ( , )) of which the triangle consists, and let denote the observations on these cells, i.e. = { : ( , ) ∈ Δ}. Accident and development periods will be assumed of equal duration, but not necessarily years. As further notation, [ ] = , [ ] = 2 . A realization of will be denoted .
Let be any random vector defined on Δ and its realization. It will sometimes be useful to vectorize these quantities, and so will denote the column vector of all listed in some defined order, and the corresponding vector of all . This paper will be concerned with forecasts produced by models fitted to the data . Forecasts are made in respect of the for ( , ) ∈ Δ * , some set of cells, disjoint from Δ, and relating to the future, i.e. + > + 1. These will now be denoted * , and the vector of these will be * .
The forecast of * by model will be denoted ̂ * ( ) and is equal to ( ; , ) for some realvalued function . As the notation indicates the forecast is -dependent but, as this notation is cumbersome, the will be suppressed and the forecast written as simply ̂ * when this does not create any ambiguity. Other quantities derived from ̂ * (e.g. immediately below) will also be notated without explicit mention of .
A model will include a likelihood function (. | ) for the observations so that the likelihood of the data vector is ( | ). Let ℓ( | ) = ( | ), the log-likelihood function.
It will be convenient to denote a value fitted by the model to a past observation by ̂= ( ; , ). According to the vector notation given above, ̂ denotes the vector of ̂ for ( , ) ∈ Δ. Similarly, let ,̂ * denote the vector of ̂ for ( , ) ∈ Δ * .
The modeller may select the model but, in practice, will rarely know whether it is a correct representation of the data. Therefore, assume that is selected from some collection ℳ of candidate models, hereafter called the model set. Assume further that the model set is equipped with a measure .
Suppose that = ( * ), for some real-valued function , and define its forecast as ̂= (̂ * ). An example is = 1 * , where 1 is a vector of the same dimension as * and with all components equal to unity. If the * denote claim payments, then is the amount of outstanding claim liability.
The forecast error associated with forecast ̂ will be defined as = −̂.
(2.1) Later sections will make use of open-ended ramp functions. These are single-knot linear splines with zero gradient in the left-hand segment and unit gradient in the right-hand segment. In a machine learning context, one of these would be referred to as a rectified linear unit ("ReLu"). Let ( ) denote the open-ended ramp function with knot at . Then For a given condition , define the indicator function ( ) = 1 when is true, and ( ) = 0 when is false. Further, define the Heaviside function ( ) = ( ≥ ).

Components of forecast error
Some regulatory regimes require that a capital margin be associated with a technical reserve such that the total of reserve plus margin equal at least the Value at Risk ("VaR") of the liability at some high percentile, such as 99.5%. These regimes may also require the evaluation of a risk margin within the capital margin. For example, the Australian prudential standards require a loss reserve to be at least equal to the 75% VaR of the associated liability (Australian Prudential Regulatory Authority, 2018).
IFRS17 requires a risk adjustment for non-financial risk. While there is no prescribed method for calculation, margins based on the Value at Risk ("VaR") are likely to be widely used.
For the sake of definiteness, this paper proceeds on the basis that the specific liability forecast under consideration is the loss reserve. For some regulatory regimes (e.g. IFRS17) it might be the one-year-ahead forecast of claim payments, with an associated VaR. In such cases, the model error methodology set out in subsequent sections translates readily to this alternative situation.
The requirement of a VaR necessitates the estimation of the distribution of a forecast liability, as opposed to a simple point estimate. Often, especially in the case of low-or mediumpercentile VaRs, it is reasonable to assume that the required distribution is characterized by its mean (the point estimate) and variance. Hence the need for examination of the variance of forecast error.
It consists of a number of identifiable components. Decomposition of forecast error is discussed by Taylor (2000), Taylor  The different authors use slightly different terminology, and so  The most comprehensive decomposition is given by Taylor (2021), on which the following is largely based. Henceforth, let the model introduced in Section 2 consist of just its distribution-free part, i.e. the algebraic form of its non-stochastic part (e.g. the linear predictor in the case of a GLM), and let denote the model distribution of the stochastic error. With a slight abuse of notation, the quantity [ ] will also be referred to as forecast error when this usage is unambiguous. The four terms on the right may be recognised as respectively model bias, process error, parameter error ("PaE") and model error.
The model error may be further decomposed with a small amount of manipulation: The two members on the right are labelled model structure error and model distribution error in Taylor (2021). If the latter is assumed away by assuming the distribution of to be concentrated in a single point, then O'Dowd, Smith and Hardy (2005) recognize that some of the variation over models ∈ ℳ relates to mis-specification of the model fitted to past data, and part to the extrapolation of that model in its forecasts.
In the case of a GLM, for example, the data will have been modelled using some linear response , and the data set augmented by the forecast experience will use an augmented linear response + + , where (3.5) and the parameter sub-vector * is not estimated from data, but consists of assumptions about the future, which are combined with the design matrix subcomponent * * to form the augmented part of the linear predictor. For example, superimposed inflation ("SI"), a diagonal effect, might be estimated at one rate in past data, but a different rate assumed for the future. The past rate would be included in , and the future in * .
The forecast of future claims experience depends on the linear response * + * * * , in which the first of the two members relates to past observations, and the second to assumptions about the future. Forecast errors in relation to these two are referred to as internal and external model errors respectively.
Hence, the first of the two members on the right side of (3.3) can be decomposed further: where , are those parts of model generating internal and external model errors.

Figure 3-1 Full decomposition of model error
The estimation of parameter and process error is covered by England and Verrall (1999, 2001, 2002, Taylor (2000) and Taylor and McGuire (2016). Model distribution error is discussed by Taylor (2021 The present paper will be concerned with model structure error, and specifically IMSE. As explained earlier in this section, this relates to any error in the selection of the algebraic form fitted to the data, and the objective here is to estimate this error on the basis of the data. External model error, on the other hand, relates to future events and influences on future claim data and, by definition, there are no available data relevant to these. The estimation of this component of error must be performed by some means other than reference to past claim data As mentioned in Section 2, one contemplates models ∈ ℳ, the model set, with probability measure π defined on ℳ. Each model defines a loss reserve estimate ̂( ) = 1̂ * ( ).
With this machinery, one is in a position to calculate summary statistics of the estimated loss reserve over candidate models. Of particular interest are Quantity (4.2) measures the dispersion of estimated loss reserve over the entire model set, and will be of interest in the formulation of model error in Section 5.2.

Model confidence set
This is discussed by Hansen, Lunde and Nason (2011), Glasserman and Xu (2014) and Schneider and Schweizer (2015). In its most basic definition, a subset ⊂ ℳ is a 100 % model confidence set relative to some reference model 0 if 0 ∈ and π( ) = . This is the approach of Hansen, Lunde and Nason (2011), but there are other possibilities. For example, Schneider and Schweizer (2015) are concerned with identification of the worst-case risk within a prescribed divergence radius of the reference model.
This approach uses the structure of Section 4.1 but with π representing a Bayesian posterior distribution on ℳ. It is supposed that a prior distribution function π is specified on ℳ.
Then application of Bayes theorem yields the posterior distribution Consider the case in which ̂( ) is any forecast of model , not necessarily a loss reserve. Its posterior mean is which is a weighted average of the forecasts across all models, with weights π ( ). The concept of averaging across models was introduced to the loss reserving context by Taylor (1985). Thus, it is supposed that the mean [ ] is defined by the parameters , and the distribution around this mean defined by the parameters .

LASSO estimation of internal
A LASSO regression estimates , of dimension , as where the dependence of ̂ on has been made explicit, is some measure of the separation between and ( ), the are the components of , and the ≥ 0, = 1, … , are fixed values selected by the modeller.
A choice of such that the regression reduces to a GLM in the event that all = 0, is the negative log-likelihood function ( ,( )) = −ℓ( | ), in which case (5.3) becomes where the model has been represented by its parameter vector .
There are a couple of fundamental derivations of the estimator (5.3), specifically as (Hastie, Tibshirani and Friedman, 2009): (a) the minimizer of the quantity subject to an upper bound on each | |; and (b) the maximizer of the a posteriori estimate of the quantity when each coefficient is random subject to a Laplace density.
Case (b) will be of interest here, and so is explained in a little more detail. It is assumed that is distributed with Laplace probability density where is a scale parameter. This is a symmetric density about a zero mean, and Then the posterior likelihood of is proportional to where the modulus operates on a vector element-wise and denotes the column -vector with the as components. This posterior likelihood may be compared with (5.4).

Implementation
The application of the LASSO to data follows the procedure set out in McGuire, Taylor and Miller (2021). The target variable denotes the observed claim payments over a 40 × 40 triangle. The model takes the same form as a GLM with log link, specifically supplementing (5.1) and (5.2)  In this formulation, the parameter from (5.1) is set equal to the , the over-dispersion parameter. This parameter requires estimation, which is discussed later in this sub-section.
A LASSO can be implemented by means of the R package glmnet but only with Poisson ( = 1) rather than Gamma distribution. The solution is to do just this but, later, in the BMA of Section 5.2, adjust all Poisson distributions to Gamma with suitable . A Gamma error term is preferred as more realistic than Poisson but, although the package glmnet v4.0 purports to provide the means for this, numerical difficulties precluded its application here.
All covariates were standardized, as described by McGuire, Taylor and Miller (2021). The orders of such basis sets are typically in the thousands or tens of thousands, depending on the dimensions of the data being used, but the LASSO will eliminate most of these from the model.
The vector is restricted to the form (0, … ,0, ′ , ′ , … , ′ ) , where ′ > 0 and the leading zeros relate to specific data features that the modeller regards as of certain existence and wishes to force into the modeltypically an intercept term is fitted without penalty. Thus, the terms of the linear response corresponding to these covariates are always included without penalty, and all other members of the linear response carry penalty parameters that are positive and equal. Details appear in McGuire, Taylor and Miller (2021).
The fit of each to data is assessed by 8-fold cross validation using the Poisson deviance as the loss function. The cross-validation loss in the -th fold is denoted by . The average of the values across is denoted by ̅̅̅̅ and the standard error of the same values by .
At this point, the parameter is estimated. The model structure (i.e. the covariates with nonzero coefficients) of the 1se model is extracted, and a GLM fitted to this structure. Henceforth, the preferred Gamma error is assumed, so the GLM is fitted using the R function glm() from the stats package with the same weights as in the LASSO fit. Since only a crude estimate of is yielded by glm(), the maximum likelihood estimate of is obtained from the gamma.dispersion() function from the MASS package.
This last value of is used in conjunction with an assumed Gamma error for all calculations henceforth. It is not re-estimated in the bootstrap replications of Section 6. The GLM itself is of no particular interest other than the estimation of .
The 1se and minCV models generate the LASSO estimates of outstanding claim liability chosen by McGuire, Taylor and Miller (2021). The extrapolation of a model of past data to the future require decisions as to the extent to which any trends identified in the past are also extrapolated, and to what extent.
For this purpose, the decisions of McGuire, Taylor and Miller (2021) have been followed, in that all past trends have been extended into the future without modification. In more precise terms, if the linear response included in the model of past data is ( , , ) (consistent of a linear combination of ramp functions and products of Heaviside functions), then the linear response for the future will also be taken as simply ( , , ).
This has the benefit of even-handedness in the generation of a model set. There is no user intervention in the exclusion of models that fail to conform with preconceptions. On the other hand, however, it does admit some models that might be viewed as dubious. For example, any payment quarter effect that is modelled to take effect in the last few payment quarters of experience will be extended indefinitely into the future. Ultimately, such exotic models may need to be pruned away (Section 5.2.5).

Model bias
The LASSO's shrinkage of parameter estimates induces bias in model estimates. For this reason, the LASSO is sometimes used just for selection of the model form, and then the selected model refitted to the data without penalty, i.e. as a GLM. The result should be approximately unbiased.
McGuire, Taylor and Miller (2021) carried out some experimentation with this procedure in relation to the example data sets of Section 7.1, but commented that the results were not encouraging in that the refit caused a substantial deterioration in model performance in some cases.
See the further comment in Section 5.2.1. 5.2. Bayesian model averaging of LASSO 5.2.1. Generation of model set The cell mean structure within the LASSO model is given by (5.2). Variation of there will create a multiplicity of models ℳ = { (. ; . , . , ), ∈ ℝ }. It is of particular interest that, as increases to a certain threshold, the LASSO forces to zero for that and all greater values. This amounts to dropping a covariate from the model, equivalent to changing the cell mean's algebraic structure. This is exactly one of the requirements of Section 4.1. In fact, the LASSO provides a rich model set that includes both variation in model algebraic structure and parametric variation within each distinct structure. The latter form of variation arises from variation in the penalty parameter, and is lost if a GLM is fitted to a LASSO model structure.

Model likelihood
As explained in Section 5.1.2, the process of BMA about to be described assumes a Gamma distribution for each observation . The log-likelihood is where = 1⁄ , = 1⁄ , with the given by (5.8) and the estimation of is described in Section 5.1.2.
In the case of (5.8), consists of columns representing all functions in the basis set defined in Section 5.1.2. In the row corresponding to cell ( , ), each of the functions is evaluated for those values of , . 5.2.3. Prior distributions of LASSO models By (5.5), assuming the random parameters are independent, and this may be viewed as a prior distribution on the models ∈ ℳ.
As appears in Section 5.1.2, the vector is defined by the scalar ′ . It is of especial relevance that, in (5.10), ′ functions as a dispersion parameter for the prior. There is a correspondence ⇔ ′ ( ) , so that models may be labelled equivalently by or ′ . As noted in Section 5.1.2, models are numbered in such a way that ′ ( ) ↓ as ↑. Decreasing ′ ( ) implies lower penalty on parameters, and so model complexity increases steadily with increasing . In the limit ′ ( ) = 0, the LASSO degenerates to a GLM with the set of covariates equal to the entire basis set, i.e. an extremely complex model.
Similarly, high values of ′ ( ) force the model toward simplicity. As ′ ( ) → ∞, the model tends to one containing only the unpenalized covariates.
Each value of ′ ( ) also induces a posterior distribution π ( | ; ′ ( ) ), where the dependence on appears twice; first, the algebraic structure of the model depends on and, second, the posterior is influenced by the prior dispersion parameter. This sets up a sequence of relations ⇔ ′ ( ) ⇒ π ( ) , where the last term is written as an abbreviation for π ( | ; ′ ( ) ).
Recall, however, that all of this applied in the context of Poisson likelihood (Section 5. ( ) may be found numerically as that value that, for any given , produces the required "equality" between π ( ) and π [ ] ( ) . This course has been followed here. (1 ) respectively. Subsequent sections will rely mainly on these as relating to reasonable priors (equivalently, selections of ′ [ ] , that may not be specific to any LASSO model ). Moreover, as each of these models will usually exhibit a high degree of compatibility with the data, the posterior distribution that follows from its prior will usually be centred somewhere in the vicinity of the posterior mode. In fact, the 1se model will be defined as the primary model, meaning that the actuary is assumed to use this as the source of an adopted loss reserve.
However, it will also be useful to investigate the extent to which ′ [ ] may credibly deviate from these selections. Therefore, two values of this dispersion parameter are identified as of special interest. These may be labelled "simple-model" and "complex-model" dispersion parameters, and are defined as follows: for some suitably small value > 0.
) is the prior dispersion parameter which, when applied to all models , is just high enough to force the aggregate prior probability of all models 1 and more complex below a threshold probability. This is the value of ′ [ ] that, when applied across all , generates models so simple that are even remotely consistent with the data and model 1 .  (1 ) would take the prior distribution associated with LASSO model 1 , and assume it to apply to all models in ℱ. When normalized across ℱ, this gives a prior distribution on ℱ, thus where is the forecast of given by model , as parameterized by the LASSO.
Though it may be reasonable, as suggested above, to select the prior 1 ( 1 ) for use in (5.13), it is of interest to consider other possibilities. Therefore, re-write (5.13) in the more general form: where 0 ( ( ) ) is constructed as for 1 ( ( ) ) except for the choice = [ ] (0) , which is some preferred but unspecified prior.
BMA will be carried out below under various selections of prior, namely (in order of increasing model complexity): The middle two of these would be regarded as yielding realistic results. The other two indicate the extremes to which results may be taken before consistency with data is totally lost.

Pruning of model set
Loss reserve forecasting is essentially an exercise in extrapolation. Certain functions are fitted to past data, and then extrapolated to forecast future data.
In the highly supervised environment of conventional loss reserving, this rarely creates difficulty. The actuary maintains strict control over the explanatory functions, often assigning physical explanations to them. Supervision of extrapolation behaviour is part of this control; if it appears unreasonable, the actuary will modify the model.
In the case of the LASSO, and other machine learning methods, on the other hand, the volume of models under test precludes this level of supervision. Moreover, the functions used to fit past data are assembled abstractly from a set of basis functions. Interpretation of the results may or may not be elementary.
It is in the nature of such exercises that strict control may be maintained over the fit of models to the data (e.g. by cross-validation error, as in Section 5.1.2), yet there be no control over the behaviour of these models in future. Hence, good quality fit of a model to data provides no assurance of reasonable forecast behaviour.

Figure 5-2 Example of high-quality fit to data but poor extrapolation
It can be seen that the candidate model coincides with the primary almost perfectly over the past payment quarters, and its posterior probability is in fact 0.67. However, its forecast quickly diverges from the primary's over future payment quarters. Ultimately, it forecasts future liability (for the last 10 accident quarters) 56% higher than the corresponding primary forecast.
It is highly unlikely that an actuary, having invested in a primary model, would regard a model that diverges to this degree as suitable for inclusion in the model set. If surveillance over inclusions were maintained, this particular case would be censored. For practical purposes, this censorship must be automated.  It is found that the dominant mode of censorship relates to breaches of the accident, rather than payment, quarter limits. This seems unsurprising as the limited experience of recent accident years provides little guide as to future development. It is easy for extrapolation functions to go awry.
The inclusion gates must be subjectively chosen, since they represent the actuary's judgement on the maximum acceptable deviation from the primary model. The gates are likely to vary from one reserving exercise to another, as the line of business, length of development tail, and other particulars change.
It is noteworthy that this introduces an unwelcome element of subjectivity into an otherwise data-driven estimation process. Unfortunately, we do not see an alternative. Failure to address the issue causes delinquent models to lead to unduly large estimates of internal model error that lack credibility.
5.2.6. Internal model structure error from the primary model Equations (4.1) and (4.2) give expressions for the mean and variance of a posterior distribution of loss reserve. Equation (4.4) is the counterpart of (4.1) when the posterior derives from BMA across LASSO models. It will often be convenient to convert the posterior variance to a coefficient of variation ("CoV").
In this context, the posterior mean may be interpreted as an estimate of the loss reserve, and the CoV as a measure of IMSE. Numerical examples of this are presented in Section 7.2.
6. Bootstrapping the LASSO 6.1. Motivation Section 5.2.6 described the estimation of IMSE for the primary model. In this sense, the problem is solved at that point. However, the estimate may be regarded as a little "thin" since the number of models with posterior probability greater than vanishingly small may be limited. It would not be unusual for this number to be in the range 10-30. This is likely to lead to underestimation of the IMSE.
Bootstrapping of the LASSO generates multiple pseudo-data sets, leading to multiple posterior distributions, and adding to the reliability of the IMSE estimate.

Form of bootstrap
A semi-parametric form of bootstrap is used. This terminology is taken from Shibata (1997), and its implementation in relation to loss reserving (in relation to a GLM in that case) is explained in Taylor and McGuire (2016).
Pseudo-data sets are formed on the basis of the LASSO standardized residuals. Standardization requires an estimate of variance for each cell. In accordance with the assumption of a Gamma distribution in Section 5.1.2, the variance of the ( , ) cell is estimated as ̂̂2 , where ̂ is the value fitted to the cell by the LASSO, and ̂ is the estimate of obtained as described at the end of Section 5.1.2.
The LASSO is known to be biased, and so the set of standardized residuals must be centred before use in the bootstrap. They are all shifted by a constant quantity such that their arithmetic average is zero. One needs to check that the resulting centred standardized residuals appear iid, and they are then re-sampled to generate multiple pseudo-data sets [ ]  Each pseudo-data set is then submitted to the entire process described in Section 5, i.e. LASSO estimation, followed by BMA, with some modifications to censor unreasonable models and to centre the bootstrapped results to the original forecasts for each prior. Centring the bootstraps leads to a change in posterior probabilities so the process followed is: • Carry out LASSO estimation to produce models . Note that it may be advisable to use a less onerous level of censoring at this level than set out in Table 5-1 to avoid eliminating models that may be acceptable after centring (next point). Hence temporary inclusion gates are obtained by scaling those in the table by a factor of 1.4. • Calculate the scaling factor required to scale the sample mean of the posterior means from each pseudo data set to the original posterior means. • Apply censoring at the desired level as discussed in Section 5.2.5.
• Drop any bootstrap replication for which posterior distribution over models [ ] is concentrated on less than 5 models in the sense of having fewer than 5 models with posterior probability masses each exceeding 0.0001.
The subsequent BMA is dealt with in Section 6.3. . Some of the rows of the matrix will be empty when all models of those rows have been censored. In machine learning parlance, the bootstrap matrix might be considered an ensemble, or a ragged array.

The bootstrap matrix
Note that, for fixed , the models [ ] vary over , i.e. a given column of the matrix might not contain the same model in all rows. Indeed, the models might be different in all rows of that column.
BMA produces an estimated liability Conventional use of the bootstrap in loss reserving enables estimation of PaE. This is achieved by means of a well-established procedure (England and Verrall, 1999;Taylor and McGuire, 2016). PaE is estimated from the variation in loss reserve estimates over bootstrap replications. Process error can be added quite separately, and does not require the bootstrap.
A parallel procedure can be followed in relation to the bootstrap matrix here. One can calculate and this is an estimate of parameter error. Once again, the estimation of process error CoV, denoted , is separate and straightforward, using sampling from its known distribution (Section 6.2), as described by England and Verrall (1999, 2001, 2002, Taylor (2000) and Taylor and McGuire (2016) .
Note that, by the law of total variance,
The combination of PaE and IMSE has squared CoV + 2 = 2 + 2 2 = 2 + 2 . (6.6) Remark 6.1. The thinness of the primary estimate of IMSE, and the possibility of underestimation, was noted in Section 6.1. Each bootstrap replication has the same stochastic properties as the primary, and so each row of the bootstrap matrix, and also the average of them, may lead to under-estimation. ∎ Remark 6.2. In view of Remark 6.1, the meanings of IMSE and PaE, as estimated by (6.2) and (6.4), merit some further consideration. The meaning of the latter is quite clear in the case of a GLM. Here there is a fixed model structure across all bootstrap replications, and it is applied to the various data sets generated by those replications. Only the parameterization of that structure varies across replications. Quite clearly, this leads to genuine estimation of PaE.
When the LASSO is bootstrapped, the situation is quite different. Variation of the data set from one replication to another, induces different sets of models. Hence, not only parameterization of models, but the models themselves, vary across replications. In this way, parameter and model errors become entangled, and the distinction between them becomes unclear. Indeed, there may even be a question as to whether the distinction is even meaningful.
From a practical viewpoint, it may be preferable to think in terms of just (6.6) as unambiguously the combination of parameter and IMSE errors without concern for the decomposition into the two components. For most practical purposes (e.g. risk margins), only the combination will be required.
The estimates and , as defined by (6.2) and (6.4), are useful for the construction of (6.6), but no particular meaning may be assigned to the components themselves. ∎ 7. Numerical results 7.1. Data sets All data submitted to analysis in this paper are synthetic. Thus, all data sets contain known features. Four synthetic data sets were constructed and analyzed. Each consisted of a 40 × 40 quarterly triangle of incremental paid claims ; thus 810 cells. The data were simulated according to exactly the same specifications as set out in McGuire, Taylor and Miller (2021).
Full details were given in that earlier paper. It suffices here to give just the briefest descriptions of the data sets, as below.
Data set 1: Cell means consist of multiplicative row and column effects, and so this data set is a simple one, consistent with a chain ladder model. Data set 2: Cell means are as for Data set 1, except that a payment quarter effect has been added. This effect represents SI at a rate that is constant over the first 12 quarters, increases over the next 12, returns to constant over the next 8, then increases once more over the final 8. For any one payment quarter, the rate of SI is constant across development quarters.
Data set 3: This is as for Data set 2, except that one further complexity has been added. This consists of a sharp increase in claim costs in development quarters 21 and later within accident quarters 17 and later. This last feature is of particularly difficulty for detection by a model, as it affects only 10 of the 810 cells of which the data set is composed.
Data set 4: This is as for Data set 2 except that greater complexity has been added to the rates of SI. The rates from Data set 2 now apply to development quarter 1 and, for any payment quarter, they decline linearly to zero at development quarter 40.
These data sets might reasonably be listed 1, 2, 3, 4 in increasing order of complexity, though estimation might be more difficult in case 3 than case 4.

Primary model
Each of the four data sets of Section 7.1 are modelled and extrapolated to produce an estimated loss reserve and IMSE on the basis of the primary model, as described in Section 5. The results are set out in Table 7-1. Since the data are synthetic, the true reserves are known, and these have also been included in the table.
The table includes two sets of forecasts. The "raw 1se" is the 1se estimate from the primary model, obtained by a single application of the LASSO to the data set. The "posterior" is the mean of the posterior distribution obtained in the BMA described in Section 5.1.2.
The "1se" and "minCV" rows are bolded as the most relevant. The "simple and "complex" rows are also included to indicate the extent to which IMSE may vary before the models taken into account by it exhibit lack of fidelity to the data.
It may be noted that the posterior means usually exceed the true liabilities. This is not a product of the BMA, as the difference between the raw 1se estimate and its Bayesian average is small in all cases. The 1se models simply over-estimate liability in these cases. The modelling of Data set 1 behaves well, with small estimates of model error, as is to be expected since the data is compatible with a simple chain ladder model, reflecting just row and column effects.
It is to be emphasized that, although Table 7-1 focuses attention on the CoV of IMSE, the full posterior distribution of loss reserves is available for each row of the table. As an example, Figure 7-1 displays the posterior corresponding to the 1se prior in the case of Data set 2. The convolution of this with the distributions of other components can be taken to obtain the distribution of full forecast error.

Bootstrap
The bootstrap procedure described in Section 6 was applied to all four data sets. In each case, the number of bootstrap replications was set to = 200, although, as will be seen in the results below, some of these are discarded by the pruning regime of Section 5.2.5. Table 7-2 displays the results corresponding to the primary model results in Table 7-1.  Table 7-2 with Table 7-1 reveals the effect of bootstrapping in supplementing the primary model. The magnitudes of these estimates of IMSE vary over the four data sets somewhat as expected; at least, the largest estimate occurs for the exotic Data set 3. It is interesting that there is little difference between the estimates for Data sets 1, 2 and 4, even though the complexity varies between them. This suggests that the LASSO models perform about equally well in these three cases.
The final column illustrates how the censorship rates increase with increasing complexity of the data set. In other words, as models are confronted with increasing data complexity, those models include additional feature that render their extrapolation more prone to delinquency. This is particularly the case for Data set 3, for which the forecasting difficulties have been discussed in Section 7.1.
As has been discussed in Section 6.3, the LASSO bootstrap provides an estimate of PaE in addition to ISME, and process error may be estimated routinely. These are displayed in Table  7-3. Only the cases of 1se and minCV priors are included as those of greatest interest.
As noted in Remark 6.2, there is some ambiguity in the separate meanings of PaE and IMSE, but not in the meaning of their combined effects. Table 7-3 therefore also contains a column displaying the CoV associated with the sum of these two components of forecast error.
The "Sub-total" error is the sum of parameter, process and IMSE. It is designated in this way to emphasize that it is not, in fact, the total forecast error because it omits the two components not discussed in this paper, namely • External model structure error; and • Model distribution error. It may be noted that the magnitude of forecast "sub-total" error usually depends little on whether it is derived from a 1se or a minCV prior distribution. There is one exception to this, Data set 4, where the difference relates mainly to estimated "parameter error".
Investigation reveals that the minCV distribution of posterior means across bootstrap replications exhibits a markedly longer right tail than its 1se counterpart. The reasons for this have not been investigated further, but note the complex structure of SI for this data set. Misestimation of SI as a result of undue model complexity could lead to volatile forecast along future diagonals.
There is some support for this explanation in Table 7-4, which repeats the parameter error CoVs from Table 7-3, but now flanked by those resulting from the simple and complex priors. The table displays a clear divide between the two simpler and the two more complex cases. 7.4. Sensitivity to inclusion gate width All numerical results reported above are governed by the inclusion gates set out in Table 5-1. Since these gates were subjectively determined, it is natural to enquire into the sensitivity of the results to the gate-widths. Table 7-5 compares the "Sub-total" forecast errors of Table 7-3 with the errors that are obtained when all gates are widened by increasing the upper limit by a factor of 1.1, and decreasing the lower limit by the same factor. Thus the first gate in

Sub-total error (CoV)
It is seen that, in very broad terms, widening the inclusion gates by a factor of 1.1 increases forecast error by a factor of about 1.5. Thus, estimated forecast error can be sensitive to the selection of these gates, and they must be selected carefully.

Reasonableness check of results
It is interesting to use a different methodology in an attempt to provide an independent check on the reasonableness of the above estimates of PaE plus IMSE. For this purpose, a GLM has been fitted to each data set, treating the model structure as known, i.e. containing the same explanatory variables as in the simulation specifications. Thus, the GLMs should have low levels of IMSE. The GLMs have then been bootstrapped to yield estimates of PaE. The GLM bootstrap forecasts were subject to the same set of inclusion gates as applied to the LASSO (see Table 5-1).
Before examining the results of these models, it is useful to consider reasonable expectations of them. As noted in Remark 6.2, some IMSE may leak into recorded PaE in the case of the LASSO. Accordingly, it would not be surprising if the PaE, as estimated by the LASSO, exceeded that estimated by the GLM.  Table 7-6 provides the comparison. It is certainly the case that the LASSO 1se estimates of PaE exceed the GLM estimates. Specific observations are that: • the very simple Data set 1 and the more complex Data set 2 both produce substantial differences; • the discrepancies for Data sets 3 and 4 are also significant, but of lower magnitude.
These results appear consistent with the expectations set earlier in this sub-section.

Conclusions
Section 1 noted that this paper set out to study model error. This is all about forecasting under model uncertainty. Sections 5 and 6 have indeed set out a LASSO-based methodology for the estimation of a specific component of model error, namely internal model structure error, and the methodology has been illustrated numerically in Section 7.
With one exception, noted in Section 5.2.5, the estimation of IMSE outlined here is rigorous, based on Bayesian model averaging over a space of admissible models. This provides an alternative to the subjective estimation found in prior literature.
The exception discussed in Section 5.2.5 is concerned with the fact that some models that fit a data set extremely well may yet forecast abominably. In such cases, there is no apparent datadriven means of excluding the delinquent models. This matter has been dealt with by the establishment of subjectively chosen inclusion gates. Models are rejected as unrealistic if they produce certain results that fail to pass through these gates.
The selection of these gates would be based on the reserving actuary's professional expertise. In one sense, the need to rely on this element of subjectivity is unfortunate. However, the issue cannot be disregarded, as delinquent models have the potential to corrupt the model error estimation completely, and produce wildly inflated estimates.
The authors see no alternative to some form of subjective decision on the admissibility of individual models. There are, however, obvious risks associated with it in that under-(over-) estimation of model error will result from unduly strict (lax) selection of the gates. Moreover, the selection requires considerable care as the estimate of model error is materially sensitive to it (Section 7.4).
The suggested estimation procedure is supplemented by bootstrapping in Section 6. This yields a bonus in that it produces an estimate of PaE as well as IMSE. However, as noted in Remark 6.2, the supposed estimates of PaE and IMSE are entangled, probably inextricably. Possibly, there is not even any clear meaning to be associated with these separate components.
But all is not lost here. Although there are not reliable estimates of PaE and IMSE (whatever they may be), Remark 6.2 points out that the estimate of their combined contribution to forecast error is indeed reliable, and this is all that is required for most practical purposes, such as loss reserve risk margins.
The estimates of PaE and IMSE obtained from the LASSO have been subjected to a sanity check in Section 7.5 by comparison with the PaE estimate obtained from a GLM, and appear generally reasonable.
It is important to note that, as pointed out in Section 3, IMSE is but one component of model error. The estimation of IMSE, PaE and process error are discussed in this paper, and illustrated numerically in Section 7, but this leaves model distribution error and external model structure error still to be brought to account.
The investigations of Taylor (2021) hint that the first of these would usually be small, but the second often would not. The only discussion of this component of forecast error in the actuarial literature occurs in O'Dowd, Smith and Hardy (2005) and Risk Margins Task Force (2008), where a subjective approach is taken. This may well be the only option in making allowance for influences that are totally divorced from past experience and data. Be that as it may, the fact remains that any estimate of external model structure error must be derived by means quite different from those used for IMSE. If a subjective estimate is required, then, once obtained, it may be combined with the other components dealt with in this paper.