A Spatial-Nonparametric Approach for Prediction of Claim Frequency in Motor Insurance

Gideon Kipngetich¹, Ananda Kube², Thomas Mageto³

¹Pan African University, Nairobi, Kenya.
²Department of Mathematics, Statistics and Actuarial Sciences, Kenyatta University, Nairobi, Kenya.
³Department of Statistics and Actuarial Sciences, JKUAT, Nairobi, Kenya.

DOI: 10.4236/ojs.2021.114031                

Abstract

Spatial modeling has largely been applied in epidemiology and disease modeling. Different methods such as Generalized linear models (GLMs) have been made available to prediction of the claim frequencies. However, due to heterogeneity nature of policies, the methods do not generate precise and accurate claim frequencies predictions; these parametric statistical methods extensively depend on limiting assumptions (linearity, normality, independence among predictor variables, and a pre-existing functional form relating the criterion variable and predictive variables). This study investigates how to derive a spatial nonparametric model estimator based on smoothing Spline for predicting claim frequencies. The simulation results showed that the proposed estimator is efficient for prediction of claim frequencies than the kernel based counterpart. The estimator derived was applied to a sample of 6500 observations obtained from Cooperative Insurance Company, Kenya for the period of 2018-2020 and the results showed that the proposed method performs better than the kernel based counterpart. It is worth noting that inclusion of the spatial effects significantly improves the estimator prediction of claim frequency.

Keywords

Nonparametric, SAR, Smoothing Spline, Claims, CIC, Spatial

1. Introduction

Recent studies on spatial modeling have been rapidly applied in many fields: epidemiology, public health, and the insurance sector. Models such as Poisson, Generalized linear models, Credibility models and Bayesian Models are the commonly used models for prediction of claim frequencies. However, from the available literature, these models appear to be relatively inflexible. Although the Generalized linear models provide accurate and fast analysis of insurance data, they fall short because they are defined based on the assumptions, and an incorrect model assumption can cause model misspecification leading to erroneous results. Nonparametric models are deemed to minimize the shortcoming of these standard parametric models since fewer assumptions are made for the model, therefore, suitable for modeling insurance data which are nonlinear as described by [1] where they concluded that the nonparametric models perform better than generalized linear models (GLMs); the only observable problem with modeling using Nonparametric models is the interpretation of some of the curves [2]. When modeling claims and risks, we need to determine their behavior and spatial dependence, and spatial heterogeneity of the data so that the insurer can determine which areas are associated with a higher riskier when determining premiums amount to be paid.

[3] proposed a Bayesian nonparametric approach for prediction of claims; here they found out that the model performs better compared to nonparametric GLMs in that it can capture the nonlinear random effects present in the data. [4] also proposed a flexible nonparametric loss model for prediction of the claims; they found out that having flexible multivariate model may allow actuaries to estimate the dependence between different risk classes and different lines of business and this topic needs to be explored further. [5] introduced the idea of using nonparametric data mining approach to modeling the claims and prediction of risk; here the approach classified the risk and predicted claim size based on the data. This study’s research idea was to be built based on the idea proposed by [6] [7] where they introduce a nonparametric spatial regression model for prediction. This study’s primary objective was to derive a spatial nonparametric estimator for the prediction of insurance claims. Therefore, this study’s main contribution was to investigate the estimator’s performance in the situation of additional covariates in the model and incorporate the aspect of spatial dependence in deriving a nonparametric estimator for the prediction of insurance claim frequencies.

The main difference between this research and [7] [8] [9] are as follows: 1) estimate a nonparametric spatial model where estimation of the unknown trend $g(\cdot )$ is based on smoothing spline; 2) spatial heterogeneity and correlation were considered simultaneously rather than assuming that correlation satisfies specific form (such as in SAR).

The paper is organized as follows: Section 2 describes the development process of the model estimator based on smoothing spline; Section 3 presents the data description and main results; simulation study and analysis of CIC insurance claims data; Section 4 presents conclusion and suggestions for further research.

2. Methods

Model Estimation

The study proposed a nonparametric regression model to predict the number of claims $Y_i,i=1,\cdots ,n$ observed in region J in order to relax restrictive assumption on the distribution of number of claims and $X_i$ covariates vector for the i^th claim. Since claims in each region, J has nonlinear relation with the covariates ${X^{\prime }}_{i}s$.

