{'fraction_non_alphanumeric': 0.05177119228118539, 'fraction_numerical': 0.024424534803583735, 'mean_word_length': 3.8790854068594487, 'pattern_counts': {'":': 1, '<': 1, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Asymptotics for in-sample density forecasting

This paper generalizes recent proposals of density forecasting models and it develops theory for this class of models. In density forecasting, the density of observations is estimated in regions where the density is not observed. Identification of the density in such regions is guaranteed by structural assumptions on the density that allows exact extrapolation. In this paper, the structural assumption is made that the density is a product of one-dimensional functions. The theory is quite general in assuming the shape of the region where the density is observed. Such models naturally arise when the time point of an observation can be written as the sum of two terms (e.g., onset and incubation period of a disease). The developed theory also allows for a multiplicative factor of seasonal effects. Seasonal effects are present in many actuarial, biostatistical, econometric and statistical studies. Smoothing estimators are proposed that are based on backfitting. Full asymptotic theory is derived for them. A practical example from the insurance business is given producing a within year budget of reported insurance claims. A small sample study supports the theoretical results.

This paper generalizes recent proposals of density forecasting models and it develops theory for this class of models. In density forecasting, the density of observations is estimated in regions where the density is not observed. Identification of the density in such regions is guaranteed by structural assumptions on the density that allows exact extrapolation. In this paper, the structural assumption is made that the density is a product of one-dimensional functions. The theory is quite general in assuming the shape of the region where the density is observed. Such models naturally arise when the time point of an observation can be written as the sum of two terms (e.g., onset and incubation period of a disease). The developed theory also allows for a multiplicative factor of seasonal effects. Seasonal effects are present in many actuarial, biostatistical, econometric and statistical studies. Smoothing estimators are proposed that are based on backfitting. Full asymptotic theory is derived for them. A practical example from the insurance business is given producing a within year budget of reported insurance claims. A small sample study supports the theoretical results.
1. Introduction. In-sample density forecasting is in this paper defined as forecasting a structured density in regions where the density is not observed. This is an electronic reprint of the original article published by the Institute of Mathematical Statistics in The Annals of Statistics, 2015, Vol. 43, No. 2, 620-651. This reprint differs from the original in pagination and typographic detail. 1 2

LEE, MAMMEN, NIELSEN AND PARK
This is possible when the density is structured in such a way that all entering components are estimable in-sample. Let us, for example, assume that we have one covariate X representing the start of something; it could be onset of some infection, underwriting of an insurance contract or the reporting of an insurance claim, birth of a new member of a cohort or an employee losing his job in the labour market. Let then Y represent the development or delay to some event from this starting point. It could be incubation period of some disease, development of an insurance claim, age of a cohort member or time spend looking for a new job. Then X + Y is the calendar time of the relevant event. This event is observed if and only if it has already happened until a calendar time, say t 0 . The forecasting exercise is about predicting the density of future events in calendar times after t 0 . The most typical example of a structured density is a simple multiplicative form studied by Mammen, Martínez-Miranda and Nielsen (2015). The multiplicative density model assumes that X and Y are independent with smooth densities f and g. When f and g are estimated by histograms, our in-sample forecasting approach could be formulated via a parametric model. This version of in-sample density forecasting is omnipresent in academic studies as well as in business forecasting; see Martínez-Miranda et al. (2013) for more details and references in insurance and in statistics of cohort models. Extensions of such parametric histogram type of models can often be understood as structured density models modelled via histograms. A structured density is defined as a known function of lower-dimensional unknown underlying functions; see Mammen and Nielsen (2003) for a formal definition of generalised structured models. Under the assumption that the model is true, our forecasts do not extrapolate any parameters or time series into the future. We therefore call our methodology "in-sample density forecasting": a structured density estimator forecasting the future without further assumptions or approximate extrapolations.
Our model is related to deconvolution, but there are two major differences. First, in our model one observes not only X + Y but also the summands X and Y . Second, X and Y are only observed if their sum lies in a certain set, for example, in an interval (0, t 0 ]. This makes X and Y be dependent and the estimation problem be an inverse problem. We will see below that the first difference leads to rates of convergence that coincide with rates for the estimation of one-dimensional functions in the classical nonparametric regression and density settings. The reason is that our model consists in a well-posed inverse problem. In contrast, deconvolution is an ill-posed inverse problem and allows only poorer rates of convergence. This paper adds three new contributions to the literature on in-sample density forecasting. First of all, we define smoothing estimators based on backfitting and we develop a complete asymptotic distribution theory for these estimators. Second, we allow for a general class of regions for which IN-SAMPLE DENSITY FORECASTING 3 the density is observed. The leading example is a triangle. A triangle arises in the above examples where the sum of two covariates is bounded by calendar time. The theoretical discussion in Mammen, Martínez-Miranda and Nielsen (2015) was restricted to this case. But there exist many other important support sets; see, for example, Kuang, Nielsen and Nielsen (2008) for a detailed discussion. Third, we generalize the forecasting model by modelling a seasonal component. This is done by introducing an additional multiplicative seasonal factor into the model. Then we have three one-dimensional density functions that enter the model and that can be estimated in sample. Seasonal effects are omnipresent: onset of some disease could be more likely in the winter than in the summer; new jobs might be less likely during the summer or they may depend on the business cycle; more auto insurance claims are reported during the winter, but they might be bigger on average in the summer; cold winters or hot summers affect mortality. When a study is running over a few years only and one or two of those years are not fully observed, data might be too sparse to leave these two years out of the study. Leaving them in might however generate bias. The inclusion of seasonality in this paper solves this type of problems and allow us in general to do well when years are not fully observed. An illustration producing a within-year budget of insurance claims is given in the application section.
Classical actuarial methodology does not include seasonal effects. Budgets are normally carried out manually by highly paid actuaries. The automatic adjustment of seasonal effects offered by this paper is therefore potentially cost saving. Insurance companies currently use the classical chain ladder technique when forecasting future claims. Classical chain ladder has recently been identified as being the above mentioned multiplicative histogram insample forecasting approach; see Martínez-Miranda et al. (2013). The seasonal adjustment suggested in this paper is therefore directly implementable to working routines and processes used by today's nonlife insurance companies.
Recent updates of classical chain ladder include Nielsen (2009), Verrall, Nielsen andJessen (2010), Martínez-Miranda et al. (2011) and Martínez-Miranda, Nielsen and Verrall (2012). These papers reinterpreted classical chain ladder in modern mathematical statistical terms. The generalised structured nonparametric model of this paper is a multiplicative density with three effects. The third seasonal effect is a function of the covariates of the first two effects. Estimation is carried out by projecting an unstructured local linear density estimator, Nielsen (1999), down on the structure of interest. The seasonal addition to the multiplicative density model of Mammen, Martínez-Miranda and Nielsen (2015) is still a generalised additive structure, a simple special case of generalised structured models. Generalised structured models have historically been more studied 4 LEE, MAMMEN, NIELSEN AND PARK in regression than in density estimation. Future developments of our insample density approach will therefore naturally be related to fundamental regression models; see Linton and Nielsen (1995), Nielsen and Linton (1998), Opsomer and Ruppert (1997), Mammen, Linton and Nielsen (1999), Jiang, Fan and Fan (2010), Park (2005, 2006), Nielsen and Sperlich (2005), Mammen and Nielsen (2003), Yu, Park and Mammen (2008), Lee, Mammen and Park (2010, 2014, Zhang, Park and Wang (2013), among others.
The paper is structured as follows. Section 2 describes our structured insample density forecasting model, and show that the model is identifiable (estimable) under weak conditions. Section 3 explains a new approach to the estimation of the model. Here, it is assumed that the data are observed in continuous time and nonparametric smoothing methods are applied. Section 4 contains the theoretical properties of our method and Section 5 considers numerical examples and discusses the performance of the new approach. The Appendix contains technical details.
2. The model. We observe a random sample {(X i , Y i ) : 1 ≤ i ≤ n} from a density f supported on a subset I of a rectangle [0, 1] 2 . The density f (x, y) of (X i , Y i ) is a multiplicative function of three univariate components, where the first two are a function of the coordinate x and y, respectively, and the third is a function of the sum of the two coordinates, x + y, and is periodic. Specifically, we consider the following multiplicative model: where m J (t) = Jmod J (t), mod J (t) = t modulo 1/J for some J > 0, that is, m J (t) = J(t − l/J) for l/J ≤ t < (l + 1)/J , j = 0, 1, 2, . . . . Here, f j are unknown nonnegative functions supported and bounded away from zero on their supports. We note that m J (t) always takes values in [0, 1) as t varies on R + , and that the third component f 3 (m J (·)) is a periodic function with period J −1 .
We will prove the identifiability of the functions f 1 , f 2 and f 3 under the constraints that 1 0 f 1 (x) dx = 1 0 f 2 (y) dy = 1. We will do this for two scenarios. In the first case, we assume that f 1 , f 2 and f 3 are smooth functions. Then identification follows by a simple argument. Our second result does not make use of smoothness conditions of the component functions. It only requires conditions on the shape of the set I. The second result is important for an understanding of our estimation procedure that is based on a projection onto the model (2.1) without using a smoothing procedure for the component functions.
Our first identifiability result makes use of the following conditions: (A1) The projections of the set I onto the x-and y-axis equal [0, 1].

5
(A2) For every z ∈ [0, 1) there exists (x, y) in the interior of I with m J (x+ y) = z. Furthermore, for every x, y ∈ (0, 1) there exist x ′ and y ′ with (x, y ′ ) and (x ′ , y) in the interior of I.
(A3) The functions f 1 , f 2 , f 3 are bounded away from zero and infinity on their supports.
We now come to our second identifiability result that does not require smoothness conditions for the functions f 1 , f 2 and f 3 . This makes use of the following conditions on the shape of the support set I. To introduce conditions on the support set I, we let I 1 (y) = {x : (x, y) ∈ I}, I 2 (x) = {y : (x, y) ∈ I} and I 3l (z) = {x ∈ [0, 1] : (x, (z + l)/J − x) ∈ I}. Below, we assume that these sets change smoothly as y, x and z, respectively, move. Here, A△B denotes the symmetric difference of two sets A and B in R, and mes(A) the Lebesgue measure of a set A ⊂ R. Recall the definition T = max{x + y : (x, y) ∈ I}, and with this define L(J) be the largest integer that is less than or equal to T J .
The assumptions (A1), (A2), (A5) and (A6) accommodate a variety of sets I that arise in real applications. Figure 1 depicts some realistic examples of the set I that satisfy the assumptions. In particular, those sets of the type in the panels (c) and (e) satisfy (A2) and (A6) if the maximal vertical or horizontal thickness of the stripe is larger than the period 1/J of the third component function f 3 (m J (·)). In the interpretation of the examples in Figure 1, we follow the equivalent discussion from Keiding (1991) and Kuang, Nielsen and Nielsen (2008). The triangle in Figure 1(a) is typical for insurance or mortality when none of the underwriting years or cohorts are fully run-off. The standard actuarial term "fully run-off" means that all events from that underwriting year or cohort have been observed. In al-7 most all practical cases of estimating outstanding liabilities, actuaries stick to the triangle format leaving out fully run-off underwriting years. While the triangle also appears in mortality studies, it is common here to leave the fully run-off cohorts in the study resulting in the support shape given in Figure 1(b). The support in Figure 1(c) arises when the data analyst only considers observations from the most recent calendar years. While this approach is omnipresent in practical actuarial science, there is no formal theory or mathematical models behind these procedures in the actuarial literature. This paper is therefore an important step toward formalising mathematically actuarial practise while at the same time improving it. The support given in Figure 1(d) and (e) arises when there is a known time transformation such that time is running at another pace for different underwriting years or cohort years. While this type of time transformations are well known in mortality studies are often coined as versions of accelerated failure time models. Time transformations are also well known in actuarial science coined as operational time. However, the academic literature of actuarial science is still struggling to find a formal definition of what operational time is. This paper offers one potential solution to this outstanding and important issue. The last Figure 1(f) is included to give an impression of the generality of support structures one could deal with inside our model approach. Data is missing in the beginning and end of the delay period, but the model is still valid and in-sample forecasts can be constructed.
The model (2.1) has taken structured density forecasting into a new territory by leaving the simple multiplicative model. If f 3 above was constant (and therefore not in the model) then our model reduces to the simple multiplicative model analysed in Martínez-Miranda et al. (2013) and Mammen, Martínez-Miranda and Nielsen (2015). These two papers point out that the simple multiplicative density forecasting model is a continuous version of a widely used parametric approach corresponding to a structured histogram version of in-sample density forecasting based on the simple multiplicative model. The in-sample density forecasting model under investigation in this paper generalizes the simple multiplicative approach in an intuitive and simple way including seasonal effects.
In the following theorem, we show that, if there are two multiplicative representations of the joint density f that agree on almost all points in I, then the component functions also agree on almost all points in [0, 1]. We will use this result later in the asymptotic analysis of our estimation procedure.
Theorem 2. Assume that model (2.1) holds with (A1)-(A3), (A5), (A6). Suppose that (g 1 , g 2 , g 3 ) is a tuple of functions that are bounded away from zero and infinity with 1 0 g 1 (x) dx = 1 0 g 2 (y) dy = 1. Let µ j = log f j − log g j . Assume that µ 1 (x)+ µ 2 (y)+ µ 3 (m J (x+ y)) = 0 a.e. on I. Then µ j ≡ 0 a.e. on [0, 1]. 8 LEE, MAMMEN, NIELSEN AND PARK 3. Methodology. We describe the estimation method for the model (2.1). We first note that the marginal densities of X, Y and m J (X + Y ) may be zero even if we assume that the joint density is bounded away from zero. For example, the marginal densities of X and Y at the point u = 1 are zero for the support set I given in Figure 1(a). We estimate the multiplicative density model on a region where we observe sufficient data. This means that we exclude the points (1, 0) and (0, 1) in the estimation in the case of Figure 1(a), and the point (1, 0) in the case of Figure 1(b). Formally, for a set S ⊂ I, let J 1 and J 2 denote versions of I 1 and I 2 , respectively, defined by J 1 (y) = {x : (x, y) ∈ S} and J 2 (x) = {y : (x, y) ∈ S}, and define J 3l (z) = {x : (x, (z + l)/J − x) ∈ S}. We take an arbitrarily small number δ > 0, and find the largest set S such that where mes(A) for a set A denotes its length. Such a set is given by . We put the following constraints on f j : Then the model (2.1) gives the following integral equations: j are chosen so that the resultingf [k] j satisfy (3.2). We note that the naive two-dimensional kernel density estimator is not consistent near the boundary region, which jeopardizes the properties of the solution of the backfitting equation (3.3) at boundaries. For a preliminary estimatorf of the joint density f , we take the local linear estimation technique. The local linear estimatorf we consider here is similar in spirit to the proposal of Cheng (1997) where (h 1 , h 2 ) is the bandwidth vector and K is a symmetric univariate probability density function. Also, definê and 0 otherwise. The local linear density estimatorf we consider in this paper is defined byη 0 , whereη = (η 0 ,η 1 ,η 2 ) is given byη Before we close this section, we give two remarks. One is that, instead of integrating the two-dimensional estimatorf , one may estimate f w,j directly from the data. In particular, one may estimate f w,j by the one-dimensional kernel density estimators Our theory that we present in the next section is valid for this alternative estimation procedure. The other thing we would like to remark is that one may be also interested in an extension of the model (2.1) that arises when one observes a covariate U i ∈ R d along with (X i , Y i ). A natural extension of the model (2.1) in this case is that the conditional density of (X, Y ) given where the constraints (B1) now applies to f 1 (·, z) and f 2 (·, z) for each z. The method and theory for this extended model are easy to derive from those we present here.
We consider a theoretical approximation off . Define a nonlinear operator by G(θ, g) = F(1 + θ, f • (1 + g)), where g 1 • g 2 denotes the entry-wise multiplication of the two function vectors g 1 and g 2 . Then G(0, 0) Let G ′−1 : R 3 × S → R 3 × S denote the inverse of G ′ , whose existence we will prove in the Appendix. We definef = where g 1 /g 2 denotes the entrywise division of the function g 1 by g 2 .
It can be seen that δ = (δ 1 , δ 2 , given as the solution of the following system of integral equations: subject to the constraints In the following theorem, we show that the approximation off byf is good enough. In the theorem, we assume thatf (x, y) − f (x, y) = O p (ε n ) uniformly on S for some nonnegative sequence {ε n } that converges to zero as n tends to infinity. For the local linear estimatorf defined by (3.5) with h 1 ∼ h 2 ∼ n −1/5 , we have ε n = n −3/10 √ log n. The theorem tells that the approximation errors off j forf j are of order O p (n −3/5 log n). In Theorem 4 below, we will show thatf j − f j have magnitude of order O p (n −2/5 √ log n) uniformly on S j . This means that the first-order properties off j are the same as those off j .
Theorem 3. Assume that the conditions of Theorem 2 hold, and that the joint density f is bounded away from zero and infinity on its support S with continuous partial derivatives on the interior of S.
Next, we present the limit distribution of (f − f )/f . In the next theorem, we assume that h 1 ∼ c 1 n −1/5 and h 2 ∼ c 2 n −1/5 for some constants c 1 , c 2 > 0. 14 LEE, MAMMEN, NIELSEN AND PARK For such constants, definẽ Also, defineμ B j for j = 1, 2, 3 asμ j at (4.2) with the local linear estimator f being replaced byf B . In the Appendix, we will show that the asymptotic mean of (f j − f j )/f j equals n −2/5 β j , where β = (β 1 , β 2 , β 3 ) is the solution of the backfitting equation (4.4) withμ being replaced byμ B . Letf A denote the centered version of the naive two-dimensional kernel density estimator.
Here and below, we write K h (u) = K(u/h)/h. Defineμ A j for j = 1, 2, 3 as µ B j withf A taking the role off B . We will also show that the asymptotic variances of (f j − f j )/f j equal those ofμ A j , respectively, and that they are given by n −4/5 σ 2 j , where where K * K denotes the two-fold convolution of the kernel K.
In the discussion of assumption (A6) in Section 2, we note that (A6) allows a finite number of jumps in I j (u) for j = 1, 2 and I 3l (u) as u changes. These jump points are actually those where the marginal densities f w,j are discontinuous. At these discontinuity points, the expression of the asymptotic distributions of the estimators is complicate. For this reason, we consider only those points in the partitions (a j k−1 , a j k ), 1 ≤ k ≤ L j , for the asymptotic distribution off j , where a j k are the points that appear in assumption (A6). We denote by S j,c the resulting subset of S j after deleting all a j k , 1 ≤ k ≤ L j − 1. Note that f w,j is continuous on S j,c due to (A6). In the theorem below, we also denote by S o j the interiors of S j , j = 1, 2, 3.
For the limit distribution off j , we put an additional condition on the support set. To state the condition, let J o 2 (u 1 ; h 2 ) be a subset of J 2 (u 1 ) such that v ∈ J o 2 (u 1 ; h 2 ) if and only if v − h 2 t ∈ J 2 (u 1 ) for all t ∈ [−1, 1]. The set J o 2 (u 1 ; h 2 ) is inside J 2 (u 1 ) at a depth h 2 . In the following assumption, a j k and κ are the points and the function that appear in assumption (A6).
(A7) There exist constants C > 0 and α > 1/2 such that the following statements hold: (i) for any sequence of positive numbers ε n , Theorem 4. Assume that (A7) and the conditions of Theorem 3 hold, and that the joint density f is twice partially continuously differentiable. Let the kernel K be supported on [−1, 1], symmetric and Lipschitz continuous. Let the bandwidths h j satisfy n 1/5 h j → c j for some constants c j > 0. Then, for fixed points u j ∈ S o j ∩ S j,c , it holds that Remark 2. In the case where the third component function f 3 is constant, that is, there is no periodic component, the above theorem continue to hold for the component f 1 and f 2 without those conditions that pertain to the set S 3 and the function f 3 .

Numerical properties.
5.1. Simulation studies. We considered two densities on I = {(x, y) : 0 ≤ x, y ≤ 1, x + y ≤ 1}. Model 1 has the components f 1 ≡ f 2 ≡ 1 on [0, 1], and f 3 (u) = c 1 (sin(2πu) + 3/2), u ∈ [0, 1], where c 1 > 0 is chosen to make f (x, y) = f 1 (x)f 2 (y)f 3 (m J (x + y)) be a density on I. Model 2 has f 1 (u) = 3/2 − u, f 2 (u) = 5/4 − 3u 2 /4 and f 3 (u) = c 2 (u 3 − 3u 2 /2 + u/2 + 1/2) for some constant c 2 > 0. We took J = 2. We computed our estimates on a grid of bandwidth choice h 1 = h 2 . For model 1, we took {0.070 + 0.001 × j : 0 ≤ j ≤ 30} in the range [0.070, 0.100], and for model 2 we chose {0.40 + 0.02 × j : 0 ≤ j ≤ 20} in the range [0.40, 0.80]. In both cases, the ranges covered the optimal bandwidths. We obtained MISE based on 100 pseudo samples. The sample sizes were n = 400 and 1000, but only the results for n = 400 are reported since the lessons are the same. thus gives some indication of how sensitive our estimators are to the choice of bandwidth. The bandwidth that gave the minimal value of MISE 1 + MISE 2 + MISE 3 was h 1 = h 2 = 0.089 in model 1, and h 1 = h 2 = 0.64 in model 2, for the sample size n = 400. The values of MISE j along with ISB j and IV j for these optimal bandwidths are reported in Table 1. Although our primary concern is the estimation of the component functions, it is also of interest to see how good the produced two-dimensional density estimator f 1 (x)f 2 (y)f 3 (m J (x + y)) behaves. For this, we include in the table the values of MISE, ISB and IV of the two-dimensional estimates computed using the optimal bandwidth h 1 = h 2 = 0.089 in model 1, and h 1 = h 2 = 0.64 in model 2. For comparison, we also report the results for the two-dimensional local linear estimates defined at (3.5). For the local linear estimator, we used its optimal choices h 1 = h 2 = 0.085 in model 1, and h 1 = h 2 = 0.48 in model 2. We found that the initial local linear estimates had a large portion of mass outside I, and thus behaved very poorly if they were not re-scaled to be integrated to one on I. The reported values in Table 1 are for the adjusted local linear estimates. Overall, our two-dimensional estimator has better performance than the local linear estimator, especially in model 2. Figure 3 depicts the true density of model 1 and our two-dimensional estimate that has the median performance in terms of ISE.

Data examples.
The original data set we analyze in this section was collected between the year 1990 to 2011 by the major global UK based nonlife insurance company RSA. The dataset-and more details about it-is publicly available via the Cass Business School web site together with the paper "Double Chain Ladder" at the Cass knowledge site. The observations were the incurred counts of large claims aggregated by months. During the 264 months, 1516 large claims were made. The dataset is provided in the form of a classical run-off triangle {N kl : 1 ≤ k, l ≤ 264, k + l ≤ 265}, where N kl denotes the number of large claims incurred in the kth month and reported in the (k + l − 1)th month, that is, with (l − 1) months delay. Since the data are grouped monthly, we need pre-smoothing of the data to apply the model (2.1) that is based on data recorded over a continuous time scale. A natural way of pre-smoothing is to perturb the data by uniform random variables. Thus, we converted each claim (k, l) on the two-dimensional discrete time scale {(k, l) : 1 ≤ k, l ≤ 264, k + l ≤ 265}, into (X, Y ) on the two-dimensional continuous time scale I = {(x, y) : 0 ≤ x, y ≤ 1, x + y ≤ 1}, by where (U 1 , U 2 ) is a two-dimensional uniform random variate on the unit square [0, 1] 2 . This gives a converted dataset {(X i , Y i ) : 1 ≤ i ≤ 1516}. We applied to this dataset our method of estimating the structured density f of (X, Y ). Since one month corresponds to an interval with length 1/264 on the [0, 1] scale, one year is equivalent to an interval with length 12/264 = 1/22 on the latter scale. We let the periodic component f 3 (m J (·)) in the model (2.1) reflect a possible seasonal effect, so that we take one year in the real time to be the period of the function. This means that we let the periodic component f 3 (m J (·)) have 1/22 as its period, and thus take J = 22. For the bandwidth, we took h 1 = h 2 = 0.01. The chosen bandwidth may be considered to be too small for the estimation of f 1 and f 2 . However, we took such a small bandwidth to detect possible seasonality. Note that the bandwidth size 0.01 corresponds to 0.01 × 12 × 22 = 2.64 months. We found that even with this small bandwidth the estimated curvef 3 was nearly a constant function, which suggests that the large claim data do not have a seasonal effect.
To see how well our method detects a possible seasonal effect in the data, we augmented the dataset by adding a certain level of seasonal effect as follows. We computed Since (k + l − 1 modulo 12) is the actual month of the claims reported, the augmented dataset has added claims in November, December, January and February. The augmentation resulted in increasing the total number of claims to 2606 from 1516. The increased counts of reported claims were 252 from 126 for November, 600 from 200 for December, 455 from 91 for January and 300 from 100 for February. In our estimation procedure, the bandwidths h 1 and h 2 control the smoothness of the local linear estimatef along the x-and y-axis, respectively. Consequently, choosing small values for h 1 and h 2 would result in nonsmooth estimates of the functions f 1 and f 2 , which we observed in the pilot study with h 1 = h 2 = 0.01. Nevertheless, in some cases setting these bandwidths to be small, relative to the scales of X and Y , might be preferred when one needs to detect possible seasonality, as is the case with the current dataset. In our dataset, the bandwidth size 1/264 = 0.0038 on the scale of [0, 1] corresponds to one month in real time. Thus, taking the bandwidths to be 0.015, for example, that corresponds to a period of four months, forces the seasonal effect to almost vanish in the estimate of f 3 .
To achieve both aims of producing smooth estimates of f 1 and f 2 , and of detecting possible seasonal effect, we applied to the augmented dataset a two-stage procedure that is based on our estimation method described in Section 3. In the first stage, we got a local linear estimatef with h 1 = h 2 = 0.01, and found an estimate of f 3 using the iteration scheme at (3.4). In the second stage, we recomputed a local linear estimatef with larger bandwidths h 1 = h 2 = 0.05, and found estimates of f 1 and f 2 using only the first two updating equations at (3.4) withf being replaced by the estimate of f 3 obtained in the first stage.
The results of applying this two-stage procedure to the augmented dataset are presented in Figure 4. Clearly, the seasonal effect of the augmented dataset was well recovered in the estimate of f 3 , and at the same time smooth estimates of f 1 and f 2 were produced. The augmented data set indicate an increased number of claims in the winter time. This is clearly reflected in the estimated results, where the first part and the last part of the estimated effect is higher than the rest of the curve. Imagine the realistic situation that a nonlife insurer on the first day of November has to produce budget expenses for the rest of the year. The classical multiplicative methodology is not able to reflect the two month perspective of such a budget. Therefore, considerable work is being done manually in finance and actuarial departments of nonlife insurance companies to correct for such effects. With our new seasonal correction, costly manual procedures can be replaced by cost saving automatic ones eventually benefitting the prices all of us as end customers have to pay for insurance products. Figure 5 depicts the resulting two-dimensional joint density. Notice that this two-dimensional density is clearly nonmultiplicative. The seasonal correction provides a visually deviation from the multiplicative shape. Also, note that while this two-dimensional density is nonmultiplicative, the nature of this deviation is not immediately clear to the eye. Whether the deviation is pure noise, a seasonal effect or some other effect is not easy to get from the full two-dimensional graph of the local linear density estimate which is also presented in Figure 5. For the local linear estimate, we used h 1 = h 2 = 0.03. We tried other bandwidth choices such as 0.01 and 0.05, but found that the smaller one gave too rough estimate and the larger one produced too smooth a surface. Our two-dimensional density estimate therefore illustrates why research into structured densities on nontrivial supports is crucial to extract information beyond the classical and simple multiplicative one. APPENDIX A.1. Proof of Theorem 1. Suppose that (g 1 , g 2 , g 3 ) is a tuple of functions that are bounded away from zero and infinity with 1 0 g 1 (x) dx = 1 0 g 2 (y) dy = 1 and f (x, y) = g 1 (x)g 2 (y)g 3 (m J (x + y)).
We now argue that µ 3 is continuous on [0, 1). To see that µ 3 is continuous at z 0 ∈ [0, 1), we choose (x 0 , y 0 ) in the interior of I such that m J (x 0 + y 0 ) = z 0 . This is possible because of assumption (A2). We can choose x 0 and y 0 such that µ 1 is continuous at x 0 and µ 2 is continuous at y 0 . Thus, we get from (A.1) that µ 3 is continuous at z 0 . Similarly, one shows that µ 1 and µ 2 are continuous functions on [0,1]. This gives that for all x, y ∈ (0, 1).
Theorem 3 basically follows from the above (i)-(iii). To prove the theorem using (i)-(iii), we note that claim (ii) with the definitions ofθ andf at (4.3) Now, from (ii) it follows that there exists a constant C > 0 such that the mapĜ ′ (θ − 1, with probability tending to one. Also, (iii) is valid for all (θ 1 , g 1 ), (θ 2 , g 2 ) ∈ B 2r (θ − 1, (f − f )/f ). Then we can argue that the solution of the equationĜ(θ, g) = 0, which is (θ − 1, (f − f )/f ), is within Cα n distance from (θ − 1, (f − f )/f ), with probability tending to one, where C > 0 is a constant and α n = Ĝ (θ − 1, (f − f )/f ) ∞ . This follows from an application of the Newton-Kantorovich theorem; see Deimling (1985) or Yu, Park and Mammen (2008) for a statement of the theorem and related applications. To compute α n , we note that . For the first equation of (A.4), we have used (iii) and the facts thatθ − 1 = O p (ε n ) and (f − f )/f = O p (ε n ). The second equation of (A.4) follows from the inequality . Claim (i) follows from the uniform convergence off to f that is assumed in the theorem: sup (x,y)∈S |f (x, y) − f (x, y)| = O p (ε n ). Below, we give the proofs of claims (ii) and (iii).
Going through the arguments in the proof of G ′ (0, 0) being one-to-one and now using the first two equations of (A.6) give η 1 = η 2 = η 3 ≡ 0. Note that the first two equations can be written as S 1 f w,1 (x)η 1 (x) dx = 0 and S 2 f w,2 (y)η 2 (y) dy = 0, and thus in the latter proof f w,j for j = 1, 2 take the roles of f j in the former proof. The foregoing arguments show that (0, 0) is the only tuple that is perpendicular to the range space of G ′ (0, 0), which implies that G ′ (0, 0) is onto.
To verify that the inverse map G ′ (0, 0) −1 is bounded, it suffices to prove that the bijective linear operator G ′ (0, 0) is bounded, owing to the bounded inverse theorem. Indeed, it holds that there exists a constant C > 0 such that G ′ (0, 0)(d, δ) ∞ ≤ C (d, δ) ∞ . This completes the proof of claim (ii).
This completes the proof of (iii).
Now, for the deterministic partf B , recall the definitions off B andμ B at (4.6) and thereafter, respectively. Let r n =μ B − n −2/5μB . According to Lemma 1, r n = o(n −2/5 ) on S ′ 1 × S ′ 2 × S ′ 3 , where S ′ j is a subset of S j with the property that mes(S j − S ′ j ) = O(n −1/5 ). We also get r n = O(n −2/5 ) on S 1 × S 2 × S 3 . This implies T(r n ) = o(n −2/5 ), so that uniformly on S 1 × S 2 × S 3 . Thus, (f B /f ,θ B − 1) equals the solution of the backfitting equation δ = d + n −2/5μB − Tδ, up to an additive term whose jth component has a magnitude of an order o(n −2/5 ) on S ′ j and O(n −2/5 ) on the whole set S j .
The asymptotic distribution of ((f j (u j ) − f j (u j ))/f j (u j ) : 1 ≤ j ≤ 3) for fixed u j ∈ S j,c ∩ S o j is then readily obtained from the above results. The asymptotic mean is given as the solution (δ j (u j ) : 1 ≤ j ≤ 3) of the backfitting equation (4.4) withμ j being replaced by n −2/5μB j , subject to the constraint (4.5). The asymptotic variances are derived from those ofμ A j , wherẽ ]. This is due to (A.9), (A.10) and the corresponding property forμ A 3 in the proof of Lemma 2 below.