Robust Estimation for Semi-Functional Linear Model with Autoregressive Errors

Abstract

It is well-known that the traditional functional regression model is mainly based on the least square or likelihood method. These methods usually rely on some strong assumptions, such as error independence and normality, that are not always satisfied. For example, the response variable may contain outliers, and the error term is serially correlated. Violation of assumptions can result in unfavorable influences on model estimation. Therefore, a robust estimation procedure of a semi-functional linear model with autoregressive error is developed to solve this problem. We compare the efficiency of our procedure to the least square method through a simulation study and two real data analyses. The conclusion illustrates that the proposed method outperforms the least square method, providing random errors follow the heavy-tail distribution.

1. Introduction

As a powerful tool in functional data analysis (FDA), functional regression modeling has been extensively studied in the past several decades. For example, Cardot et al. [1,2] introduced the functional linear model in which a scalar response variable is explained by a functional predictor. Yao et al. [3] proposed a functional linear regression model to analyse sparse longitudinal data. Hall and Horowitz [4] obtained the optimal convergence rates of estimators based on functional principal components analysis. Further, Li et al. [5] recently extended classical clusterwise linear regression to incorporate multiple functional predictors and proposed clusterwise functional linear regression models. For the functional nonparametric model, Ferraty and Vieu [6] extended the kernel regression to functional data. Baíllo and Grané [7] proposed a local linear regression and Burba et al. [8] used the k-Nearest Neighbour method to deal with the functional nonparametric models, respectively. Moreover, to improve the power of prediction or interpretation of the functional regression model, some functional semiparametric models have also been developed. For instance, Aneiros-Pérez and Vieu [9] first introduced a semi-functional partial linear regression model, and Aneiros-Pérez and Vieu [10] used this model to predict time series. Further, Shin [11] proposed new estimators of a partial functional linear model which explore the relationship between a scalar response variable and mixed-type predictors. Zhou and Chen [12] introduced a semi-functional linear model by combining the feature of a functional linear regression model and a nonparametric regression model.

It is well-known that the traditional functional regression models are mainly based on the least square or likelihood method. However, small proportions of outliers or the heavy-tailed nature of the error distribution could dramatically deteriorate the efficiency of estimators. To tackle this problem, many robust estimations have been introduced into the functional regression model. For example, Yu et al. [13] proposed a robust estimation method based on exponential square loss for the semi-functional linear regression model. Subsequently, Yu et al. [14] proposed a robust estimation procedure of a partial function linear model based on modal regression. Cai et al. [15] constructed a robust estimation with a modified Huber’s loss for partial functional linear models. Cao and Sun [16] studied the robust estimation of partial functional linear regression models based on functional principal component analysis and weighted compound quantile regression. Tang [17] introduced a robust estimation method for the semi-functional linear model, and Boente et al. [18] constructed new robust estimators for semi-functional linear regression models with a bounded loss function and a preliminary residual scale estimator. Recently, Boente and Daniela [19] considered the robust estimation for functional quadratic regression models. Moreover, Pannu and Billor [20] proposed a functional adaptive group LASSO variable selection method based on the weighted least absolute deviation. In the aforementioned literature, the random error is usually assumed to be independent and identically distributed. However, this assumption may be inappropriate in practice, especially when the observations are collected sequentially over time. However, as far as we know, only a few works have focused on robust estimation for regression models with serially correlated error. For example, Sanjoy et al. [21] considered robust estimation of nonlinear regression with autoregressive errors. Riazoshams et al. [22] studied the performance of a robust two-stage estimator in nonlinear regression with an autocorrelated error. Recently, Serenay and Baris [23] proposed a novel hybrid robust tapering approach in nonlinear regression with the presence of autocorrelation and outliers. Motivated by these previous works, this paper is mainly concerned with a robust semi-functional linear model with autoregressive errors (SFLAR). The rest of the paper is organized as follows. Section 2 introduces robust estimation and algorithm. Some simulation studies are presented in Section 3. Section 4 illustrates the good performance of the robust estimation method with two real data. Section 5 is the conclusion.

2. Model and Estimation

Suppose that ${\left\{x_i\left(t\right),y_i,z_i\right\}}_{i=1}^n$ satisfies the following semi-functional linear model with autoregressive errors (SFLAR)

$$\begin{array}{c}{y}_i={\int }_{\mathcal{T}}\beta \left(t\right)x_i\left(t\right)\mathrm{d}t+g\left(z_i\right)+{\epsilon }_i,\\ {\epsilon }_i={\displaystyle \sum _{l=1}^q}{a}_l{\epsilon }_{i-l}+e_i,\end{array}$$

(1)

where $y_i$ is a real-valued random response variable, and $z_i$ is a real-valued covariate defined on a compact interval $z=[\underline{z},\overline{z}]$. The functional covariate $x_i\left(t\right)$ defined on an interval $\mathcal{T}\subset R$ is a zero mean, second-order (i.e., ${\mathbb{E}\left|x\left(t\right)\right|}^2<\infty$ for all $t\in \mathcal{T}$ ) stochastic process defined on $(\mathrm{\Omega },B,P)$ with sample paths in the Hilbert space $L^2\left(\mathcal{T}\right)$. The inner product of $L^2\left(\mathcal{T}\right)$ is defined by $\langle u,v\rangle ={\int }_{\mathcal{T}}u\left(t\right)v\left(t\right)\mathrm{d}t$ for any $u,v\in L^2\left(\mathcal{T}\right)$ and norm $\parallel u\parallel ={\langle u,u\rangle }^{1/2}$. The unknown slope function $\beta \left(t\right)$ belongs to $L^2\left(\mathcal{T}\right)$. Without losing generality, throughout this paper, we assume $\mathcal{T}=[0,1]$. The nonparametric function $g\left(z\right)$ is an unknown smooth function; $e_i$ is a random error with zero mean, finite variance, and independent of $(x_i\left(t\right),z_i)$.