The nonparametric form of the model is given by the general form [10] [11].

$y_i=g\left(x_i\right)+Z_i^{\text{T}}b+{\epsilon }_i$

$g(\cdot )$ is unknown nonparametric function used to model fixed effects, $Z_i^{\text{T}}b$ and ${\epsilon }_i$ cater for random effects.

Since the form of $Z_i^{\text{T}}b=R_i$ for $R_i$ is unknown. The main work of this study is to estimate the form of $R_i$ then establish the functional form of $Z_i^{\text{T}}b$ that captures the spatial effects.

Let $n={\left(n_1\mathrm{,}\cdots \mathrm{,}{n}_N\right)}^{\text{T}}$ be two N dimensional vectors. We can make assumption about the spatial model as

$Y_i=g\left(X_i\right)+R_i,var\left(R\right)=\Sigma ,i\in {\Lambda }_n=\left\{1,\cdots ,n_1\right\}\times \cdots \times \left\{1,\cdots ,n_N\right\}$ (1)

where $i=\left(i_1\mathrm{,}\cdots \mathrm{,}{i}_N\right)$ in ${\Lambda }_n$ will be referred to as site, $R_i$ cater for the spatial effects (Random effects) and the cardinality of ${\Lambda }_n$ is $\left|{\Lambda }_n\right|=n$ [12].

Spatial data is modelled as finite realization of vector stochastic process indexed by $i\in {\Lambda }_n$, $R={\left(R_1\mathrm{,}\cdots \mathrm{,}{R}_n\right)}^{\text{T}}$ is assumed to follow a joint Gaussian distribution where $E\left(R_i\right)=0$, is known $\forall i\in {\Lambda }_n$, $\Sigma =\left[\rho \left(R_i\mathrm{,}{R}_j\right)\right]$ is the unknown correlation coefficient matrix (need to be estimated). The vector $X_i=\left(X_{i1}\mathrm{,}\cdots \mathrm{,}{X}_{id}\right)\in {\Re }^d$, $Y_i\in \Re$ and $g(\cdot )$ is the unknown trend function.

The aim is to estimate $g\left(x\right)$ for some given $x=\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_d\right)\in {\Re }^d$, the response variable $Y_i$ is claim frequency and $X_i$ is six dimensional vector consisting of the following explanatory variables: gender, claim amount, age of the policyholder, gender, vehicle age, model of the vehicle and age category of the policyholder.

Estimating $g\left(x\right)$ at some point $x\in {\Re }^d$, for $X_i$ in the neighbourhood of x, g can be approximated using smoothing spline [13] [14].

To estimate the smoothing spline estimator $\hat{g}(\cdot )$ of $g(\cdot )$, the study considers minimizing the equation

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left(Y_i-g\left(x_i\right)\right)}^2+\lambda {\displaystyle \int }{\left(g^{″}\left(x\right)\right)}^2\text{d}x$ (2)

over the function g This criterion trades-off least squares error of g over $\left(x_i,y_i\right),i=1,\cdots ,n$, with a regularization term that grows large when the second derivative of g is wiggly. The coefficients are chosen to minimize Equation (3) which is a simplified form of Equation (2)

$\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left\{Y_i-g\left(X_i;\beta \right)\right\}}^2+\lambda {\beta }^{\text{T}}\Omega \beta$ (3)

which can be represented as

${\Vert Y_i-G\beta \Vert }^2+\lambda {\beta }^{\text{T}}\Omega \beta$

where $G\in {ℝ}^{n\times n}$ is basis matrix defined as

$G_{ij}={\psi }_j\left(x_i\right),i,j=1,\cdots ,n$

where ${\psi }_1\mathrm{,}\cdots ,{\psi }_n$ are the truncated power basis functions with knots at $x_1\mathrm{,}\cdots ,x_n$ which is evaluated at the data values

