Evaluation of parameters for the sustainable supply chain management: a Taiwanese fresh-fruit sector

1.   Introduction

Functional data analysis (FDA) (e.g., ^[1]) has drawn considerable attention over recent years, owing to a great deal of flexibilities and universal applications in handling high-dimensional data sets. A fundamental and important tool for FDA is functional linear models.

There are a lot of researches in literature on the inference of functional linear models and their extensions, see, among others, ^[2,3,4] for earlier works, and ^{[5,6,7,8,9,10]} for recent works. As is well known, in the estimation of regression models, the choice of loss function is essential to obtain a highly efficient and robust estimator. Most of earlier works employed the square loss function and obtained ordinary least squares (OLS) estimators. In recent years, many other loss functions have been considered in the estimation of functional linear models and their extensions. Kato ^[6], Tang and Cheng ^[11] studied the quantile regression (QR) with functional linear models and partial functional linear models, respectively. Yu et al. ^[12] proposed a robust exponential squared loss estimation procedure (ESL) and established the asymptotic properties of the proposed estimators. Cai et al. ^[13] introduced a new robust estimation procedure by employing a modified Huber function, whose tail function is replaced by the exponential squared loss (H-ESL) in the partial functional linear model.

It is well known that the square loss function pays attention to reflect the distributional features of the entire distribution, whereas QR, ESL and H-ESL methods focus on bounding the influence of outliers when the data are heavy-tailed, respectively. Thus, developing a method, which can both reflect distributional features and bound outliers effectively, is highly desirable in data analysis. We note that, in the context of principal component analysis (PCA), Lim and Oh ^[14] proposed a new approach using a weighted linear combination of asymmetric Huber loss functions to demonstrate the distributional features of data as well as keep robust to outliers. The asymmetric Huber loss functions are defined as

with $c^* = 1.345$, and $\tau \in (0, 1)$ being a parameter to control the degree of skewness. The function $\rho_{\tau}(\cdot)$ is equivalent to the Huber loss function (see, ^[15]) when $\tau$ is equal to 0.5 and is most exactly the same as the quantile loss function when $c^*$ is small enough.

Motivated by the appealing characteristics of the asymmetric Huber functions, in this paper, we first investigate a new asymmetric Huber regression (AHR) estimation procedure to analyse skewed data for the partial functional linear model, based on the functional principal component analysis. To improve the estimation accuracy for single AHR estimation, we develop a weighted composite asymmetric Huber regression (WCAHR) estimation by combining the strength across multiple asymmetric Huber regression models. A practical algorithm for WCAHR estimators based on pseudo data is developed to implement the estimation method. The asymptotic properties of the proposed estimators are also derived. Extensive simulations are carried out to show the superiority of the proposed estimators.

Finally, we apply the proposed methods to two data sets. In the first example, we analyze the electricity data. Figure 1 presents the estimated density of the residuals and the residual diagnostic plot obtained by fitting the model (4.1) in Section 4.1 via the OLS method. The distribution of the residuals is skewed, bimodal, and there are some outliers in the dataset. Given that the WCAHR can effectively manage such data, we use the proposed method to conduct an analysis to this data set. Another example in Section 4 considers the Tecator data set. Similarly, Figure 2 presents the density of the residuals and the residual diagnostic plot obtained by fitting the model (4.2) in Section 4.2 via the OLS method, which demonstrates that the distribution of the residuals is skewed and far from normality. Undoubtedly, WCAHR regression is applicable to analyzing this data set on account of its appealing features.

To our knowledge, it is the first to discuss the asymmetric Huber regression problems under functional models framework. The proposed WCAHR method possesses advantages that include the robustness to outliers as well as reflecting the relationships between potential explanatory variables and the entire distribution of response. It retains the advantages in analysing skewed data and the obtained estimators rely on the shape of the entire distribution rather than merely on the data nearby a specific quantile level or skewness level of the asymmetric Huber loss, thereby avoiding the limitations of these methods. These advantages are revealed by both theoretical conclusions and numerical results. The relevant algorithm is data-adaptive, and capable of reflecting the distributional features of the data without prior information, and is robust to outliers.

The rest of this paper is organized as follows. In Section 2, we formally describe the estimation procedures, and develop a new algorithm. We also establish the asymptotic behaviors of the proposed estimators as well as a list of technical assumptions needed in the theoretical research. In Section 3, the finite performances of the proposed estimators are evaluated through simulations. Section 4 illustrates the use of the proposed methods in the analyses of electricity data and Tecator data. Brief conclusions on the proposed methods are made in Section 5. All technique proofs are provided in Section A.

2.   Methodology and main results

2.1.   Proposed methods

Let $Y$ be a real value random variable, $\mathit{\boldsymbol{Z}} = {(Z_1, \cdots, Z_p)^{\rm T}}$ be a $p$ dimensional random vector with zero mean and finite second moment. Let $\{X(t): t \in \mathscr{T}\}$ be a zero mean, second-order stochastic process with sample paths in $L_{2}(\mathscr{T})$, which consists of square integrable functions with inner product $\langle x, y\rangle = \int_{\mathscr{T}}x(t)y(t)dt$ and norm $\|x\| = \langle x, x\rangle^{1/2}$, respectively, here $\mathscr{T}$ is a bounded closed interval. Without loss of generality, we suppose $\mathscr{T} = [0, 1]$ throughout the paper. The dependence between $Y$ and $(X, \mathit{\boldsymbol{Z}})$ is expressed by the partial functional linear regression as following,

Here, random error $e$ is assumed to be independent of $\mathit{\boldsymbol{Z}}$ and $X$, $\mathit{\boldsymbol{\alpha}} = \big(\alpha_{1}, \cdots, \alpha_{p}\big)^{\rm T}$ is an unknown $p$-dimensional parameter vector, and the slope function $\beta(\cdot)$ is an unknown square integrable function on $[0, 1]$.

Let $(\mathit{\boldsymbol{Z}}_i, X_i(\cdot), Y_i), i = 1, \cdots, n,$ be independent observations generated by model ($2.1$) and let $e_i = Y_i-\mathit{\boldsymbol{Z}}_i^{\rm T} \mathit{\boldsymbol{\alpha}}-\int_0^1\beta(t)X_i(t)dt, i = 1, \cdots, n$. The covariance and empirical covariance functions for $X(\cdot)$ are defined as $c_X(t, s) = {\rm{Cov}}(X(t), X(s))$, $\hat{c}_X(t, s) = \frac{1}{n}\sum_{i = 1}^nX_i(t)X_i(s)$ respectively. Based on the Mercer's Theorem, $c_X$ and $\hat{c}_X$ can be represented as following,