Since the unknown functions $\beta \left(t\right)$ and $g\left(z\right)$ are infinite-dimensional parameters, we cannot obtain their estimators directly via an optimization method. Therefore, some dimension reduction methods need to be applied. Specifically, suppose the covariance operator $C_x$ of the functional variable x is denoted by

$$\begin{array}{c}\hfill C_x(s,u)=Cov[x\left(s\right),x\left(u\right)],\end{array}$$

(2)

where $C_x(s,u)$ is continuous on ${[0,1]}^2$, and the eigenvalues of $C_x$, denoted by ${\left\{{\lambda }_j\right\}}_{j=1}^{\infty }$, satisfy ${\lambda }_1>{\lambda }_2>\cdots >0$. By Mercer’s Theorem, $C_x(s,u)$ can be expanded by

$$\begin{array}{c}\hfill C_x(s,u)=\sum _{j=1}^{\infty }{\lambda }_j{\varphi }_j\left(s\right){\varphi }_j\left(u\right),\end{array}$$

(3)

where ${\varphi }_j$ is the continuous orthogonal eigenfunction, and ${\lambda }_j$ is the corresponding eigenvalue. Obviously, the eigenfunction sequence ${\left\{{\varphi }_j\right\}}_{j=1}^{\infty }$ forms an orthonormal basis for $L^2[0,1]$. Then, random function $x_i\left(t\right)$ and the slope function $\beta \left(t\right)$ can be respectively expressed as

$$x_i\left(t\right)=\sum _{j=1}^{\infty }{\xi }_{ij}{\varphi }_j\left(t\right),\beta \left(t\right)=\sum _{j=1}^{\infty }{b}_j{\varphi }_j\left(t\right),$$

(4)

where ${\xi }_{ij}={\int }_{\mathcal{T}}{x}_i\left(t\right){\varphi }_j\left(t\right)\mathrm{d}t$ are called the principle component scores satisfying $\mathbb{E}{\xi }_{ij}=0$, and $\mathbb{E}{\xi }_{ij}{\xi }_{ik}={\lambda }_j\mathit{I}(j=k)$, and $b_j={\int }_{\mathcal{T}}\beta \left(t\right){\varphi }_j\left(t\right)\mathrm{d}t$.

Moreover, let $S_{r,N_n}$ be the space of polynomial splines defined on interval $[\underline{z},\overline{z}]$ with degree $r-1$ and knot sequence $\underline{z}<z_1<\cdots <z_{N_n}<\overline{z}$. The space $S_{r,N_n}$ is a K-dimensional linear space, $K=K\left(n\right)=N_n+r$. Following the arguments of de Boor [24], we can conclude that, if it is sufficiently smooth, the nonparametric function $g\left(z\right)$ can be approximately expressed as

$$g\left(z\right)\approx \sum _{k=1}^K{c}_k{B}_k\left(z\right),$$

(5)

where $B_k,k=1,2,\cdots ,K$ are B-spline basis functions in $S_{r,N_n}$. For more details on spline function, we refer to de Boor [24]. Substituting Equation (4) and Equation (5) into Equation (1) yields

$$y_i\approx \sum _{j=1}^J{\xi }_{ij}{b}_j+\sum _{k=1}^K{c}_k{B}_k\left(z_i\right)+{\epsilon }_i.$$

(6)

However, because the covariance function $C_x$ is unknown, eigenfunctions ${\varphi }_j$ are also unknown, and variables ${\xi }_{ij}$ are unobservable. To tackle this problem, we define the empirical version of covariance function $C_x$ by

$$\begin{array}{c}\hfill {\widehat{C}}_x(s,u)=\frac{1}{n}\sum _{i=1}^n{x}_i\left(s\right)x_i\left(u\right),\end{array}$$

(7)

define the estimated eigenelements of ${\widehat{C}}_x$ by

$$\begin{array}{c}\hfill \int {\widehat{C}}_x(s,u){\widehat{\varphi }}_j\left(u\right)\mathrm{d}u={\widehat{\lambda }}_j{\widehat{\varphi }}_j\left(s\right),(1<j<n).\end{array}$$

(8)

By Mercer’s Theorem, we have

$$\begin{array}{c}\hfill {\widehat{C}}_x(s,u)=\sum _{j=1}^n{\widehat{\lambda }}_j{\widehat{\varphi }}_j\left(s\right){\widehat{\varphi }}_j\left(u\right),\end{array}$$

(9)

where $({\widehat{\lambda }}_j,{\widehat{\varphi }}_j)$ are (eigenvalue, eigenfunction) pairs for the linear operator with kernel ${\widehat{C}}_x$, ordered such that ${\widehat{\lambda }}_1\ge {\widehat{\lambda }}_2>\cdots \ge 0$. We regard $({\widehat{\lambda }}_j,{\widehat{\varphi }}_j)$ as an estimator of $({\lambda }_j,{\varphi }_j)$. Moreover, Equation (6) can be written as

$$y_i\approx {\widehat{\mathit{U}}}_i^{\top }\mathit{b}+{\mathit{B}}_i^{\top }\mathit{c}+\sum _{l=1}^q{a}_l\left(y_{i-l}-{\widehat{\mathit{U}}}_{i-l}^{\top }\mathit{b}-{\mathit{B}}_{i-l}^{\top }\mathit{c}\right)+e_i,q+1\le i\le n.$$

(10)

Furthermore, denote