${\psi }_j\left(x\right)=(\begin{array}{l}{x}_i^j\left(0\le n\le p\right)i,j=1,\cdots ,n\\ {\left(x_i-N_{j+1-p}\right)}_{+}^p\left(p+1\right)\le j\le N\end{array}$ (4)

${\left(x-N_{j+1-p}\right)}_{+}^p=\max {\left(0,x_i-N_j\right)}^p,j\in \varphi$ where $\varphi$ is compact interval. p is the degree of the spline and $j_i<\cdots <j_{N-p}$ are fixed points or knots in $\varphi$.

$\Omega \in {ℝ}^{n\times n}$ is the penalty matrix defined as

${\Omega }_{ij}={\displaystyle \int {g^{″}}_i\left(x\right){{\psi }^{″}}_j\left(x\right)\text{d}x}\mathrm{,}i\mathrm{,}j=\mathrm{1,}\cdots ,n$

Given the optimal coefficients $\hat{\beta }$ minimizing (3) through penalized least squares, the smoothing spline estimator at x is therefore defined as

$\hat{g}\left(x_i\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\hat{\beta }}_j{\psi }_j\left(x\right)$ (5)

The term affects shrinking the components of estimation $\hat{\beta }$ towards zero. The parameter $\lambda \ge 0$ is the smoothing parameter.

Each computed coefficient ${\hat{\beta }}_j$ corresponds to a particular basis function ${\psi }_j$. The term ${\beta }^{\text{T}}\Omega \beta$ in (3) imparts more shrinkage on the coefficients ${\hat{\beta }}_j$ that correspond to wigglier functions ${\psi }_j\left(x\right)$. Hence, as we increase $\lambda$, we are shrinking away from the wiggler basis functions.

Similar to least squares regression, the coefficients $\hat{\beta }$ minimizing (3) is

$\hat{\beta }={\left(G^{\text{T}}G+\lambda \Omega \right)}^{-1}{G}^{\text{T}}Y={\left(X^{\text{T}}X+n\lambda D\right)}^{-1}{X}^{\text{T}}Y$

where X is a design matrix with entries $x_i$ for $i=1,\cdots ,n$, Y is a vector of the response variables, D is a diagonal matrix with $p+1$ zeros on the diagonal followed by N ones and $n\lambda D$ is a penalty term.

Smoothing splines can be seen as a linear smoother, where $k\left(x\right)=\left({\psi }_1\left(x_1\right),\cdots ,{\psi }_n\left(x_n\right)\right)$. Therefore, Equation (5) can be represented as

$\hat{g}\left(x\right)=k{\left(x\right)}^{\text{T}}\hat{\beta }=k{\left(x\right)}^{\text{T}}{\left(X^{\text{T}}X+n\lambda D\right)}^{-1}{X}^{\text{T}}Y$ (6)

which is linear combination of the points $y_i,i=1,\cdots ,n$, $\lambda$ is estimated using Generalized Cross Validation (GCV) method given by

$\text{GCV}\left(\lambda \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left(\frac{Y\left(z_i\right)-{\hat{Y}}_{\lambda }^{-i}\left(z_i\right)}{1-\left(p+tr\left(S_{\lambda }\right)\right)/n}\right)}^2$ (7)

where $Y\left(z_i\right)$ is the observation in point $z_i$, $Y_{\lambda }^{-i}\left(z_i\right)$ is the predicted value from a fitted smoothing spline model from the data less the i^th data and $S_{\lambda }$ is the degree of the smoother.

As proposed by [6] [7], R² is used to assess the performance of predictor function, given by

$R^2=1-\frac{{\displaystyle {\sum }_{i=1}^n{\left[g\left(x_i\right)-\hat{g}\left(x_i\right)\right]}^2}}{{\displaystyle {\sum }_{i=1}^n{\left[g\left(x_i\right)-\bar{g}\right]}^2}}$ (8)

where $\bar{g}$ is the sample mean of $g\left(x_i\right),i=1,\cdots ,n$.

After estimating the function $g(\cdot )$, then from (1) $R_i$ is estimated as ${\hat{R}}_i=Y_i-\hat{g}\left(X_i\right)$. Since $\Sigma$ in model Equation (1) is unknown, we assume that $R_i,i=1,2,\cdots ,n$ is 2^nd-order stationary and isotopic process (does not depend on direction).

Before prediction can be performed on spatial data sets, the variogram is usually estimated at various lags and a nonparametric model is fitted to those estimates.

Then let $C\left(h\right)$ and $2\gamma \left(h\right)$ be covariogram and variogram of the process where h represents the distance between 2 points at which the process is obtained [12] [15]. The two quantities are related by

$C\left(h\right)=C\left(0\right)-\gamma \left(h\right)$ (9)

where $C\left(0\right)={\sigma }^2=var\left(Y\left(z\right)\right)$, $Y\left(z\right)$ is the value of the process at spatial location z within region C.

$\underset{h\to \infty }{lim}C\left(h\right)=0$

implies

$\underset{h\to \infty }{lim}\gamma \left(h\right)=Var\left(Y\left(z\right)\right)=C(\; 0\; )$

for validity of variogram the condition that

$\underset{h\to \infty }{lim}\frac{2\gamma \left(h\right)}{h^2}=0$

must be met [16].

$\Sigma =\left[\rho \left(R_i\mathrm{,}{R}_j\right)\right]=\left[C\left(\Vert z_i-z_j\Vert \right)/{\sigma }^2\right]$, while $z_i$ and $z_j$ are the spatial locations associated with the error values $R_i$ and $R_j$ thus to estimate $\Sigma$ it is sufficient to estimate $\gamma \left(h\right)$ [16] [17].

$2\hat{\gamma }\left(h\right)={\displaystyle \underset{S\left(h\right)}{\sum }}{\left[z_i-z_j\right]}^2/N\left(h\right)$ (10)

$S\left(h\right)=\left\{\left(z_i,z_j\right):\left|z_i-z_j\right|=h\right\}$, $h\in {\Re }^d$, $N\left(h\right)$ is a number of distinct pairs in $S\left(h\right)$ since $r\left(z_i\right)$ the error at location $z_i$ is unobserved, the quantity is to be estimated as well.

Since we have to estimate the variogram $\hat{\gamma }\left(h\right)$ in Equation (10) in nonparametric approach [18] [19], then $\gamma \left(h\right)$ can be estimated as

$\gamma \left(h\right)={\displaystyle {\int }_0^{\infty }}\left(1-{\omega }_d\left(ht\right)\right)\text{d}M\left(t\right)$ (11)

$M\left(t\right)$ is nonnegative bounded nondecreasing function for nodes(or location of the jumps) $t\ge 0$ and ${\omega }_d$ is a basis for functions in ${ℝ}^d$ (d is the dimension of the spatial domain D) given by

${\omega }_d\left(ht\right)={\left(2/ht\right)}^{\left(d-2\right)/2}\Gamma \left(d/2\right)J_{\left(d-2\right)/2}(\; h\; t\; )$

$\Gamma \left(d/2\right)$ is the gamma function, and $J_{\mathrm{<mo\; stretchy="false">\; (\; </mo>}\cdot \mathrm{<mo\; stretchy="false">\; )\; </mo>}}$ is the Bessel function of the first kind. Some familiar examples of ${\omega }_d$ are ${\omega }_1\left(ht\right)=\cos \left(ht\right)$, ${\omega }_2\left(ht\right)=J_0\left(ht\right)$, and here ${\omega }_3\left(ht\right)=\frac{\sin \left(ht\right)}{ht}$ is chosen which yields a non-parametric estimate which is conditionally negative definite for spatial data from 1 - 3 dimensions.

The characteristics of the estimator (11) are estimated using Integrated square error [20], given by

$\text{ISE}\left(\gamma \right)={\displaystyle {\int }_{h_1}^{h_k}}{\left\{\hat{\gamma }\left(h\right)-\gamma \left(h\right)\right\}}^2\text{d}h$ (12)

where $h_1$ and $h_k$ are the smallest and largest distances for which variogram estimates are available [17].

Model (1) can therefore be represented as

$Y\left(z_i\right)=g\left(X_i\left(z_i\right)\right)+R\left(z_i\right)\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}n$ (13)