where $\lambda_1 > \lambda_2 > \cdot\cdot\cdot > 0 \quad {\rm{and}}\quad\hat{\lambda}_1\geq\hat{\lambda}_2\geq\cdots\geq\hat{\lambda}_{n+1} = \cdot\cdot\cdot = 0$ are each the ordered eigenvalue sequences of the covariance operator $C_X$ and its estimator $\hat{C}_X$ with kernels $c_X$ and $\hat{c}_X$, which are defined by $C_Xf(s) = \int_{0}^{1} c_X(t, s)f(t)dt$ and $\hat{C}_Xf(s) = \int_{0}^{1} \hat{c}_X(t, s)f(t)dt$ with $C_X$ being assumed strictly positive, and $\{v_i(\cdot)\}$ and $\{\hat{v}_i(\cdot)\}$ are the corresponding orthonormal eigenfunction sequences. Besides, $(\hat{v}_i(\cdot), \hat{\lambda}_i)$ is treated as an estimator of $(v_i(\cdot), \lambda_i)$.

Similarly, we can define ${c}_{YX}(\cdot) = {\rm {Cov}}(Y, X(\cdot))$, ${c}_{\mathit{\boldsymbol{Z}}} = {\rm {Var}}(\mathit{\boldsymbol{Z}}) = {E[{{\mathit{\boldsymbol{Z}}}}{{\mathit{\boldsymbol{Z}}}}^{\rm T}]}$, ${c}_{\mathit{\boldsymbol{Z}}Y} = {\rm {Cov}}(\mathit{\boldsymbol{Z}}, Y)$, ${c}_{\mathit{\boldsymbol{Z}}X}(\cdot) = {\rm {Cov}}(\mathit{\boldsymbol{Z}}, X(\cdot)) = ({c}_{{Z}_1X}(\cdot), \cdots, {c}_{{Z}_pX}(\cdot))^{\rm T}$. And the corresponding empirical counterparts defined below can be used as their estimators,

By the Karhunen-Loève representation, $X_i(t)$ and $\beta(t)$ can be expanded into

here $\gamma_j = \langle \beta(\cdot), v_j(\cdot)\rangle = \int_{0}^{1}\beta(t)v_j(t)dt$, and $\xi_{ij} = \langle X_i(\cdot), v_j(\cdot)\rangle$.

Owing to the orthogonality of $\{v_1(\cdot), ..., v_m(\cdot)\}$ and Eq ($2.2$), Model ($2.1$) can be transformed into

where ${\mathit{\boldsymbol{U}}}_i = ({\xi}_{i1}, ..., {\xi}_{im})^{\rm T}$, ${\mathit{\boldsymbol{\gamma}}} = (\gamma_1, ..., \gamma_m)^{\rm T}$, $\widetilde{e}_i = \sum_{j = m+1}^{\infty}\gamma_j\xi_{ij}+{e}_i$, and the tuning parameter $m$ may increase with the sample size $n$.

Replacing $v_j(\cdot)$ with its estimator $\hat{v}_j(\cdot)$, the $\tau$th AHR estimators $\bar{\mathit{\boldsymbol{\alpha}}}$ and $\bar{\beta}(t) = \sum_{j = 1}^{m}\bar{\gamma}_j\hat{v}_j(t)$ can be obtained by minimizing the loss function over $\mathit{\boldsymbol{\alpha}}$, $\mathit{\boldsymbol{\gamma}}$ and $b_\tau$ as follows:

where the asymmetric Huber loss function $\rho_{\tau}(u)$ is defined in ($1.1$), and $\hat{\mathit{\boldsymbol{U}}}_i = (\hat{\xi}_{i1}, \cdots, \hat{\xi}_{im})^{\rm T}$ with $\hat{\xi}_{ij} = \langle X_i(\cdot), \hat{v}_j(\cdot)\rangle$. Here the true value of $b_{\tau}$ is defined as the solution that minimizes the loss function ${E}\left\{\rho_{\tau}(e-\theta) \right\}$ over $\theta \in \mathbb{R}$, and we call it the $\tau$th number of $e$ with respect to (1.1).

Remark 1. In model (2.1), we suppose the intercept term is zero. In fact, if there is an intercept, we then may absorb it into the distribution of $e$. Thus, the main impact of model (2.1) is finding the contribution of the predictors to the response, and the zero mean assumption for $e$ is not needed.

Noting that the regression coefficients are the same across different skewness asymmetric Huber regression models, and being inspirited by ^[14] and ^[16], we combine the strength across multiple asymmetric Huber regression models and propose a WCAHR method. Specifically, the WCAHR estimators $\hat{\mathit{\boldsymbol{\alpha}}}$ and $\hat{\beta}(t) = \sum_{j = 1}^{m}\hat{\gamma}_j\hat{v}_j(t)$ can be obtained by minimizing the following loss function with respect to $(\mathit{\boldsymbol{b}}, \mathit{\boldsymbol{\alpha}}, \mathit{\boldsymbol{\gamma}})$:

where $\left\{\tau_{k}\right\}_{k = 1}^{K}$ are predetermined levels over (0, 1), $b_k = b_{\tau_k}$ for brevity, $\mathit{\boldsymbol{b}} = \left(b_1, \cdots, b_K\right)^{\rm T}$, and the weights $w_{1}, \ldots, w_{K}$, which control the contribution of each loss function, are positive constants satisfying $\sum_{k} w_{k} = 1$.

Remark 2. Generally speaking, one can choose the equidistant levels as $\tau_{k} = k /(K+1)$ for $k = 1, 2, \ldots, K$ for a given $K$, similar to what often has been done in composite quantile regression. Although one can also apply data-adaptive methods, such as cross validation, to select $K$, we do not pursue this topic here. As for the weights, we consider two choices. The first is using the equal weights $w_{1} = \cdots = w_{K} = 1/K$. The obtained estimators are called composite asymmetric Huber regression (CAHR) estimators. As the second choice in this study, to preferably describe the distribution information of the data, we consider a $K$-dimensional weight vector $(w_{1}, \cdots, w_{K}) = \big({f}({b}_{01}), \cdots, {f}({b}_{0K})\big)\Big/\sum_{{k} = 1}^K{f}({b}_{0k})$, where $\mathit{\boldsymbol{b}}_0 = (b_{01}, \cdots, b_{0K})$ is the true value vector of $\mathit{\boldsymbol{b}}$, and $f(\cdot)$ is the density function of the random error. In practice, we estimate $f(\cdot)$ by kernel density estimation method.

Denote $\mathcal{S} = \big\{(Y_i, \mathit{\boldsymbol{Z}}_i, X_i(\cdot)):1\leq i\leq n\big\}$, and given a new copy of $(\mathit{\boldsymbol{Z}}, X)$, namely the predictor variables $(\mathit{\boldsymbol{Z}}_{n+1}, {X_{n+1}(\cdot)})$, once we gain the estimated $\mathit{\boldsymbol{\alpha}}$ and $\beta(t)$, the mean squared prediction error (MSPE) can be obtained, take asymmetric Huber regression for example,

where $\mathit{\boldsymbol{\alpha}}_0$ and $\beta_0(t)$ are the true values of $\mathit{\boldsymbol{\alpha}}$ and $\beta(t)$, respectively. The MSPEs of CAHR and WCAHR have analoguos definitions, and denoted by ${\rm{MSPE}}_{CAHR}$ and ${\rm{MSPE}}_{WCAHR}$, respectively.

2.2.   Computational algorithm