$$L_n(\mathit{a},\mathit{b},\mathit{c})=\frac{1}{n-q}\sum _{i=q+1}^n\rho \left(y_i-{\widehat{\mathit{U}}}_i^{\top }\mathit{b}-{\mathit{B}}_i^{\top }\mathit{c}-\sum _{l=1}^q{a}_l\left(y_{i-l}-{\widehat{\mathit{U}}}_{i-l}^{\top }\mathit{b}-{\mathit{B}}_{i-l}^{\top }\mathit{c}\right)\right),$$

(11)

where $\rho (·)$ is a proper robust loss function, $\mathit{a}={(a_1,\cdots ,a_q)}^{\top }$, $\mathit{b}={(b_1,\cdots ,b_J)}^{\top }$, $\mathit{c}={(c_1,\cdots ,c_K)}^{\top }$, ${\widehat{\mathit{U}}}_i={({\widehat{\xi }}_{i1},{\widehat{\xi }}_{i2},\cdots ,{\widehat{\xi }}_{iJ})}^{\top }$, and ${\mathit{B}}_i={\left(B_1\left(z_i\right),\cdots ,B_K\left(z_i\right)\right)}^{\top }$. Therefore, $\widehat{\mathit{a}}$, $\widehat{\mathit{b}}$ and $\widehat{\mathit{c}}$ can be obtained by solving the following minimization problem:

$$\underset{\mathit{a},\mathit{b},\mathit{c}}{min}{L}_n(\mathit{a},\mathit{b},\mathit{c}).$$

(12)

Finally, the estimator of functional coefficient $\beta \left(t\right)$ and the spline estimator of the nonparametric function $g\left(z\right)$ are given by $\widehat{\beta }\left(t\right)={\sum }_{j=1}^J{\widehat{b}}_j{\widehat{\varphi }}_j\left(t\right)$, $\widehat{g}\left(z_i\right)={\sum }_{k=1}^K{\widehat{c}}_k{B}_k\left(z_i\right)$, respectively.

Various loss functions $\rho (·)$ have been proposed for the robust estimation when outliers are present in the data. For example, the least absolute deviations (LAD) loss, that is, $\rho \left(r\right)=\left|r\right|$, can be regarded as a median regression. Although the LAD is not sensitive to outliers, its disadvantage is that the LAD is not smooth at the zero point, which will lead to some problems in model training. Therefore, considering the requirements of robustness and differentiability of the loss function, the Huber loss function (HB) [25] is considered in this paper, which is convex, smooth, and symmetric about 0. Hence, the optimization of Equation (12) for estimators is well-defined, and the local minimum problem is avoided.

To implement the proposed method, we take B-spline functions with equally spaced knots and the fixed degree three to approximate the unknown function $g\left(z\right)$. We only need to choose the numbers of eigenfunctions J and B-spline functions K. Many methods, such as the Akaike information criterion (AIC) and Bayesian information criterion (BIC), can be used to select the numbers of eigenfunctions J and B-spline functions K. In this paper, we choose the truncation parameters J and K by minimizing the following Bayesian information criterion (BIC) criteria:

$$\begin{array}{c}\hfill BIC\left(J,K\right)=log\left\{\sum _{i=q+1}^n\rho \left(y_i-\sum _{j=1}^J{\widehat{b}}_j{\widehat{\xi }}_{ij}-\sum _{k=1}^K{\widehat{c}}_k{B}_k\left(z_i\right)\right)\right\}+\frac{logn}{2n}\left(J+K\right),\end{array}$$

(13)

which is a function of J and K, the smaller the value of BIC, the better the fitting.

Furthermore, denote $\mathit{\theta }={({\mathit{b}}^{\top },{\mathit{c}}^{\top })}^{\top }$, ${\widehat{\mathit{H}}}_i^{\top }={({\widehat{\mathit{U}}}_i^{\top },{\mathit{B}}_i^{\top })}^{\top }$. Then, Equation (11) can be written as the following equation:

$$L_n(\mathbf{a},\mathit{\theta })=\frac{1}{n-q}\sum _{i=q+1}^n\rho \left(y_i-{\widehat{\mathit{H}}}_i^{\top }\mathit{\theta }-\sum _{l=1}^q{a}_l\left(y_{i-l}-{\widehat{\mathit{H}}}_{i-l}^{\top }\mathit{\theta }\right)\right).$$

(14)

Moreover, Equation (11) can be rewritten as:

$$L_n(\mathit{a},\mathit{\theta })=\frac{1}{n-q}\sum _{i=q+1}^n\rho \left(V_i-{\mathit{D}}_i^{\top }\mathit{\theta }\right),$$

(15)

where $V_i=y_i-{\sum }_{l=1}^q{a}_l{y}_{i-l}$, ${\mathit{D}}_i={\widehat{\mathit{H}}}_i-{\sum }_{l=1}^q{a}_l{\widehat{\mathit{H}}}_{i-l}^{\top }$. Then, the estimators of parameter vector $\mathit{a}={(a_1,a_2,\cdots ,a_q)}^{\top }$ and $\mathit{\theta }={({\mathit{b}}^{\top },{\mathit{c}}^{\top })}^{\top }$ can be obtained by minimizing Equation (15). However, the objective function $L_n(\mathit{a},\mathit{\theta })$ has no analytic solution for the unknown parameters. Therefore, this paper uses a two-step iteratively reweighted least-squares (TSIRLS) to obtain the local optimal solution $\widehat{\mathit{a}}$ and $\widehat{\mathit{\theta }}$ of Equation (15). The TSIRLS algorithm is summarized as follows:

Initialization: given initial value: ${\widehat{\mathit{\theta }}}^{\left(0\right)}$, ${\widehat{\mathit{a}}}^{\left(0\right)}$ using the LS estimators and let $k=0$.

Step 1. Update ${\widehat{\epsilon }}_i^{\left(k\right)}$, ${\mathit{W}}^{\left(k\right)}=diag\{w_1^{\left(k\right)},\cdots ,w_n^{\left(k\right)}\}$:

$$\begin{array}{c}\hfill {\widehat{\epsilon }}_i^{\left(k\right)}=y_i-{\mathit{H}}_i^{\top }{\widehat{\mathit{\theta }}}^{\left(k\right)},w_i^{\left(k\right)}=\frac{\psi \left({\widehat{\epsilon }}_i^{\left(k\right)}\right)}{{\widehat{\epsilon }}_i^{\left(k\right)}};\end{array}$$

Step 2. Update ${\tilde{\mathit{V}}}^{\left(k\right)}$, ${\tilde{\mathit{D}}}^{\left(k\right)}$:

$$\begin{array}{c}\hfill {\tilde{\mathit{V}}}^{\left(k\right)}=\left(\begin{array}{c}{\widehat{\epsilon }}_{q+1}^{\left(k\right)}\\ ⋮\\ {\widehat{\epsilon }}_n^{\left(k\right)}\end{array}\right),{\tilde{\mathit{D}}}^{\left(k\right)}=\left(\begin{array}{ccc}{\widehat{\epsilon }}_q^{\left(k\right)}& \cdots & {\widehat{\epsilon }}_1^{\left(k\right)}\\ ⋮& & ⋮\\ {\widehat{\epsilon }}_{n-1}^{\left(k\right)}& \cdots & {\widehat{\epsilon }}_{n-q}^{\left(k\right)}\end{array}\right);\end{array}$$

Step 3. Update ${\widehat{e}}_i^{\left(k\right)}$, ${\tilde{\mathit{W}}}^{\left(k\right)}=diag\{{\tilde{w}}_1^{\left(k\right)},\cdots ,{\tilde{w}}_n^{\left(k\right)}\}$:

$$\begin{array}{c}\hfill {\widehat{e}}_i^{\left(k\right)}={\widehat{\epsilon }}_i^{\left(k\right)}-\sum _{l=1}^q{\widehat{a}}_l^{\left(k\right)}{\widehat{\epsilon }}_{i-l}^{\left(k\right)},{\tilde{w}}_i^{\left(k\right)}=\frac{\psi \left({\widehat{e}}_i^{\left(k\right)}\right)}{{\widehat{e}}_i^{\left(k\right)}};\end{array}$$

Step 4. Update ${\widehat{\mathit{a}}}^{(k+1)}$:

$$\begin{array}{c}\hfill {\widehat{\mathit{a}}}^{(k+1)}={\left({\tilde{\mathit{D}}}^{\left(k\right)\top }{\tilde{\mathit{W}}}^{\left(k\right)}{\tilde{\mathit{D}}}^{\left(k\right)}\right)}^{-1}{\tilde{\mathit{D}}}^{\left(k\right)\top }{\tilde{\mathit{W}}}^{\left(k\right)}{\tilde{\mathit{V}}}^{\left(k\right)};\end{array}$$

Step 5. Update $V_i^{\left(k\right)}$, ${\mathit{D}}_i^{\left(k\right)}$:

$$\begin{array}{c}{V}_i^{\left(k\right)}=y_i-\sum _{l=1}^q{a}_l^{(k+1)}{y}_{i-l},{\mathit{D}}_i^{\left(k\right)}={\widehat{\mathit{H}}}_i-\sum _{l=1}^q{a}_l^{(k+1)}{\widehat{\mathit{H}}}_{i-l};\end{array}$$

Step 6. Update ${\widehat{\mathit{\theta }}}^{(k+1)}$:

$$\begin{array}{c}{\widehat{\mathit{\theta }}}^{(k+1)}={\left({\mathit{D}}^{\left(k\right)\top }{\mathit{W}}^{\left(k\right)}{\mathit{D}}^{\left(\mathbf{k}\right)}\right)}^{-1}{\mathit{D}}^{\left(k\right)\top }{\mathit{W}}^{\left(k\right)}{\mathit{V}}^{\left(k\right)},\end{array}$$

where

$$\begin{array}{c}\hfill {\mathit{V}}^{\left(k\right)}=\left(\begin{array}{c}{V}_{q+1}^{\left(k\right)}\\ ⋮\\ V_n^{\left(k\right)}\end{array}\right),{\mathit{D}}^{\left(k\right)}=\left(\begin{array}{c}{\mathbf{D}}_{q+1}^{\left(k\right)}\\ ⋮\\ {\mathit{D}}_n^{\left(k\right)}\end{array}\right);\end{array}$$

Step 7. Update $k=k+1$. Repeat the above steps until $\parallel {\widehat{\mathit{\theta }}}^{(k+1)}-{\widehat{\mathit{\theta }}}^{\left(k\right)}{\parallel }_1+{\parallel {\widehat{\mathit{a}}}^{(k+1)}-{\widehat{\mathit{a}}}^{\left(k\right)}\parallel }_1\le {10}^{-8}$. Then, the estimation of the parameter vectors $\widehat{\mathit{a}}={({\widehat{a}}_1,{\widehat{a}}_2,\cdots ,{\widehat{a}}_q)}^{\top }$, $\widehat{\mathit{b}}={({\widehat{b}}_1,{\widehat{b}}_2,\cdots ,{\widehat{b}}_J)}^{\top }$, and $\widehat{\mathit{c}}={({\widehat{c}}_1,{\widehat{c}}_2,\cdots ,{\widehat{c}}_K)}^{\top }$ are obtained.

3. Simulation

In this section, we conduct some simulation studies to investigate the finite sample performance of our proposed method. The simulation data $\{(x_i\left(t\right),y_i,z_i),1\le i\le n\}$ are generated from the following SFLAR model:

$$\begin{array}{c}\hfill y_i={\int }_{\mathcal{T}}\beta \left(t\right)x_i\left(t\right)dt+g\left(z_i\right)+{\epsilon }_i,{\epsilon }_i=\sum _{l=1}^q{a}_l{\epsilon }_{i-l}+e_i,\end{array}$$

(16)

where $z_i$ is distributed uniformly on $[-1,1]$, and the functional predictor $x_i\left(t\right)={\sum }_{j=1}^{50}{\xi }_{ij}{\varphi }_j\left(t\right)$, ${\xi }_{ij}$ is distributed as independent normal with mean 0 and variance ${\lambda }_j={\left((j-0.5)\pi \right)}^{-1}$ and ${\varphi }_j\left(t\right)=\sqrt{2}sin\left((j-0.5)\pi t\right)$. For the unknown functions $\beta (·)$ and $g(·)$ and autoregressive coefficients $(a_1,\cdots ,a_q)$, given $\beta \left(t\right)=\sqrt{2}sin\left(0.5\pi t\right)+3\sqrt{2}sin\left(1.5\pi t\right)$, $g\left(z\right)=z+3cos\left(2\pi z\right)$, the first-order autoregressive error $AR\left(1\right):{\epsilon }_t=0.8{\epsilon }_{t-1}+e_t$. We evaluate the performance of various loss functions with different sample sizes $n=100,300$ and 500. Each experiment has been repeated 1000 times. Similar to Cai et al. [15], the following four scenarios error distributions are considered:

 (1) 

$e_i\sim N(0,1)$;

 (2) 

$e_i\sim \pi N(0,1)+(1-\pi )N(0,100)$, where $\pi$ = 0.8;

 (3) 

$e_i\sim t\left(2\right)$;

 (4) 

$e_i\sim \mathit{Cauchy}(0,1)$.

Moreover, according to one of the reviewers’ suggestions, we also consider different values of $\pi$ in Scenario 2. Specifically, the values of $\pi$ are taken as 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 and 0.9.

To evaluate the performances of our estimators, we define the integrated mean square errors (IMSE) of the estimators of function $\beta \left(t\right)$, $g\left(z\right)$,

$$\begin{array}{c}\hfill {\mathrm{IMSE}}_1={\int }_0^1{[\widehat{\beta }\left(t\right)-\beta \left(t\right)]}^2\mathrm{d}t,{\mathrm{IMSE}}_2={\int }_0^1{[\widehat{g}\left(z\right)-g\left(z\right)]}^2\mathrm{d}z,\end{array}$$

(17)

and the mean square prediction error (MSPE):

$$\begin{array}{c}\hfill \mathrm{MSPE}=\frac{1}{n}\sum _{i=q+1}^n{\left(\widehat{y_i}-g\left(z_i\right)-{\int }_0^1\beta \left(t\right)x_i\left(t\right)\mathrm{d}t\right)}^2,\end{array}$$

(18)

where ${\widehat{y}}_i=\widehat{g}\left(z_i\right)+{\int }_{\mathcal{T}}\widehat{\beta }\left(t\right)x_i\left(t\right)\mathrm{d}t$. Further, in order to obtain $\widehat{\beta }\left(t\right)$ and $\widehat{g}\left(z_i\right)$, the two tuning parameters J and K are selected by minimizing the BIC value. Specifically, the range of J is selected within 1 to 6, and the range of K is selected within 3 to 6.

We compare the finite sample performance of the SFLAR model with ordinary least squares (OLS), least absolute deviations (LAD) loss, and Huber (HB) loss function. Further, according to one of the reviewers’ suggestions, we also consider the exponential square (EXP) loss function (Yu et al. [13]). To implement the robust method considered in Yu et al. [13], the tuning parameters h, J and K need to be selected. In this paper, the tuning parameters are chosen by cross-validation and BIC. Specifically, for some fixed h, the tuning parameters J and K are selected by BIC. Then, under the selected J and K, like Jiang [26], the tuning parameter h is chosen by $5-$fold cross-validation (CV) within the possible grids points $h=0.5\times 1.{02}^l$ ($l=0,1,\cdots ,100$). We repeat the process until the tuning parameters are invariant.

Table 1. The Bias and MSE of $\widehat{a}$ under different error distributions.

n	Methods	$\mathit{N}(0,1)$		MN		$\mathit{t}\left(2\right)$		$\mathit{Cauchy}(0,1)$
		Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE
	OLS	−0.0106	0.0043	−0.0084	0.0040	−0.0081	0.0037	−0.0225	0.0038
100	LAD	−0.1424	0.0264	−0.0613	0.0068	−0.1036	0.0164	−0.0232	0.0017
	EXP	−0.0240	0.0059	−0.0025	0.0008	−0.0084	0.0021	−0.0007	0.0300
	HB	−0.0106	0.0044	−0.0019	0.0006	−0.0049	0.0019	−0.0011	0.0002
	OLS	−0.0051	0.0014	−0.0067	0.0012	−0.0051	0.0012	−0.0097	0.0011
300	LAD	−0.0400	0.0033	−0.0140	0.0004	−0.0226	0.0012	−0.0026	0.0000
	EXP	−0.0143	0.0019	−0.0016	0.0001	−0.0041	0.0005	0.0000	0.0094
	HB	−0.0053	0.0014	−0.0021	0.0002	−0.0031	0.0005	−0.0002	0.0000
	OLS	−0.0045	0.0008	−0.0042	0.0007	−0.0017	0.0007	−0.0053	0.0005
500	LAD	−0.0234	0.0015	−0.0002	0.0001	0.0099	0.0004	−0.0011	0.0000
	EXP	−0.0119	0.0010	−0.0006	0.0001	−0.0018	0.0002	−0.0001	0.0058
	HB	−0.0044	0.0008	−0.0009	0.0001	−0.0008	0.0002	−0.0002	0.0000