where $Y\left(z_i\right)\mathrm{<mo>\; :\; </mo>}i=\mathrm{1,}\cdots \mathrm{,}n$ is the observations (claims) in region $z_i$ associated with independent variables $X_i\left(z_i\right)$ in region $z_i$, $R\left(z_i\right)$ is the unobserved error in region $z_i$ and $g(\cdot )$ is the estimated function in (6).

To evaluate performance of the proposed method we used $R^2$ to assess prediction accuracy of the method

$R^2=1-\frac{{\displaystyle {\sum }_{i=1}^n{\left[Y\left(z_i\right)-\hat{Y}\left(z_i\right)\right]}^2}}{{\displaystyle {\sum }_{i=1}^n{\left[Y\left(z_i\right)-\bar{Y}\right]}^2}}$ (14)

3. Main Results

3.1. Data Description

The study used motor third party liability data for 2018-2020 from the insurance company Cooperative Insurance Company (CIC). The data include 6500 policies, out of which many policies have total claim sizes other than zero, and an appropriate number of policies without any claims were taken. The following policy data were used: the region where the policy was taken, age, gender, type of vehicle, number of claims per policy, years of policy ownership, claim amount, insured cases number for a user, and average claim size. In the process of preparation, data was cleaned, and imputation of data will be done; age is categorized into old (over the age of 50), Young (up to the age of 25), and Middle (aged 25 - 50) age. Policies with extremely low and extremely high average claim sizes are removed; categorical variables with multiple categories were replaced with dummy (indicator) variables.