Note that the minimization problems for AHR and CAHR estimators are special cases of WCAHR method. Here, we just present the practical algorithm for WCAHR based on pseudo data. A similar argument can be found in ^[14] to implement the data-adaptive principal component analysis. The algorithm is simply described as following.

Step 1. Given initial estimators $\hat{\boldsymbol{\alpha}}^{(0)}$ and $\hat{\boldsymbol\gamma}^{(0)}$ for $\boldsymbol{\alpha}_{0}$ and $\gamma_{0}$, respectively.

Step 2. Iterate, until convergence, following these three steps for $L = 0, 1, \ldots$

(a) Compute the residuals as $\hat{e}_{i}^{(L)} = Y_{i}-{{\mathit{\boldsymbol{Z}}}}_{i}^{\rm T} \hat{\boldsymbol\alpha}^{(L)}-\hat{{\mathit{\boldsymbol{U}}}}_{i}^{\rm T} \hat{\boldsymbol\gamma}^{(L)}$. (b) Calculate the empirical pseudo data vector $\mathit{\boldsymbol{G}}^{(L)} = ({G}_1^{(L)}, \cdots, {G}_n^{(L)})^{\rm T}$ in the element-wise way, $G_{i}^{(L)} = {{\mathit{\boldsymbol{Z}}}}_{i}^{\rm T} \hat{\boldsymbol\alpha}^{(L)}+\hat{{\mathit{\boldsymbol{U}}}}_{i}^{\rm T} \hat{\boldsymbol\gamma}^{(L)}+ \sum_{k = 1}^{K} w_{k} \psi_{\tau_{k}}\big(\hat{e}_{i}^{(L)}-\hat{b}_k^{(L)}\big),$ for given weights $\left(w_{1}, \cdots, w_{K}\right)$ and $\hat{b}_k^{(L)} = {\rm argmin}_\mu\sum_{i = 1}^n \rho_{\tau_k}\left(\hat{e}_{i}^{(L)}-\mu\right)$ at each $k$. Here $\psi_{\tau_k}(u) = \rho'_{\tau_k}(u) = ({\tau_k}-1)I(u < -c^*)+\frac{(1-{\tau_k})}{c^*}uI(-c^*\leq u < 0)+\frac{\tau_k}{c^*}uI(0\leq u < c^*)+\tau_kI(u\geq c^*)$. (c) Obtain next iterative estimates $\hat{\boldsymbol\alpha}^{(L+1)}$ and $\hat{\boldsymbol\gamma}^{(L+1)}$ by using the OLS method for response variable $\tilde{Y}_{i} = G_{i}^{(L)}$ and covariates ${{\mathit{\boldsymbol{Z}}}}_{i}, \hat{{\mathit{\boldsymbol{U}}}}_{i}$.

2.3.   Asymptotic properties

In this section, we will establish the asymptotic properties of the estimators defined in the previous section. We shall first present some notations, suppose $\mathit{\boldsymbol{\gamma}}_0 = (\gamma_{01}, ..., \gamma_{0m})^{\rm T}$ is the true values of $\mathit{\boldsymbol{\gamma}}$, $F(\cdot)$ is the cumulative distribution function of the random error. In addition, the notation $\|\cdot\|$ represents the $\mathcal{L}^2$ norm of a function or the Euclidean norm of a vector, and $a_n\sim b_n$ indicates that $a_n/b_n$ is bounded away from zero and infinity as $n \rightarrow \infty$. For simplicity, in this paper, $C$ represents a generic positive constant whose value may change from line to line. Next, to obtain the asymptotic properties, some technical assumptions are listed as follows.

C1. The random process $X(\cdot)$ and score $\xi_j = \langle X(\cdot), v_j(\cdot)\rangle$ satisfy the following condition: ${\rm {E}}\|X(\cdot)\|^4 < \infty$ and ${\rm{E}}[\xi_j^4]\leq C\lambda_j^2, \ j \geq 1$.

C2. The eigenvalues of $C_X$ and the score coefficients fulfil the conditions below:

${(a)}$ There exist constants $C$ and $a > 1$ such that $C^{-1}j^{-a}\leq\lambda_j\leq Cj^{-a}, \lambda_j-\lambda_{j+1}\geq Cj^{-a-1}, \ j \geq 1$;

${(b)}$ There exist constants $C$ and $b > a/2+1$ such that $|\gamma_j|\leq Cj^{-b}, \ j \geq 1$.

C3. The random vector $\mathit{\boldsymbol{Z}}$ satisfies ${\rm {E}}\|\mathit{\boldsymbol{Z}}\|^4 < \infty$.

C4. There is some constant $C$ such that $\mid \langle {c}_{ {Z}_lX}, v_j\rangle \mid \leq Cj^{-(a+b)}, \; l = 1, \cdots, p, \ j \geq 1$.

C5. Let $\eta_{il} = {Z}_{il}-\langle g_l, X_i \rangle$ with $g_l = \sum\limits_{j = 1}^\infty\lambda_j^{-1}\langle { c}_{{Z}_lX}, v_j\rangle v_j$, and $\boldsymbol \eta_i = (\eta_{i1}, \cdots, \eta_{ip})^{\rm T}$, then ${\boldsymbol\eta}_{1}, \cdots, {\boldsymbol\eta}_{n}$ are i.i.d random vectors. We further assume that ${{\rm {E}}[{\boldsymbol\eta}_{i}|{X_i}] = 0, \; \; {\rm {E}}[{\boldsymbol \eta_i}{\boldsymbol \eta_i}^{\rm T}| X_i]{\overset{a.s.}{ = }}\mathit{\boldsymbol{\Sigma}}}$ is a constant positive definite matrix.

C6. $b_\tau$ is the unique solution of ${\rm {E}}\left[\rho_{\tau}^{\prime}\left(e-b_{\tau}\right)\right] = 0$, and $h_{\tau}(b_{\tau}) = (1-\tau)f(b_{\tau}-c^*)+\frac{1-\tau}{c^*}(F(b_{\tau})-F(b_{\tau}-c^*)) +\frac{\tau}{c^*}(F(b_{\tau}+c^*)-F(b_{\tau}))+\tau f(b_{\tau}+c^*)$ is continuous at $b_{\tau}$. Furthermore, we suppose that $h_{\tau}(b_{\tau}) > 0$.

C6'. $b_{0k}$ is the unique solution of $E\left[\rho_{\tau_{k}}^{\prime}\left(e_i-b_{0 k}\right)\right] = 0$, $h_{\tau_k}(b_{0k}) = (1-\tau_k)f(b_{0k}-c^*)+\frac{1-\tau_k}{c^*}(F(b_{0k})-F(b_{0k}-c^*)) +\frac{\tau_k}{c^*}(F(b_{0k}+c^*)-F(b_{0k}))+\tau_kf(b_{0k}+c^*)$ is continuous at $b_{0k}$, $k = 1, \cdots, K$. Furthermore, there exist some positive constants $C_1, \; C_2$ such that $0 < C_1\leq\underset{1\leq k\leq K}\min h_{\tau_k}(b_{0k})\leq\underset{1\leq k\leq K}\max h_{\tau_k}(b_{0k})\leq C_2 < +\infty$.