reports the bias and mean square error (MSE) of estimator $\widehat{a}$ under different error distributions.

Table 2. The mean (SD) of ${\mathrm{IMSE}}_1$ for $\widehat{\beta }\left(t\right)$ under different error distributions and sample sizes.

n	Methods	$\mathit{N}(0,1)$	MN	$\mathit{t}\left(2\right)$	$\mathit{Cauchy}(0,1)$
		Mean (SD)	Mean (SD)	Mean (SD)	Mean (SD)
	OLS	0.2664 (0.1266)	1.4068 (1.2305)	0.8923 (4.2805)	-    (-)
100	LAD	0.3298 (0.1663)	0.4360 (0.2541)	0.4730 (0.3153)	0.4996 (0.3083)
	EXP	0.3243 (0.1714)	0.4247 (0.7672)	0.3786 (0.2029)	0.4602 (0.2534)
	HB	0.2711 (0.1305)	0.3963 (0.2124)	0.3800 (0.2211)	0.4734 (0.2721)
	OLS	0.1099 (0.0445)	0.4753 (0.3015)	0.4589 (2.0454)	-    (-)
300	LAD	0.1272 (0.0526)	0.1870 (0.0788)	0.1438 (0.0638)	0.1964 (0.0830)
	EXP	0.1663 (0.0712)	0.1787 (0.0745)	0.1792 (0.0759)	0.1982(0.0815)
	HB	0.1110 (0.0451)	0.1813 (0.0764)	0.1357 (0.0569)	0.1986 (0.0805)
	OLS	0.0836 (0.0305)	0.3248 (0.1985)	0.7908 (17.3387)	-    (-)
500	LAD	0.0936 (0.0352)	0.0928 (0.0346)	0.1026 (0.0417)	0.1557 (0.0573)
	EXP	0.1394 (0.0537)	0.1469 (0.0553)	0.1463 (0.0546)	0.1570(0.0583 )
	HB	0.0844 (0.0308)	0.1474 (0.0550)	0.0980 (0.0380)	0.1567 (0.0584)

and

Table 3. The mean (SD) of ${\mathrm{IMSE}}_2$ for $\widehat{g}(·)$ under different error distributions and sample sizes.

n	Methods	$\mathit{N}(0,1)$	MN	$\mathit{t}\left(2\right)$	$\mathit{Cauchy}(0,1)$
		Mean (SD)	Mean (SD)	Mean (SD)	Mean (SD)
	OLS	0.3509 (0.3710)	7.1999 (8.9383)	7.0279 (-)	-    (-)
100	LAD	0.4821 (0.7918)	1.4814 (2.1051)	1.1504 (4.6860)	1.9271 (2.9633)
	EXP	0.4132 (0.4437)	1.0579 (2.0157)	0.7482 (0.8314)	1.2962 (1.5780)
	HB	0.3713 (0.3992)	0.8917 (1.0224)	0.7717 (0.8362)	1.3913 (1.7224)
	OLS	0.1140 (0.1166)	2.2063 (2.6667)	2.0490 (12.9143)	-    (-)
300	LAD	0.1597 (0.1726)	0.3016 (0.3629)	0.2406 (0.2751)	0.3699 (0.4413)
	EXP	0.1385 (0.1356)	0.2434 (0.2624)	0.2229 (0.2740 )	0.3681 (0.4213)
	HB	0.1189 (0.1242)	0.2388 (0.2597)	0.2074 (0.2528)	0.3629 (0.4087)
	OLS	0.0677 (0.0754)	1.3708 (1.5933)	3.2699 (-)	-    (-)
500	LAD	0.1043 (0.1177)	0.1046 (0.1065)	0.1410 (0.1513)	0.1945 (0.2207)
	EXP	0.0914 (0.0891)	0.1414 (0.1428)	0.1376 (0.1440)	0.2129 (0.2366)
	HB	0.0722 (0.0787)	0.1485 (0.1636)	0.1266 (0.1406)	0.2085 (0.2365)

present the sample mean and standard deviation (SD) of ${\mathrm{IMSE}}_1$ and ${\mathrm{IMSE}}_2$ under different error distributions.

Table 4. The mean (SD) of MSPE under different error distributions and sample sizes.

n	Methods	$\mathit{N}(0,1)$	MN	$\mathit{t}\left(2\right)$	$\mathit{Cauchy}(0,1)$
		Mean (SD)	Mean (SD)	Mean (SD)	Mean (SD)
	OLS	0.3588 (0.3542)	7.3662 (8.8633)	7.3858 (-)	-    (-)
100	LAD	0.4677 (0.4338)	1.4368 (1.7128)	0.9803 (0.9559)	1.9208 (2.7565)
	EXP	0.4426 (0.4447)	1.0267 (1.6935)	0.7628 (0.8031)	1.3085 (1.3842)
	HB	0.3792 (0.3800)	0.8862 (0.9175)	0.7665 (0.7987)	1.4059 (1.6041)
	OLS	0.1266 (0.1173)	2.3559 (2.6711)	2.1143 (12.5808)	-    (-)
300	LAD	0.1767 (0.1725)	0.3310 (0.3619)	0.2613 (0.2739)	0.4033 (0.4421)
	EXP	0.1607 (0.1349)	0.2703 (0.2626)	0.2502 (0.2712)	0.4017 (0.4206)
	HB	0.1317 (0.1246)	0.2662 (0.2597)	0.2274 (0.2509)	0.3969 (0.4090)
	OLS	0.0775 (0.0751)	1.4681 (1.5899)	3.6659 (-)	-    (-)