Table 1 shows that there is a very large number of observations with no claims in claims dataset where the maximum number of claims made in a region was 4 in an observation.

Table 1. Summary of claim frequencies in the data.

3.2. Simulation Study

This section describes the simulation and their analysis results of the proposed method, we simulate spatial data with a length of n = 100 observations. This is to ensure that the simulated data mimic the real claims dataset so that the results can be inferred to evaluate the performance of our method in data analysis. 65 spatial sampling locations were selected randomly and denoted by $z_1\mathrm{,}\cdots \mathrm{,}{z}_n$. The responses $Y\left(z_i\right)$ for $i=\mathrm{1,}\cdots \mathrm{,}n$ are the observations and were simulated from the spatial nonparametric model (13) with $p=2$

$Y\left(z\right)=g\left(x_i\left(z\right)\right)+R(\; z\; )$

$R\left(z\right)$ is the term for spatial effects $z_i$ and $z_j$ in 2-dimensional space [16] with mean 0 and covariance given by (11). The covariates, $x_i\left(z_i\right)$ for $i=\mathrm{1,}\cdots \mathrm{,}n$, were generated as iid $N\left(\mathrm{0,1}\right)$ and are independent of each other, before the simulations i.e., the variables were treated as fixed terms when $Y\left(z_i\right)$ were generated repeatedly. Within each simulation, the spatial random effects $R\left(z_i\right)$ were generated from a Gaussian process with mean zero and the covariance function (11), for $i=\mathrm{1,}\cdots \mathrm{,}n$.

The $\text{ISE}\left(\gamma \right)$ defined as $\text{ISE}\left(\gamma \right)={\displaystyle {\int }_{h_1}^{h_k}}{\left\{\hat{\gamma }\left(h\right)-\gamma \left(h\right)\right\}}^2\text{d}h$ was approximated numerically from simulated data for the proposed estimator (11) and NW kernel. Table 2 present the mean values of ISE. From the results in Table, the proposed estimator (11) offers a better performance compared to NW kernel estimator.

Assessing how well the proposed method performs, we compare the proposed method under which $R={\displaystyle {\int }_0^{\infty }}\left(1-{\omega }_d\left(ht\right)\right)\text{d}M\left(t\right)$ with the method under which the spatial component (R) is based on kernel estimation, we calculated the MSE and the R² of the estimators from 100 simulations and present the results in Table 3. From the table it was found that MSE for the proposed method is smaller ranging from 0.0221 to 0.0102 in all the sample sizes taken while the MSE of the kernel based estimator ranges between 0.308 to 0.0176, in addition, the R² for the proposed method were larger ranging between 0.7003 to 0.99963 in all the sample sizes compared those of kernel based estimator which ranges from 0.6751 to 0.9694. Thus the results demonstrate the superior performance of the proposed method compared to the kernel based estimator.

The results in Table 3 were visualized in Figure 1 and Figure 2.

Based on performance of the proposed method, the method was applied to the simulated data to check its performance in prediction of future values. Table 4 describes the distribution of the predicted values out of 100 simulation.

Table 2. Mean values of the standardized ISE from the estimators (Sample size of n = 1000).

Table 3. MSE and R² for the model over 100 simulation under different sample sizes.

Figure 1. R-squared plot.

Figure 2. Mean squared plot of both K (kernel) and proposed estimator.

Figure 3. Histogram for the predicted values.