C1 is the commonly used condition for establishing the consistency of the the empirical covariance operator of $X$ in functional linear model and partial functional regression model. For example, it has been adopted in ^[3,17,18] (mean regression), ^[6,11] (quantile regression), ^[12,13] (robust estimation procedure), among others. C2(a) is used to identify the slope function $\beta(t)$ via preventing the spacings between eigenvalues being too small, and C2(b) ensures the sufficiently smooth of slope function $\beta(t)$. C3–C5 are needed to deal with the vector-type covariate in the model ($2.1$) (see ^[19]). More concretely, C3 is for the asymptotic behaviors of $\hat{{c}}_{{{\mathit{\boldsymbol{Z}}}}X}$ and $\hat{{c}}_{{{\mathit{\boldsymbol{Z}}}}}$. C4 is used to ensure the effect of truncation on the estimator of ${\beta}(\cdot)$ be sufficiently small. C5 is a commonly used condition in the literature on partial functional regression models (see for example, ^[4,20,21]). The assumptions on ${\rm {E}}[{\boldsymbol\eta}_{i}|{X_i}]$ and ${\rm {E}}[{\boldsymbol \eta_i}{\boldsymbol \eta_i}^{\rm T}|X_i]$ are slightly strong, and are used to fix the identifiability of ${\boldsymbol \alpha}$ and simplify the proof of the theorems. It is easy to see that $\langle g_l, X_i \rangle$ is the projection of $Z_l$ onto $X_i$, and ${{\rm{E}}}({\boldsymbol \eta_i}) = 0$, ${\rm Cov}\left(\eta_{il}, \langle g_l, X_i \rangle\right) = 0$, ${\rm E}\|\boldsymbol \eta_i\|^4 < \infty$ even without the assumptions. The facts can be used to obtain similar results to the following theorems with more complicated technics (see, for example, ^[6]) and more conditions to ensure the identifiability. Other type conditions on ${\boldsymbol \eta_i}$ can be found in ^[11,22]. C6 and C6$^\prime$ are specific to the AHR and WCAHR (CAHR) cases respectively, which are primarily used to ensure the asymptotic behaviors of our estimators.

The following theorems discuss the convergence rate of the estimated ${\beta}(\cdot)$, the asymptotic normality of the estimated ${\boldsymbol\alpha}$ and the convergence rate of the mean squared prediction error. To obtain this, we further assume $(\mathit{\boldsymbol{Z}}_{n+1}, {X_{n+1}(\cdot)})$ is independent of $\mathcal{S}$ in this paper.

The next theorem establishes the large sample properties of the AHR estimators.

Theorem 1. Suppose that the Conditions C1–C6 are satisfied, and the tuning parameter $m\sim n^{1/(a+2b)}$, then

where $\stackrel{d}{\longrightarrow}$ represents convergence in distribution, and $V = \mathit{E}\left[\psi_{\tau}^2\left(e-b_{\tau}\right)\right]$ with $\psi_{\tau}(u) = \rho'_{\tau}(u) = ({\tau}-1)I(u < -c^*)+\frac{(1-{\tau})}{c^*}uI(-c^*\leq u < 0)+\frac{\tau}{c^*}uI(0\leq u < c^*)+\tau I(u\geq c^*)$.

The asymptotic properties of the proposed WCAHR estimators are presented in the following theorem.

Theorem 2. Under the Conditions C1–C5 and C6$'$, if the tuning parameter is taken as $m\sim n^{1/(a+2b)}$, then

where ${{\mathit{\boldsymbol{w}}}} = \left(w_{1}, \ldots, w_{K}\right)^{T}$ and $\mathit{\boldsymbol{V}} = \left({V_{kl}}\right)_{1 \leq k, \; l \leq K}$, here $V_{k l} = \mathit{E}\left[\psi_{\tau_{k}}\left(e-b_{0 k}\right)\psi_{\tau_{l}}\left(e-b_{0 l}\right)\right]$ with $1 \leq k, l \leq K$.

Remark 3. The results illustrate that the slope function estimator has the same convergence rate as the estimators in ^[6] and ^[4], which are optimal in the minimax sense. Note that it is similar to quantile regression that no moment condition on error term is needed here. In addition, we notice that the rate attained in predicting ${Y}_{n+1}$ is faster than the rate attained in estimating $\beta(t)$. Trace its root and use $MSPE_{AHR}$ as an example, it is for the integral operator providing additional smoothness in computing $\int^{1}_{0}\bar{\beta}(t)X_{n+1}(t)dt$ from $\bar\beta(t)$.

Remark 4. If all $w_{k}$s are equal, then Theorem 2.2 reduces to the asymptotic properties of the CAHR estimators. Taking $\tau_{1} = \tau$, it is easy to see that there is a weight vector ${{\mathit{\boldsymbol{w}}}}$ such that $\frac{{{\mathit{\boldsymbol{w}}}}^{T} {{\mathit{\boldsymbol{V}}}} {{\mathit{\boldsymbol{w}}}}}{\left\{\sum_{k = 1}^{K}w_{k}h_{\tau_{k}}\left(b_{0 k}\right)\right\}^2} < \frac{V}{\left\{h_{\tau}\left(b_{01}\right)\right\}^2}$. Note that the right hand side of the inequality is just the asymptotic variance given in Theorem 1.

3.   Simulation study

In this section, a Monte Carlo simulation is used to investigate the finite sample properties of the proposed estimation approaches. The data sets used in the simulation are generated from the following model,

where the slope function $\beta(t) = \sqrt 2 \sin(\pi t/2)+3 \sqrt 2 \sin(3 \pi t/2)$, and $X(t) = \sum\limits_{j = 1}^{50} \xi_j \phi_j(t)$, here $\phi_j(t) = \sqrt 2 \sin((j-0.5)\pi t)$, and $\xi_j$s are mutually independent normal random variables with mean 0 and variance $\lambda_j = ((j-0.5)\pi)^{-2}$. The true values of parameters are set as $\mathit{\boldsymbol{\alpha}} = (\alpha_1, \alpha_2)^{\rm T} = (10, 5)^{\rm T}$, and ${\mathit{\boldsymbol{Z}}}\sim N\left(0, \Sigma_{{\mathit{\boldsymbol{Z}}}}\right)$ with $\left(\Sigma_{{\mathit{\boldsymbol{Z}}}}\right)_{i, j} = 0.75^{|i-j|}$ for $i, j = 1, 2$.

Five different distributions for $e$ are considered as follows: (a) standard normal distribution $N(0, 1)$; (b) positively skewed normal distribution ${SN}(0, 1, 15)$; (c) positively skewed $t$-distribution ${St}(0, 1, 5, 3)$; (d) mixture of normals $(MN)$ $0.95 N(0, 1)+0.05 N\left(0, 10^{2}\right)$, which produces a distribution with outliers of response; (e) bimodal distribution ($Bimodal$) ${\tilde{\eta}}N(-1.2, 1) + (1-\tilde{\eta})N(1.2, 1)$ with $\tilde{\eta}\sim Binomial(1, 0.5)$. The multiplier $\sigma(\mathit{\boldsymbol{Z}}, X)$ can be generated from either of the following two models:

(A) (homoscedastic) $\sigma(\mathit{\boldsymbol{Z}}, X) = 1$;

(B) (heteroscedastic) $\sigma(\mathit{\boldsymbol{Z}}, X) = { \left|1+0.1\left({Z_{1}\alpha_1^*+Z_{2}\alpha_2^*+\int_0^1 X(t)\beta^*(t) dt}\right)\right|}$, where $\alpha_1^* = \alpha_2^* = 1$, and $\beta^*(t) = \sqrt 2 \sin(\pi t/2)+\sqrt 2 \sin(3 \pi t/2)$.

Implementing the proposed estimation method requires the predetermined levels over (0, 1), i.e., $\left\{\tau_{k}\right\}_{k = 1}^{K}$. Similar to the setting in ^[14], we take $K = 19$, and choose the equidistant levels $\tau_{k} = k /(K+1)$, $k = 1, 2, \ldots, K$. In addition, for the WCAHR estimator, we employ the adaptive weights given in Remark 2.

For comparison, we also calculate the OLS estimator, the least absolute deviation (LAD) estimator, the ESL estimator, the H-ESL estimator, the Huber estimator (which corresponds to the case of AHR estimator at $\tau = 0.5$), the CAHR estimator, and the AHR estimators at $\tau = 0.25$ and $0.75$. In this study, the sample size $n$ is set as 200 or 400.

To implement these methods, we need to choose the tuning parameter $m$. In this paper, $m$ is selected by the cumulative percentage of total variability (CPV) method, that is,

where $\delta$ equals $85\%$. Other criterion, such as AIC, BIC, can be employed.

For each setting and different methods, bias (Bias), standard deviation (Sd) of the estimated $\alpha_1$ and $\alpha_2$, and the mean squared error (MSE) of the estimated $\mathit{\boldsymbol{\alpha}}$ with ${\rm{MSE}} = \frac{1}{S}\sum_{s = 1}^S\sum_{j = 1}^2(\hat{\alpha}_j^s-{\alpha}_j)^2$, as well as the mean integrated squared error (MISE) of the estimated $\beta(t)$ over S = 500 repetitions are summarized, where ${\rm{MISE}} = {\left\{\frac{1}{100S}\sum_{s = 1}^S\sum_{i = 1}^{100} {\left(\hat{\beta}^s(t_i)-{\beta}(t_i)\right)}^2\right\}}$ with $t_i$s being 100 equally spaced grids in [0, 1], here $\hat{\alpha}_j^s$, $\hat{\beta}^s(\cdot)$ are the estimates of $\hat{\alpha}_j$ and $\hat{\beta}(\cdot)$ from the $s$th sampling, $j = 1, 2$.

Table 1 presents the results in the homoscedastic case. From Table 1, we can see the following facts: (a) The Sd, MSE and MISE decrease as the sample size $n$ increases from $200$ to $400$. (b) The proposed estimators are almost unbiased, which further illustrates the consistency combining with the fact (a). (c) The proposed adaptively weighted estimator performs similarly to the OLS estimator under the normal error, and is comparable to the corresponding H-ESL estimator for the mixture of normal distributions, but is significantly better than the other estimators considered when the distribution of model error is skewed or bimodal, and still enjoys the favoured being robust to outliers. This demonstrates that the proposed WCAHR estimator can well adapt to different error distributions, thus is more useful in practice.

Table 2 presents the results in the more challenged heteroscedastic case, which violates the condition in this paper. The proposed WCAHR estimator outperforms the other estimators considered for the normal, skewed normal and skewed $t$ error distributions, and is comparable to the corresponding H-ESL estimator for the mixture of normal distribution and bimodal distribution. This further illustrates that the proposed WCAHR estimator may be more applicable. Although the simulation results show the appealing performance for the considered heteroscedastic errors, general theoretical results are still challenging, and more conditions on the conditional moments of $e$ may be helpful.

In order to detect the effect of the level choice to the performance of the WCAHR estimators, especially for the skewed error distributions, we also change in the simulation the number $K_1$ of levels over (0, 0.5) given the total level number $K$. Specifically, for the given $K = 19$ and different values of $K_1$, we set $\tau_i = \frac{i}{2K_1}, \; \text{for}\; i = 1, \cdots, K_1$ and $\tau_i = \frac{K+1-2K_1+i}{2(K+1-K_1)}, \; \text{for}\; i = K_1+1, \cdots, K.$ Table 3 presents the estimation results. We find from the results that the choice of the levels does not destroy the performance of the WCAHR estimators, although less number of levels over (0, 0.5) leads to slightly larger MSE and MISE for the positively skewed error distributions. In addition, the MISEs and MSEs decrease as $K_1$ increases, and stabilize eventually. This may motivate that one can take more levels appropriately over (0, 0.5) in dealing with the positively skewed error distributions.

4.   Real data application

In this section, we use the proposed estimation methods to the Electricity data and the Tecator data set, and the competing methods mentioned in Section 3. In the applications, all the observations are centralized before the regression analysis.

4.1.   Electricity data

The data set consists of the daily average hourly electricity spot prices of the German electricity market ($Y$), the hourly values of Germany's wind power infeed ($X(t)$), the precipitation height ($Z_1$) and the sunshine duration ($Z_2$). Here we only consider the working days span from January 1, 2006 to September 30, 2008. The hourly values of Germany's wind power infeed curves are shown in the left panel of Figure 3. The data set can be obtained from the online supplements of Liebl ^[23]. Now we adopt the following partial functional linear regression model to fit the data:

Firstly, the OLS method is applied to fit model (4.1). Then Shapiro-Wilk test is applied to test the normality of the residuals and the $p$-value is less than $2.2 \times 10^{-16}$. In addition, we also give the estimated density of the residuals and the residual diagnostic plot (see Figure 1). Both the test and plot clearly indicate that $e$ follows a skewed distribution with outliers. Notice that the density of the residuals is similar to the error distribution ${Bimodal}$ discussed in Section 3, and the simulation results illustrate that the proposed method can be applied and provide reliable inference for this kind of data.

To evaluate the predictions obtained with different estimation methods, we randomly divide data set into a training sample of size 478 subjects and a testing sample with the remaining 160 subjects (indexed by $\mathcal{J}$). The data are split for $N = 100,200,400$ times, respectively. We use the median quadratic error of prediction ($\rm{MEDQEP}$) defined below as the criterion to evaluate the performances:

The left 3 columns of Table 4 present the MEDQEPs of different methods mentioned above. According to the results of calculation, the WCAHR method is uniformly superior to the other estimators.

Table 5 (the first 2 columns) presents the estimates of the parametric part by the estimation methods based on the whole data set. According to the results, both the precipitation height and the sunshine duration have negative effects on the daily average hourly electricity spot prices. In addition, Figure 4(a) plots the estimated slope function obtained by the WCAHR method, the estimates for $\beta(\cdot)$ obtained by other methods mentioned above exhibit similar patterns and thus omitted here. From the figure, we can see the prices have a larger (in the absolute value) linkage with the wind power infeed in the daytime, which reflects the economic phenomena of price sensitivity, and more specifically, the market is active during the daytime and thus there is more correlation between the prices and the wind power infeed in the daytime. Secondly, the Germany's wind power infeed has a negative effect on the daily average hourly electricity spot prices, which reflects supply-demand balance, that is, more wind infeed creates the oversupply of electricity and thus reduces the price.