500	LAD	0.1165 (0.1175)	0.1172 (0.1067)	0.1563 (0.1508)	0.2207 (0.2209)
	EXP	0.1105 (0.0889)	0.1635 (0.1429)	0.1601 (0.1448)	0.2399 (0.2368)
	HB	0.0821 (0.0785)	0.1706 (0.1644)	0.1409 (0.1411)	0.2353 (0.2365)

lists the MSPE results of model prediction under different error distributions for different estimation methods. The symbol “-” indicates that these values are much higher than others, so they are not listed in the table. Figure 1 displays the sample means of $\mathrm{IMSE}={\mathrm{IMSE}}_1+{\mathrm{IMSE}}_2$ for functional parameters $\beta \left(t\right)$, $g\left(z\right)$ and MSE ($\times {10}^3$) for autoregressive coefficient a when $n=500$, and $\pi$ takes values in the set $\{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9\}$ for Scenario (2).

From the simulation results of different estimation methods under different error distributions in Table 1, Table 2, Table 3 and Table 4 and Figure 1, we can conclude that:

$\left(\mathbf{1}\right)$ For the autoregressive coefficient vector a, the estimation results of all methods based on OLS, LAD, EXP and HB seem equally good, whenever the error distribution is normal or heavy-tailed, which indicates the estimator $\widehat{a}$ is not sensitive to the distribution of random error. At the same time, we can see that the values of Bias and MSE of HB and EXP are less than OLS. Hence, it can be concluded the robust estimators based on HB and EXP outperform the OLS method.

$\left(\mathbf{2}\right)$ For the slope function $\beta \left(t\right)$ and nonparametric function $g\left(z\right)$, the OLS estimator behaves badly when the error distribution is non-normal, particularly when the random error is heavy-tailed. In contrast, the estimators of LAD, EXP and HB are much better than the OLS. Even if the expectation and variance of the error do not exist, the robust estimation method is still feasible.

$\left(\mathbf{3}\right)$ Whenever the distribution of random error is normal or heavy-tailed, the bias and MSE of $\widehat{a}$, mean and SD of ${\mathrm{IMSE}}_1$, ${\mathrm{IMSE}}_2$ and $\mathrm{MSPE}$ based on LAD, EXP and HB methods decrease as the sample size increases. Hence, the proposed robust procedure is effective.

4. Applications to Real Data

In this section, we illustrate the proposed robust estimation procedure for a semi-functional linear model with autoregressive errors by two real data analyses. The first example is the electricity consumed data, and the second one is the Tecator data.

4.1. Electricity Consumption Data

This subsection aims to compare the proposed robust estimation method with the OLS method, that is, whether the proposed robust estimation method is as good as the OLS estimation method in the absence of outliers in the data set. For this purpose, we consider electricity consumption data, which is available on the website http://www.eia.doe.gov/emeu/aer, (accessed on 31 December 2022). The data set includes the electricity consumption $c_i$ (KWH) from January 1973 to January 2016 (517 months) (see Figure 2) and annual average retail price $z_i$ (per KWH, including tax) (43 years).

As shown in Figure 2, the time series obviously shows some linear trends and some heterogeneity in the variance structure. Similar to Aneiros and Vieu [10] and Yu et al. [27], we eliminate the heteroscedasticity and the linear trend of the electricity retail sales data by differencing the log(electricity data) and obtain the time series (see Figure 3): ${\mathcal{X}}_i=log\left(c_{i+1}\right)-log\left(c_i\right),i=1,2,\cdots ,516$. Let $x_i\left(s\right)=\left\{{\mathcal{X}}_{12(i-1)+s},i=1,2,\cdots ,43,s=1,2,\cdots ,12\right\}$ be the monthly log (difference electricity data).

There is only one observation per month, so we use 10 cubic B-spline basis functions to convert the discrete monthly electricity data into functional continuous smooth annual curve predictors $\{x_i\left(t\right),t\in [1,12],i=1,2,\cdots ,43\}$ (see Figure 4). The response variable is $y_i\left(s\right)={\mathcal{X}}_{12i+s}$, which represents the electricity consumption data in the s-month of the i-year, where $i=1,2,\cdots ,42,s=1,2,\cdots ,12$. The electricity data set is divided into the training sample ${\left\{y_i\left(s\right),z_i,x_i\left(t\right)\right\}}_{i=1}^{41}$ and the test sample $\left\{y_{42}\left(s\right),z_{42},x_{42}\left(t\right)\right\}$.

To ascertain whether there are outliers in the data set, the boxplot of the electricity consumption data from January to December is displayed in Figure 5. As shown in Figure 5, there are some outliers in March and November. Hence, the robust method may perform better than OLS. Meanwhile, to determine whether the random errors are serially correlated, we first display the autocorrelation function plot (see Figure 6) of the residual sequences ${\left\{{\widehat{\epsilon }}_i\left(s\right)\right\}}_{i=1}^{41}$, $s=1,\cdots ,12$ of the following semi-functional linear model (SFLM) based on the HB estimation procedure:

$$y_i\left(s\right)={\int }_1^{12}\beta \left(t\right)x_i\left(t\right)\mathrm{d}t+g\left(z_i\right)+{\epsilon }_i\left(s\right),i=1,2,\cdots ,41;s=1,\cdots ,12.$$

(19)

From Figure 6, it can be seen that there exist serial correlations for random error sequences ${\epsilon }_i\left(8\right)$ and ${\epsilon }_i\left(9\right)$. We further use the Ljung–Box (LB) test statistics to test whether first-fourth order serial correlations exist. The corresponding $p-$values are 0.0049 and 0.0013, respectively. Under the given significance level of 0.05, we can also conclude that the random error sequences ${\epsilon }_i\left(8\right)$ and ${\epsilon }_i\left(9\right)$ are the fourth-order serial correlation, which is consistent with Figure 6. Hence, we consider using the robust method to estimate the following SFLAR model based on the training set:

$$y_i\left(s\right)={\int }_1^{12}\beta \left(t\right)x_i\left(t\right)\mathrm{d}t+g\left(z_i\right)+{\epsilon }_i,{\epsilon }_i=\sum _{l=1}^q{a}_l{\epsilon }_{i-l}+e_i,i=1,2,\cdots ,41,$$

(20)

and then the test set is used to predict the values of $y_{42}\left(s\right)$.

We also compare the robust estimation procedure (LAD, EXP and HB) with OLS to illustrate the proposed mehtod. Further, to implement our approach, the unknown function $g(·)$ is approximated by a cubic B-spline function with equispaced knots. The truncation parameters J and K are selected by the BIC method and calculated according to Equation (13) in Section 2. Similar to Aneiros and Vieu [10] and Yu et al. [27], the following two error criteria, mean quadratic error (MQE):

$$\begin{array}{c}\hfill \mathrm{MQE}=\frac{1}{12}\sum _{s=1}^{12}{\left(y_{42}\left(s\right)-{\widehat{y}}_{42}\left(s\right)\right)}^2,\end{array}$$

(21)

and mean relative quadratic error (MRQE):

$$\begin{array}{c}\hfill \mathrm{MRQE}=\frac{1}{12}\sum _{s=1}^{12}\frac{{\left(y_{42}\left(s\right)-{\widehat{y}}_{42}\left(s\right)\right)}^2}{Var\left(y\right(s\left)\right)}\end{array}$$

(22)

are used to evaluate the prediction ability of the model, where $Var\left(y\right(s\left)\right)$ is the empirical variance of ${\left\{y_i\left(s\right)\right\}}_{i=1}^{41}$.

Table 5. MQE and MRQE of three estimation methods under the SFLAR model.

Criterions	Methods
	OLS	LAD	HB	EXP
MQE	0.00017	0.00016	0.00017	0.00017
MRQE	0.34483	0.27929	0.31347	0.34435

shows the MQE and MRQE of LAD, HB, EXP and OLS estimation. It can be seen from Table 5 that the robust estimation procedure performs better than OLS when there are some outliers.

4.2. Tecator Data

In this subsection, we further apply the robust procedure to analyze the Tecator spectral data set, which originated from http://lib.stat.cmu.edu/datasets/tecator, (accessed on 31 December 2022). This data set consists of 215 meat samples. Each sample records a spectrometric curve corresponding to a spectrum of absorbance measured at 100 wavelengths between 850 and 1050 nm by the Near-Infrared Transmission (NIT) principle, and the contents of moisture, protein and fat measured in percentages. Chen et al. [28], Aneiros and Vieu [29], Lin et al. [30] and Sang et al. [31] have conducted extensive research on these data, aiming to predict the content of certain components (such as fat) by using spectral absorbance. Recently, Cai et al. [15] predicted fat content (Y) according to water content ($Z_1$), protein content ($Z_2$) and spectral curve $X\left(t\right)$.

This paper is interested in whether the accuracy of LAD, EXP and HB is better than OLS in predicting fat content (y) when the distribution of the response is not normal, and the random errors are correlated. For this purpose, we first use the Shapiro–Wilk statistic to test whether the fat content distribution is normal. The $p-$value of the test statistic is 0.000, which indicates that the distribution is not normal. Further, the skewness value (0.796) shows that the data have a right-skewed distribution. The least squares estimation may not be effective, so the robust method is adopted. Then, like in Section 4.1, we apply the LB test statistics to test whether there are first-order serial correlations based on the HB method. The corresponding $p-$value is 0.0222, which is less than the given significance level of 0.05. Therefore, we also apply the robust method to estimate the following SFLAR model:

$$y_i={\int }_{850}^{1050}\beta \left(t\right)x_i\left(t\right)\mathrm{d}t+g\left(z_i\right)+{\epsilon }_i,{\epsilon }_i=a_1{\epsilon }_{i-1}+e_i,i=1,2,\cdots ,215,$$

(23)

where the response variable $y_i$ is the fat content of the i-th sample, and the explanatory variables $z_i$ and $x_i\left(t\right)$ stand for the protein and the spectrometric curve of the i-th sample, respectively.

Similar to Cai et al. [15], the two criteria, mean absolute prediction error (MAPE):

$$\begin{array}{c}\hfill \mathrm{MAPE}=\frac{1}{214}\sum _{i=2}^{215}\left|y_i-{\widehat{y}}_i\right|,\end{array}$$

(24)

and the median prediction error (MPE):

$$\begin{array}{c}\hfill \mathrm{MPE}=\underset{i\in \{2,3,\cdots ,215\}}{\mathrm{median}}\left|y_i-{\widehat{y}}_i\right|\end{array}$$

(25)

are used to evaluate the prediction effect of different methods. From the values of MAPE and MPE of LAD, HB, EXP and OLS displayed in

Table 6. MAPE and MPE of three estimation methods under SFLAR.

Criterions	Methods
	OLS	LAD	HB	EXP
MAPE	1.862	1.771	1.783	1.891
MPE	1.528	1.345	1.388	1.138

, it can be seen that the robust methods have better performance than the OLS estimation when the distribution of the response is not normal.

5. Conclusions

This paper considers the estimation problem of a semi-functional linear model with autoregressive errors when the random error follows a heavy-tailed distribution, or there are some outliers. The traditional estimation methods, such as the least square or likelihood method, are susceptible to outliers or heavily tailed errors, resulting in the corresponding estimation no longer being effective. Therefore, this paper introduces the robust estimation procedure. Simulation research and two real data analyses show that if the random error follows the normal distribution, the robust estimation has the same good statistical properties as the OLS. However, if the random error is heavily tailed, the robust estimation outperforms the OLS method.