Table 4. Summary of predicted claim frequencies from simulation.

The predicted values generated by the proposed method as presented in Table 4 were graphically presented and the prediction intervals superimposed on the distributional histogram of the predicted values as shown in Figure 3, the prediction interval (in red dotted lines) showed that a larger number of predicted values lies between 1 - 2 this means that there are higher chances of getting future values as 1 and 2.

3.3. Analysis for Claims Data

The study considered claims data from CIC insurance observed in different parts of 7 counties of Kenya to exhibit the performance of the proposed method. The main interest of this study was predicting claims frequencies, the study considers a set of 6500 observations. Let $Y_i$ denote the claim frequency, and $X_i={\left(X_1\mathrm{,}\cdots \mathrm{,}{X}_6\right)}^{{\rm T}}$ be a vector which consists of the following explanatory variables: gender, claim amount, age of the policyholder, gender, vehicle age, model of the vehicle and age category of the policyholder. Using the estimated model (13) we predict claim frequencies. The observations were from random process over a countable sample of spatial locations. The claim data at a particular location typically represent the entire region (Figure 4).

Using the proposed method future claim frequencies were predicted and the results were presented in Table 5.

Figure 5 shows graphical representation of the predicted claims as described in Table 5 with the prediction interval (red dotted vertical lines), the future values will lie between 1 and 4.

From the prediction results, R² values using Equation (14) were evaluated to access the performance of two methods, the results presented in Table 6.

From the results in Table 6 the R² for N(kernel) was 0.543, and that from the proposed method is 0.566, this showed that the proposed method for prediction has a higher prediction accuracy than the kernel based estimator. Therefore the

Table 5. Summary of predicted claim frequencies.

Table 6. R² for the estimators.

Figure 4. Hotspots locations in Nakuru, Nairobi, Kajiado, Muranga, Kiambu, Machakos, Makueni counties.

Figure 5. Predicted number of claims.

study concluded that the proposed method is more efficient than N (kernel) model, this implies that the predicted value was more likely to be more identically equal to the observed claims.

4. Conclusions and Suggestions

The idea of deriving an appropriate estimator in predicting frequency claims in the insurance industry has gained more interest in finance and statistical research. Many researchers heavily rely on parametric estimators; however, the insurance datasets have some aspect of non-linearity. Hence, researchers in statistics and econometrics are currently developing nonparametric models incorporating spatial effects to improve on the prediction based on the existing parametric models such as aggregate claim models and GLMs which are rather more restrictive on their transformed mean of the response; the nonparametric methods provide a more flexible method for prediction. The study proposed a spatial nonparametric (based on splines) estimator for predicting claim frequencies in motor insurance.

The simulation study showed that the proposed method performs better than the kernel based estimator; here the Mean Squared Error values of the proposed method were smaller than those of the kernel estimator which also implies a higher value of R-squared, particularly in presence of spatial dependence. Case study findings also showed that the proposed method performs better than the kernel based estimator on predicting the future claim frequencies. Therefore, the proposed method compared to kernel based estimator provides a more efficient prediction method for motor insurance claim data and ultimately leads to more accurate predictions.

Suggestions

Some additional exogenous variables such as environmental among other institutional factors may have effect on claim frequencies therefore, more robust spatial estimator need to be constructed using the proposed idea to investigate how these factors may affect claim frequencies. Further research can also be done on the theoretical properties of this proposed model estimator. In addition, this study made the assumption that the errors were correlated, for this reason future studies could consider a case of uncorrelated error structure.

Acknowledgements