4.2.   Tecator data

The Tecator data set consists of 215 meat samples. For each sample, moisture, fat and protein are recorded in percent, and a 100-channel spectrum of absorbances is measured by the spectrometer. The data set is available from the R package fda.usc (see ^[24]). The right panel of Figure 3 shows the spectral trajectories. In this paper, the objective is to investigate the effects of the spectral trajectories $X(t)$, water content $Z_1$ and protein content $Z_2$ on the fat content $Y$ by fitting the following model:

The density of the residuals and the residual diagnostic plot in Figure 2 illustrate the error follows a skewed distribution with outliers. Similarly, to assess the prediction accuracy, the 215 meat samples are randomized into training set with 180 subjects and testing set with 35 subjects. We also randomly split the data set for $N = 100,200,400$ times and use MEDQEP as criterion to evaluate the finite sample performances of different estimation procedures. The comparison results are shown in the right 3 columns of Table 2, from which we know the proposed method performs better than the competing estimation procedures in view of prediction accuracy.

The estimated coefficients $\hat{\alpha}_{1}$, $\hat{\alpha}_{2}$ using various methods based on the whole data set are also shown in the last 2 columns of Table 3. Both the protein content and water content have negative effects on the fat content. Next, the right panel of Figure 4 demonstrates the estimated slope function curves based on the WCAHR method. It is obvious that the spectrum curve of absorbance has negative impact on the quantity of fatty. In addition, the estimated slope functions by other methods mentioned above exhibit similar patterns and thus omitted here.

5.   Conclusions

In this paper, we study the WCAHR estimation in the partial functional linear regression model. We use the functional principal component basis to approximate the functions, and obtain the estimators of the unknown parameter vector and the slope function through minimizing the weighted asymmetric Huber loss function. The asymptotic normality of the estimated parameter vector and the convergence rate of the estimated slope function are presented.

The proposed approach is designed for automatically reflecting distributional features as well as bounding outliers effectively without requiring prior information of the data. Simulation results show that the proposed method is almost as efficient as OLS when the error follows a normal distribution, while keeps robust to outliers and still works well when the error follows skewed or bimodal distributions. That is to say, the method is adaptive to the distribution of the error in the regression model. The analyses of two examples further illustrate that the utility of the proposed methods in modelling and forecasting.

The novelty of the method is that it focuses on the extraction of major features as well as shielding the estimator from outliers. The proposed WCAHR estimation procedure can be extended to more general situations, including dependent functional data, sparse modeling, partially observed functional data, and high dimension setting. In addition, we still need to try objective way to select $K$.

Acknowledgments

Authors are highly grateful to anonymous referees and editor in chief for their valuable comments and suggestions for improving this paper. This research was supported by National Natural Science Foundation of China (Grant No. 11771032 and No. 11971045), Natural Science Foundation of Beijing (Grant No. 1202001), and Natural Science Foundation of Shanxi Province, China (Grant No. 20210302124530).

Conflict of interest

The authors declare that they have no competing interests.

A.   Proofs of the Theorems

We just prove Theorem 2. The Theorem 1 is a special case of Theorem 2.

Proof of Theorem 2:

Let $\delta_n = \sqrt{\frac{m}{n}}$, ${P}_{n} = \frac{1}{n}\sum_{i = 1}^{n}\mathit{\boldsymbol{\eta}}_i\mathit{\boldsymbol{\eta}}_i^{\rm T}$, $\mathit{\boldsymbol{V}}_n = \delta_n^{-1}{P}_{n}^{{ 1/2}}(\mathit{\boldsymbol{\hat{\alpha}}}-\mathit{\boldsymbol{\alpha}}_0)$, $\tilde{\boldsymbol\eta}_{i} = P_{n}^{-1 / 2} {\boldsymbol\eta}_{i}$, $\Lambda = \operatorname{diag}\left\{\lambda_{1}, \cdots, \lambda_{m}\right\}$, ${{{\mathit{\boldsymbol{A}}}}}_i = {\Lambda}^{-{ 1/2}}{{{\mathit{\boldsymbol{U}}}}_i}$, ${\hat{{\mathit{\boldsymbol{A}}}}}_i = {\Lambda}^{-{ 1/2}}\hat{{{\mathit{\boldsymbol{U}}}}_i}$, $H_{m} = \left(\lambda_{1}^{-1}\langle {c}_{{\mathit{\boldsymbol{Z}}}, X}, v_{1}\rangle, \ldots, \lambda_{m}^{-1}\langle {c}_{{\mathit{\boldsymbol{Z}}}, X}, v_{m}\rangle\right)^{\mathrm{T}}$, $\mathit{\boldsymbol{W}}_n = \delta_n^{-1}\Lambda^{ 1/2}[(\mathit{\boldsymbol{\hat{\gamma}}}-\mathit{\boldsymbol{\gamma}}_0)+H_{m}(\mathit{\boldsymbol{\hat{\alpha}}}-\mathit{\boldsymbol{\alpha}}_0)]$, $r_i = \int_{0}^{1}\beta_0(t)X_i(t)dt-{\hat{\mathit{\boldsymbol{U}}}}_i^{\rm T}\mathit{\boldsymbol{\gamma}}_0$, $B_i = H_{m}^{\rm T}({{\mathit{\boldsymbol{U}}}}_i-\hat{{\mathit{\boldsymbol{U}}}}_i)+\sum\limits_{j = m+1}^\infty\lambda_j^{-1}\langle {c}_{\mathit{\boldsymbol{Z}}X}, v_j\rangle \xi_{ij}$, $\tilde{B}_i = \delta_n^{-1}(\mathit{\boldsymbol{\hat{\alpha}}}-\mathit{\boldsymbol{\alpha}}_0)^{\rm T}B_i$, $S_{nk} = \delta_{n}^{-1}\left(\hat{b}_{k}-b_{0 k}\right)$, $\boldsymbol{S}_{n} = \left(\sqrt{w_1}S_{n 1}, \ldots, \; \sqrt{w_K}S_{n K}\right)^{\rm T}$, $\mathcal{F}_n = \Big\{(\mathit{\boldsymbol{V}}_n, \mathit{\boldsymbol{W}}_n, \mathit{\boldsymbol{S}}_{n}):\big\|(\mathit{\boldsymbol{V}}_n^{\rm T}, \mathit{\boldsymbol{W}}_n^{\rm T}, \mathit{\boldsymbol{S}}_{n}^{\rm T})^{\rm T}\big\|\leq L\Big\}$, $T_n = \left\{\left(\mathit{\boldsymbol{Z}}_1, X_1(\cdot)\big), ..., \big(\mathit{\boldsymbol{Z}}_n, X_n(\cdot)\right)\right\}$.

Next, we will show that, for any given $\eta > 0$, there exists a sufficiently large constant $L$, such that

This means that there is a local minimizer $\mathit{\boldsymbol{\hat{\alpha}}}$, $\hat{\mathit{\boldsymbol{b}}}$ and $\mathit{\boldsymbol{\hat{\gamma}}}$ in the ball $\Big\{(\mathit{\boldsymbol{V}}_n, \mathit{\boldsymbol{W}}_n, \mathit{\boldsymbol{S}}_{n}):\big\|(\mathit{\boldsymbol{V}}_n^{\rm T}, \mathit{\boldsymbol{W}}_n^{\rm T}, \mathit{\boldsymbol{S}}_{n}^{\rm T})^{\rm T}\big\|\leq L\Big\}$ such that $\|\mathit{\boldsymbol{\hat{\alpha}}}-\mathit{\boldsymbol{\alpha}}_0\| = O_p(\delta_n)$, $\left|\hat{b}_{k}-b_{0 k}\right| = O_p(\delta_n)$, $\|\Lambda^{1/2}(\mathit{\boldsymbol{\hat{\gamma}}}-\mathit{\boldsymbol{\gamma}}_0)\| = O_p(\delta_n)$ with a probability at least $1-\eta$.

First, by $\|v_j-\hat{v}_j\|^2 = O_p(n^{-1}j^2)$ (see, ^[4]), one has

For $\rm{D}_1$, using Conditions C1, C2, and the Hölder inequality, one can obtain

As for $\rm{D}_2$, due to $\mathrm{E}\Big \{\sum_{j = m+1}^{ \infty}\langle X_i, v_j\rangle\gamma_{0j}\Big \} = 0, \; \mathrm{Var}\Big \{\sum_{j = m+1}^{ \infty}\langle X_i, {v}_j\rangle\gamma_{0j}\Big \} = \sum_{j = m+1}^{ \infty}\lambda_j{\gamma_{0j}}^2 \leq C\sum_{j = m+1}^{ \infty}j^{-(a+2b)} = O(n^{-\frac{a+2b-1}{a+2b}}),$ one has $D_2 = O_p(n^{-\frac{a+2b-1}{a+2b}}) = O_p(\delta_n^2).$ To sum up, we have $|r_i|^2 = O_p(\delta_n^2)$.

Now consider $B_i$. Due to

by Conditions C1, C2, C4, and the Hölder inequality, we obtain

In addition, noting that

together with the above inequality, one has

Recall that $\psi_{\tau_k}(u) = \rho'_{\tau_k}(u) = ({\tau_k}-1)I(u < -c^*)+\frac{(1-{\tau_k})}{c^*}uI(-c^*\leq u < 0)+\frac{\tau_k}{c^*}uI(0\leq u < c^*)+\tau_kI(u\geq c^*),$ and denote $\widetilde{Q}_{n}\big(\boldsymbol{V}_{n}, \boldsymbol{W}_{n}, \boldsymbol{S}_{n}\big) = \sum_{i = 1}^{n}\sum_{k = 1}^{K}w_{k}\rho_{\tau_k} \Big(Y_i-b_{0k}-\mathit{\boldsymbol{Z}}_i^{\rm T} \mathit{\boldsymbol{\alpha}}_0-\hat{\mathit{\boldsymbol{U}}}_i^{\rm T}\mathit{\boldsymbol{\gamma}}_0-\delta_n\big(\mathit{\boldsymbol{V}}_n^{\rm T} \tilde{\boldsymbol\eta}_{i}+\boldsymbol{W}_{n}^{\rm T}\hat{\mathit{\boldsymbol{A}}}_i+{\tilde B}_i+S_{n k}\big)\Big)$ $-\sum_{i = 1}^{n}\sum_{k = 1}^{K}w_{k}\rho_{\tau_k} \Big(Y_i-b_{0k}-\mathit{\boldsymbol{Z}}_i^{\rm T} \mathit{\boldsymbol{\alpha}}_0-\hat{\mathit{\boldsymbol{U}}}_i^{\rm T}\mathit{\boldsymbol{\gamma}}_0\Big)$. Then $\widetilde{Q}_{n}\big(\boldsymbol{V}_{n}, \boldsymbol{W}_{n}, \boldsymbol{S}_{n}\big)$ can be transformed into

where

Consider $D_1^*$. According to (A.2), we have

The proof of Theorem 3.1 in ^[6] indicates that $\left|\frac{1}{n}\sum\limits_{i = 1}^{n}\left[\left({\mathit{\boldsymbol{A}}}_{i}^{{\rm T}}h^m\right)^{2}\right]-1\right| = O_{P}\left(n^{-1 / 4} m^{1 / 2}(\log n)^{1 / 2}\right) = o_{P}(1)$ with $h^{m} = \mathit{\boldsymbol{W}}_n/\|\mathit{\boldsymbol{W}}_n\|$, which leads to

Observe that $\sum_{i = 1}^{n} h_{\tau_k}(b_{0k})\mathit{\boldsymbol{V}}_n^{\rm T}\tilde{\boldsymbol\eta}_{i} {\mathit{\boldsymbol{A}}}_i^{\rm T}\boldsymbol{W}_{n} = \mathit{\boldsymbol{V}}_n^{\rm T} P_{n}^{-1 / 2} \sum_{i = 1}^{n} h_{\tau_k}(b_{0k}) {\boldsymbol\eta}_{i} {{\mathit{\boldsymbol{U}}}}_{i}^{{\rm T}} \Lambda^{-1 / 2}\mathit{\boldsymbol{W}}_n,$ then by Conditions C1–C3, C5, $E \left[\sum_{i = 1}^{n} h_{\tau_k}(b_{0k}) {\eta}_{il} {{\mathit{\boldsymbol{U}}}}_{i}^{{\rm T}} \Lambda^{-1 / 2}\mathit{\boldsymbol{W}}_n\right] = 0$ and $E\left(\left[\sum_{i = 1}^{n} h_{\tau_k}(b_{0k}) {\eta}_{il} {{\mathit{\boldsymbol{U}}}}_{i}^{{\rm T}} \Lambda^{-1 / 2}\mathit{\boldsymbol{W}}_n\right]^{2}\right) = O(n m)$. Hence,

Similarly, $\sum_{i = 1}^{n} h_{\tau_k}(b_{0k})\mathit{\boldsymbol{V}}_n^{\rm T}\tilde{\boldsymbol\eta}_{i}S_{nk} = O_{p}\left(n^{1 / 2}\right), \; \sum_{i = 1}^{n} h_{\tau_k}(b_{0k}) {\mathit{\boldsymbol{A}}}_i^{\rm T}\boldsymbol{W}_{n}S_{nk} = O_{p}\left(n^{1 / 2} m^{1 / 2}\right).$ Then, together with formulas (A.3)–(A.5), we have

As for $D_2^*$, due to the continuity of $\psi_{\tau_{k}}(\cdot)$, we have

then

from which we get

For the term $D_3^*$, it is easy to show that

and

Combining with the equation (A.4), we can obtain