Sincere thanks to my supervisors Dr. Kube Anada and Dr. Thomas Mageto for their professional contribution and performance, and special thanks to my parents for their moral support and rare attitude of high quality.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Doupe, P., Faghmous, J. and Basu, S. (2019) Machine Learning for Health Services Researchers. Value in Health, 22, 808-815. https://doi.org/10.1016/j.jval.2019.02.012
[2]	Islamiyati, A., Chamidah, N., et al. (2020) Penalized Spline Estimator with Multi Smoothing Parameters in Bi-Response Multi-Predictor Nonparametric Regression Model for Longitudinal Data. Songklanakarin Journal of Science & Technology, 42, 897-909.
[3]	Fellingham, G.W., Kottas, A. and Hartman, B.M. (2015) Bayesian Nonparametric Predictive Modeling of Group Health Claims. Insurance: Mathematics and Economics, 60, 1-10. https://doi.org/10.1016/j.insmatheco.2014.10.011
[4]	Hong, L. and Martin, R. (2017) A Flexible Bayesian Nonparametric Model for Predicting Future Insurance Claims. North American Actuarial Journal, 21, 228-241. https://doi.org/10.1080/10920277.2016.1247720
[5]	Kascelan, V., Kascelan, L. and Buric, M.N. (2016) A Nonparametric Data Mining Approach for Risk Prediction in Car Insurance: A Case Study from the Montenegrin Market. Economic Research—Ekonomska istrazivanja, 29, 545-558. https://doi.org/10.1080/1331677X.2016.1175729
[6]	Wang, H.X., Wang, J.D. and Huang, B. (2012) Prediction for Spatiotemporal Models with Autoregression in Errors. Journal of Non-Parametric Statistics, 24, 217-244. https://doi.org/10.1080/10485252.2011.616893
[7]	Wang, H.X., Lin, J.G. and Wang, J.D. (2016) Nonparametric Spatial Regression with Spatial Autoregressive Error Structure. Statistics, 50, 60-75. https://doi.org/10.1080/02331888.2015.1094068
[8]	Li, Y., Qin, Y. and Li, Y. (2021) Empirical Likelihood for Nonparametric Regression models with Spatial Autoregressive Errors. Journal of the Korean Statistical Society, 50, 447-478. https://doi.org/10.1007/s42952-020-00088-z
[9]	Xu, G.Y. and Bai, Y. (2020) Estimation of Nonparametric Additive Models with High Order Spatial Autoregressive Errors. Canadian Journal of Statistics, 49, 311-343.
[10]	Rice, J.A. and Wu, C.O. (2001) Nonparametric Mixed Effects Models for Unequally Sampled Noisy Curves. Biometrics, 57, 253-259. https://doi.org/10.1111/j.0006-341X.2001.00253.x
[11]	Karcher, P. and Wang, Y.D. (2001) Generalized Nonparametric Mixed Effects Models. Journal of Computational and Graphical Statistics, 10, 641-655. https://doi.org/10.1198/106186001317243377
[12]	Wang, H.X., Wu, Y.H. and Chan, E. (2017) Efficient Estimation of Nonparametric Spatial Models with General Correlation Structures. Australian & New Zealand Journal of Statistics, 59, 215-233. https://doi.org/10.1111/anzs.12192
[13]	Wang, Y.D. (2019) Smoothing Splines: Methods and Applications. Chapman and Hall/CRC, Boca Raton, FL.
[14]	Tait, A. and Woods, R. (2007) Spatial Interpolation of Daily Potential Evapotranspiration for New Zealand Using a Spline Model. Journal of Hydrometeorology, 8, 430-438. https://doi.org/10.1175/JHM572.1
[15]	Laslett, G.M. (1994) Kriging and Splines: An Empirical Comparison of Their Predictive Performance in Some Applications. Journal of the American Statistical Association, 89, 391-400. https://doi.org/10.1080/01621459.1994.10476759
[16]	Cressie, N. (2015) Statistics for Spatial Data. John Wiley & Sons, Hoboken.
[17]	Huang, C.F., Hsing, T. and Cressie, N. (2011) Nonparametric Estimation of the Variogram and Its Spectrum. Biometrika, 98, 775-789. https://doi.org/10.1093/biomet/asr056
[18]	Fernández-Casal, R., Castillo-Páez, S. and Francisco-Fernández, M. (2018) Nonparametric Geostatistical Risk Mapping. Stochastic Environmental Research and Risk Assessment, 32, 675-684. https://doi.org/10.1007/s00477-017-1407-y
[19]	Qadir, G.A. and Sun, Y. (2020) Semiparametric Estimation of Cross-Covariance Functions for Multivariate Random Fields. Biometrics, 77, 547-560. https://doi.org/10.1111/biom.13323
[20]	Yu, K.M., Mateu, J. and Porcu, E. (2007) A Kernel-Based Method for Nonparametric Estimation of Variograms. Statistica Neerlandica, 61, 173-197. https://doi.org/10.1111/j.1467-9574.2007.00326.x