From the results about $D_1^*, D_2^*$ and $D^*_3$, it is easy to obtain that $\widetilde{Q}_{n}\left(\boldsymbol{V}_{n}, \boldsymbol{W}_{n}, \boldsymbol{S}_{n}\right)$ is dominated by the positive quadratic term $Cn\delta_n^2\left(\left\|{{\mathit{\boldsymbol{V}}}}_{n}\right\|^{2}+\left\|{{\mathit{\boldsymbol{W}}}}_{n}\right\|^{2}+\left\|{{\mathit{\boldsymbol{S}}}}_{n}\right\|^{2}\right)$. Hence, Eq (A.1) is established, and there exists local minimizer $\hat{\mathit{\boldsymbol{\gamma}}}$ such that

Now we consider the convergence rate of $\hat{\beta}$. Since

and based on the basic Condition C2, we have

Note that

Based on the Condition C2, Eq (A.7) and the orthogonality of $\{\hat{v}_j\}$, as well as $\|v_j-\hat{v}_j\|^2 = O_p(n^{-1}j^2)$, we can obtain

These lead to

Next, we turn to the asymptotic normality of $\hat{\alpha}$. Note that $Q_{n}(\boldsymbol{\alpha}, \boldsymbol{\gamma}, \boldsymbol{b})$ attains the minimal value at $(\hat{\boldsymbol{\alpha}}, \hat{\boldsymbol{\gamma}}, \hat{\boldsymbol{b}})$ with probability tending to 1 as $n$ tends to infinity. Then, we have the following score equations

Further, we can write (A.9) as $H_{n}+\sum_{k = 1}^{K}w_k B_{n 1}^{(k)}+\sum_{k = 1}^{K} w_k B_{n 2}^{(k)} = 0$, with

By simple calculations, we have $B_{n 1}^{(k)} = -\frac{1}{n} \sum_{i = 1}^{n} h_{\tau_k}(b_{0k})\left[{{\mathit{\boldsymbol{Z}}}}_{i}{{\mathit{\boldsymbol{Z}}}}_{i}^{\rm T}\left(\hat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}_{0}\right)+{{\mathit{\boldsymbol{Z}}}}_{i} \hat{\boldsymbol{U}}_{i}^{\rm T}\left(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{0}\right)\right](1+o_p(1)).$

Through calculating the mean and variance directly, we can obtain $B_{n 2}^{(k)} = o_{p}\left(\delta_{n}\right)$. Then, Eq (A.9) can be written as

Similarly, Eq (A.10) can be changed into

Now let $\Phi_{n} = \frac{1}{n} \sum_{i = 1}^{n} \hat{\mathit{\boldsymbol{U}}}_{i} \hat{\mathit{\boldsymbol{U}}}_{i}^{\rm T}, \; \Gamma_{n} = \frac{1}{n} \sum_{i = 1}^{n} \hat{{\mathit{\boldsymbol{U}}}}_{i} {{\mathit{\boldsymbol{Z}}}}_{i}^{\rm T}, \; \Gamma = E({\mathit{\boldsymbol{U}}}_{i} \mathit{\boldsymbol{Z}}_{i}^{\rm T}), \;$ $\Sigma_{n} = \frac{1}{n} \sum_{i = 1}^{n} \widetilde{\mathit{\boldsymbol{Z}}}_{i} \widetilde{\mathit{\boldsymbol{Z}}}_{i}^{\rm T}, \; \widetilde{\mathit{\boldsymbol{Z}}}_{i} = \mathit{\boldsymbol{Z}}_{i}-\Gamma_{n}^{\rm T} \Phi_{n}^{-1} \hat{\mathit{\boldsymbol{U}}}_{i}, \; \Upsilon_{n k} = \frac{1}{n} \sum_{i = 1}^{n} \hat{\boldsymbol{U}}_{i}\psi_{\tau_{k}}\left(e_{i}-b_{0 k}\right).$ Then, above equation allows $\sum_{k = 1}^{K}w_kh_{\tau_k}(b_{0k})(\hat{\mathit{\boldsymbol{\gamma}}}-\mathit{\boldsymbol{\gamma}}_{0}) = \sum_{k = 1}^{K}w_k\left(\Phi_{n}+o_{p}(1)\right)^{-1}[\Upsilon_{nk} +h_{\tau_k}(b_{0k})\Gamma_{n}(\mathit{\boldsymbol{\alpha}}_{0}-\hat{\mathit{\boldsymbol{\alpha}}})]$. Substituting it into (A.11), we obtain

Noting that

then, it is easy to conclude that

Note that

Then we have

It is easy to see that $\Phi_{n} = \Lambda+o_p(1), \; \Gamma_{n} = \Gamma+o_p(1)$. And based on Lemma 1 in ^[8] and the Condition $\mathrm{C}5,$ we can obtain $\frac{1}{n} \sum_{i = 1}^{n} \widetilde{\mathit{\boldsymbol{Z}}}_{i} \widetilde{\mathit{\boldsymbol{Z}}}_{i}^{\rm T} \stackrel{p}{\longrightarrow} \boldsymbol{\Sigma}\; (n \rightarrow \infty)$.

Through some simple calculations, one has

where $\boldsymbol{V} = \left(V_{k, l}\right)_{1 \leq k, l \leq K}$ with $V_{k l} = E[\psi_{\tau_{k}}\left(e_{i}-b_{0 k}\right)\psi_{\tau_{l}}\left(e_{i}-b_{0 l}\right)]$ and ${{\mathit{\boldsymbol{w}}}} = \left(w_{1}, \ldots, w_{k}\right)^{\rm T}$. Then, according to Eqs (A.12) and (A.13), the Slutsky's theorem, and the properties of multivariate normal distributions, we can obtain

Lastly, we prove the third conclusion of Theorem 2. By the definition of ${\rm{MSPE}}_{WCAHR}$, we have

Firstly, according to the previous proof process, we have $F_1 = \|{\Lambda}^{1/2}(\mathit{\boldsymbol{\hat{\gamma}}}-\mathit{\boldsymbol{\gamma}}_0)\|^2 = O_p(\frac{m}{n})$. As for $F_2$, based on triangle inequality and $C_R$ inequality,

By Eq (A.8), $F_{21} = m\sum\limits_{j = 1}^{m}{\gamma}_{0 j}^2\|v_j-\hat{v}_j\|^2 = O_p\left(\frac{m}{n}\right)$ by. As for $F_{22}$, it is easy to know ${\sum\limits_{j = 1}^{m}\lambda_j^{-1}(v_j-\hat{v}_j)^2 = o_p(1)}$, then $F_{22} = \|{\Lambda}^\frac{1}{2}(\mathit{\boldsymbol{\hat{\gamma}}}-\mathit{\boldsymbol{\gamma}}_0)\|^2\sum\limits_{j = 1}^{m}\lambda_j^{-1}(v_j-\hat{v}_j)^2 = o_p\left(\frac{m}{n}\right).$ Next, by (A.8), we have $F_3 = \sum\limits_{j = m+1}^{\infty}{\gamma}_{0 j}^2\lambda_j = O(m^{-a-2b+1})$. By (A.6), we know $F_4 = O_p\left(\frac{m}{n}\right)$ and $F_5 = O_p\left(\frac{m}{n}\right)$. Then, we have

The proof of Theorem 2 is complete.