# Least Absolute Deviation Estimate for Functional Coefficient Partially Linear Regression Models

Authors: Yanqin Feng; Guoxin Zuo; Li Liu
Journal: Journal of Probability and Statistics
 (2012)

---

## Abstract

The functional coefficient partially linear regression model is a useful generalization of the nonparametric
model, partial linear model, and varying coefficient model. In this paper, the local linear technique and theL1 method are employed to estimate all the functions in the functional coefficient partially linear regression model. The asymptotic properties of the proposed estimators are studied. Simulation studies are conducted to show the validity of the estimate procedure.


---

## Body

## 1. Introduction

In this paper, we are concerned with a functional coefficient partially linear regression model (FCPLR), that is,(1.1)Y=a0(X)+∑j=1p‍aj(U)Zj+ε,
where X and U are random explanatory variables, Z=(Z1,…,Zp)′ is a random vector, and aj(·) is some measurable function from R to R for j=0,…,p. We call a0(·) the intercept function, and {aj(·)}, j=1,…,p, the coefficient functions. As usual, ε denotes the errors with zero-mean and fixed variance.The FCPLR model, first introduced by Wong et al. [1], is a generalization of the nonparametric model, partial linear model, and varying coefficient model. Zhu et al. [2] studied a similar functional coefficient model, a functional mixed model, using a new Bayesian method. Model (1.1) reduces to a varying coefficient regression mode if the intercept function a0(·) is a constant function and a partially linear regression model when the coefficient functions {aj(·)}j=1p are constants. Many researchers, for example, Aneiros-Pérez and Vieu [3], made contributions to studying this kind of model. When X and U in model (1.1) coincide, the model becomes a semiparametric varying coefficient model discussed by Ahmad et al. [4]. Since the FCPLR model combines the nonparametric and functional coefficient regression model, its flexibility makes it attractable in various regression problems [1].Statistical inference for the FCPLR model mainly includes the estimations of the intercept functiona0(·) and the coefficient functions {aj(·)}j=1p. To estimate the unknown functions in the nonparametric/semiparametric regression models, many statistical inference methods have been proposed over the past decades, such as the kernel estimate method [5–7], spline smoothing [8, 9], and two-step estimation method [10, 11]. Wong et al. [1] employed local linear regression method and integrated method to give the initial estimates of all functions in the FCPLR model. All the papers mentioned above used the least-squares technique to obtain the estimators of the unknown coefficient functions. The least-squares estimators (L2 method), of course, have some good properties, especially for the normal random errors case. It is well known that, however, the least-squares method will perform poor when the random errors have a heavy-tailed distribution in that it is highly sensitive to extreme values and outliers. This motivates us to find more robust estimation methods instead of the aforementioned inference methods for model (1.1).Local linear approximation is a good method for nonparametric regression problems [12], and the L1 method based on the least absolute deviations overcomes the sensitivity caused by outliers. As noted in Wang and Scott [13] and Fan and Gijbels [12], among many robust estimation methods, the L1 method based on the local least absolute deviations behaves quite well. In this paper, we adopt the L1 method, accompany with the local linear technique and the integrated method to estimate all the unknown functions in model (1.1). Furthermore, the estimating problem can be reduced to a linear programming problem, and the numerical solutions are obtained quickly by some available softwares subsequently (e.g., Matlab is very useful for this kind of problems). The main difficulty encountered in the proof of the asymptotic normalities is that the L1 estimates have no closed form. This paper shows the asymptotic normalities of L1 estimators through a method completely different from those based on the L2 method, and the simulation results show that the L1 method is a robust method indeed.The rest of this paper is organized as follows. In Section2, we describe the estimation method and the associated bandwidth selection procedure. Section 3 gives the the asymptotic theories of the estimators. Simulation studies are conducted in Section 4. A real application is given in Section 5. Section 6 gives the proofs of the main results.


## 2. Least Absolute Deviation Estimate

This section gives the main idea of the proposed estimation method; that is, local linear polynomials are used to approximate the nonparametric function and the functional coefficients, and the least absolute deviation technique is used to find the best approximation. Bandwidth selection technique is also discussed in this section. Throughout this paper, we suppose that{(Yi,Xi,Ui,Zi')′}i=1n is an i.i.d. sample from model (1.1) and assume the following conditions.Assumptions
(1)
Ω
=
E
(
Z
Z
'
∣
X
=
x
0
,
U
=
u
0
) is a positive definite matrix, E(Z∣X=x0,U=u0)=ω.
(2)
Bandwidthh subjects to h~n-1/6.
(3)
Random errorε, with zero mean and zero median, is independent of Z conditional on (X,U). The conditional probability density g(·|x,u) of ε given random variables X and U is continuous in a neighborhood of the point 0, and g(0|x,u)>0. ε1,…,εn are independent and identically distributed.
(4)
The density functionsf(x,u), f1(x), f2(u) of (X,U), X, and U are continuous in neighborhoods of (x0,u0), x0 and u0, and f(x0,u0)>0.
(5)
All functionsaj(u), j=1,…,p and a0(x) are twice continuously differentiable in neighborhoods of u0 and x0, respectively.
(6)
Kernel functionsKk(·), k=1,2 are bounded, nonnegative, and compactly supported.
(7)
max
1
≤
i
≤
n
∥
Z
i
∥
=
o
p
(
n
1
/
3
).
To simplify typesetting, we introduce the following symbols:(2.1)μi(k)=∫‍viKk(v)dv,γi(k)=∫‍viKk2(v)dvfork=1,2.i=1,2,….


### 2.1. Local Linear Estimate Based on Least Absolute Deviation

The main idea is to approximate the functional coefficientsaj(·) by linear functions for j=0,…,p, that is, a0(x) can be approximated by
(2.2)a0(x)≈a0(x0)+a0′(x0)(x-x0)
for x in a neighborhood of x0 within the closed support of X, and aj(u) by
(2.3)aj(u)≈aj(u0)+aj′(u0)(u-u0),j=1,…,p
for u in a neighborhood of u0 within the closed support of U. For simplicity, denoting a0(x0), a0′(x0) as a0, b0, and aj(u0), aj′(u0) as aj, bj for j=1,…,p. The local linear least absolute deviation estimate (L1 estimate) of the unknown parameters a0, b0, a:=(a1,…,ap)′ and b:=(b1,…,bp)′, denoted, respectively, by a^0, b^0, a^, b^, is the optimal solution of the minimization problem as follows:
(2.4)min∑i=1n‍|Yi-a0-b0(Xi-x0)-a′Zi-b′(Ui-u0)Zi|K1,h(Xi-x0)K2,h(Ui-u0),
where Kk,h(·)=Kk(·/h), k=1,2 are the given kernel functions and h is the chosen bandwidth. The optimization problem is equivalent to the following linear programming problem:
(2.5)min∑i=1n‍(ei++ei-)K1,h(Xi-x0)K2,h(Ui-u0)s.t.a0+b0(Xi-x0)+a′Zi+b′(Ui-u0)Zi+ei+-ei-=Yi,t.Zi+ei+-ei-ei-=Yiei+≥0,ei-≥0,i=1,…,n.
There are many algorithms available for the optimal solution of problem (2.5); for example, the feasible direction method can be directly used to compute the optimal solution [14], and the numerical solution of (2.5) can be quickly computed by a series of Matlab functions.By the integrated method [1], the estimator of the intercept function a0(x) is defined by
(2.6)a~0(x0)=1n∑k=1n‍a^0(x0,Uk),
and the estimators of the coefficient functions aj(u0), j=1,…,p are defined by
(2.7)a~j(u0)=1n∑k=1n‍a^j(Xk,u0),j=1,…,p.
We focus our main task in establishing the asymptotic distributions of the estimators a~0(x0) and a~j(u0) for j=1,…,p.


### 2.2. Selection of Bandwidth

It is well known that the choice of the bandwidth strongly influences the adequacy of the estimators. We use an automatic bandwidth choice procedure in this paper; that is, the absolute cross-validation (ACV) method, and the ACV bandwidthhACV is defined as
(2.8)hACV=argminACV(h)=1n∑k=1n‍|Yk-a^0-k(Xk)-∑j=1p‍a^j-k(Uk)Zkj|,
where a0-k(Xk) and aj-k(Uk), j=1,…,p are constructed based on observations with size n-1 by leaving out the kth observation ((Yk,Xk,Uk)′,Zk′)′. According to Wang and Scott [13], this bandwidth is better than cross-validation (CV) bandwidth hCV. The latter one is suggested by Rice and Silverman [15] and often used in curve regression as in Hoover et al. [6] and Wu et al. [7].


## 2.1. Local Linear Estimate Based on Least Absolute Deviation

The main idea is to approximate the functional coefficientsaj(·) by linear functions for j=0,…,p, that is, a0(x) can be approximated by
(2.2)a0(x)≈a0(x0)+a0′(x0)(x-x0)
for x in a neighborhood of x0 within the closed support of X, and aj(u) by
(2.3)aj(u)≈aj(u0)+aj′(u0)(u-u0),j=1,…,p
for u in a neighborhood of u0 within the closed support of U. For simplicity, denoting a0(x0), a0′(x0) as a0, b0, and aj(u0), aj′(u0) as aj, bj for j=1,…,p. The local linear least absolute deviation estimate (L1 estimate) of the unknown parameters a0, b0, a:=(a1,…,ap)′ and b:=(b1,…,bp)′, denoted, respectively, by a^0, b^0, a^, b^, is the optimal solution of the minimization problem as follows:
(2.4)min∑i=1n‍|Yi-a0-b0(Xi-x0)-a′Zi-b′(Ui-u0)Zi|K1,h(Xi-x0)K2,h(Ui-u0),
where Kk,h(·)=Kk(·/h), k=1,2 are the given kernel functions and h is the chosen bandwidth. The optimization problem is equivalent to the following linear programming problem:
(2.5)min∑i=1n‍(ei++ei-)K1,h(Xi-x0)K2,h(Ui-u0)s.t.a0+b0(Xi-x0)+a′Zi+b′(Ui-u0)Zi+ei+-ei-=Yi,t.Zi+ei+-ei-ei-=Yiei+≥0,ei-≥0,i=1,…,n.
There are many algorithms available for the optimal solution of problem (2.5); for example, the feasible direction method can be directly used to compute the optimal solution [14], and the numerical solution of (2.5) can be quickly computed by a series of Matlab functions.By the integrated method [1], the estimator of the intercept function a0(x) is defined by
(2.6)a~0(x0)=1n∑k=1n‍a^0(x0,Uk),
and the estimators of the coefficient functions aj(u0), j=1,…,p are defined by
(2.7)a~j(u0)=1n∑k=1n‍a^j(Xk,u0),j=1,…,p.
We focus our main task in establishing the asymptotic distributions of the estimators a~0(x0) and a~j(u0) for j=1,…,p.


## 2.2. Selection of Bandwidth

It is well known that the choice of the bandwidth strongly influences the adequacy of the estimators. We use an automatic bandwidth choice procedure in this paper; that is, the absolute cross-validation (ACV) method, and the ACV bandwidthhACV is defined as
(2.8)hACV=argminACV(h)=1n∑k=1n‍|Yk-a^0-k(Xk)-∑j=1p‍a^j-k(Uk)Zkj|,
where a0-k(Xk) and aj-k(Uk), j=1,…,p are constructed based on observations with size n-1 by leaving out the kth observation ((Yk,Xk,Uk)′,Zk′)′. According to Wang and Scott [13], this bandwidth is better than cross-validation (CV) bandwidth hCV. The latter one is suggested by Rice and Silverman [15] and often used in curve regression as in Hoover et al. [6] and Wu et al. [7].


## 3. Asymptotic Theory

This section gives the asymptotic distribution theories of the estimators. Using Taylor’s expansion, for|Xi-x0|≤Ch and |Ui-u0|≤Ch, we have
(3.1)a0(Xi)=a0(x0)+a0′(x0)(Xi-x0)+12a0′′(ξi)(Xi-x0)2,aj(Ui)=aj(u0)+aj′(u0)(Ui-u0)+12aj′′(ηij)(Ui-u0)2,forj=1,…,p,
where ξi is between x0 and Xi, and ηij is between Ui and u0. Let a(u0)=(a1(u0),…,ap(u0))′, b(u0)=(a1′(u0),…,ap′(u0))′, c(u0)=(a1′′(u0),…,ap′′(u0))′, then we have
(3.2)(a^0,b^0,a^′,b^′)′=argmin∑i=1n‍|Yi-a0-b0(Xi-x0)-a′Zi-b′(Ui-u0)Zi|K1,h(Xi-x0)K2,h(Ui-u0)=argmin∑i=1n‍{{|12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi|}|[12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi]nh2[(a0-a0(x0))+(b0-a0′(x0))(Xi-x0)]1nh2=argmin+∑i=1n‍+nh2[(a-a(u0))′Zi+(b-b(u0))′(Ui-u0)Zi]1nh2=argmin+∑i=1n‍-[12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi]|=argmin∑i=1n‍-|12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi|}×K1,h(Xi-x0)K2,h(Ui-u0).The aim of this paper is to study the asymptotic behavior ofnh(a~0-a0(x0)) and nh(a~-a(u0)). Combining some technique reasons, we first introduce the new variables α0=nh2(a0-a0(x0)), α=nh2(a-a(u0)), β0=nh2[(b0-a0′(x0))h], and β=nh2[(b-b(u0))h] [16] and form a new equivalent problem as follows:
(3.3)(α^0,β^0,α^′,β^′)′=argmin∑i=1n‍{{[12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi]}|[12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi][α0+β0(Xi-x0h)]1nh2+[α′Zi+β′(Ui-u0h)Zi]1nh2argm2-in-[12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi]|argmin-|12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi|}×K1,h(Xi-x0)K2,h(Ui-u0).Let(3.4)Fn=E(Sn),Rn=Sn-Fn+∑i=1n‍{[α0+β0(Xi-x0h)]1nh2Sn-FnFn+∑i=1n5‍+[α′Zi+β′(Ui-u0h)Zi]1nh2}K1,h(Xi-x0)K2,h(Ui-u0)sgn(εi):=Sn-Fn+∑i=1n‍Lni,
where Sn is the objective function of the equality above and sgn(·) is the sign function.Since theL1 estimators have no closed forms, we first give the limit form of the function Fn, which is critical to obtain the asymptotic properties of the estimators.Theorem 3.1.
Suppose Assumptions (1)–(7) hold, andn→∞, then for any fixed α0, β0, α, β,Fn converges to F(α0,β0,α,β), which is defined as
(3.5)g(0∣x0,u0)f(x0,u0){(α02+2α0β0μ1(1)+β02μ2(1))α′α+2μ1(2)+1+(α′Ωα+2μ1(2)α′Ωβ+μ2(2)β′Ωβ)+2α0(α′ω+μ1(2)β′ω)α′α+2μ1(2)+1+2β0μ1(1)(α′ω+μ1(2)β′ω)α′α+2μ1(2)+1-a′′(x0)(α0μ2(1)+β0μ3(1)+μ2(1)α′ω+μ2(1)μ1(2)β′ω)α′α+2μ1(2)+1-c(u0)′(α0μ2(2)ω+β0μ1(1)μ2(2)ω+μ2(2)Ωα+μ3(2)Ωβ)(α02+2α0β0μ1(1)+β02μ2(1))}.Remark 3.2.
If the kernel functionsK1(·), K2(·) are symmetric about zero and Lipschitz continuous, the limit form of F(α0,β0,α,β) can be simplified as
(3.6)F(α0,β0,α,β)=g(0∣x0,u0)f(x0,u0){α02+μ2(1)β02+α′Ωα+μ2(2)β′Ωβ=g(0∣x0,u0)f(x0,u0)1+2α0α′ω-μ2(1)a′′(x0)(α0+α′ω)-μ2(2)c(u0)′(α0ω+Ωα)(1)β02}.
Now we are in the position to state the asymptotic properties of the estimators.Theorem 3.3.
Suppose Assumptions (1)–(7) hold andn→∞, then one has
(3.7)nh(a~0(x0)-a0(x0)-h2[μ22(1)-μ3(1)μ1(1)]a′′(x0)2[μ2(1)-μ12(1)])→dN(0,σα02),(3.8)nh(a~(u0)-a(u0)-h2[μ22(2)-μ1(2)μ3(2)]c(u0)2[μ2(2)-μ12(2)])→dNp(0,Σα),
where
(3.9)σα02={(1-ω′Ω-1ω)γ0(1)r22(x0)4+μ12(1)[μ12(1)γ0(1)-2μ1(1)γ1(1)+γ2(1)]r22(x0)4g2(0∣x0,u0)f2(x0,u0)[μ12(1)-μ2(1)]2σα02+=1-(1-ω′Ω-1ω)μ1(1)[μ1(1)γ0(1)-γ1(1)]r22(x0)2g(0∣x0,u0)f(x0,u0)[μ12(1)-μ2(1)](1-ω′Ω-1ω)γ0(1)r22(x0)4},Σα={(Ω-ωω′)-1γ0(2)r12(u0)4+μ12(2)[μ12(2)γ0(2)-2μ1(2)γ1(2)+γ2(2)]r12(u0)Ω-14g2(0∣x0,u0)f2(x0,u0)[μ12(2)-μ2(2)]2∑α+5=-μ1(2)[μ1(2)γ1(2)-γ1(2)]r12(u0)2g(0∣x0,u0)f(x0,u0)[μ12(2)-μ2(2)](Ω-ωω′)-1γ0(2)r12(u0)4}
with r12(u0)=∫‍f12(u)f(u,u0)du and r22(x0)=∫‍f22(u)f(x0,u)du.Remark 3.4.
If the kernel functionsK1(·),K2(·) are symmetric about zero and Lipschitz continuous, the results of Theorem 3.3 can be simplified as
(3.10)nh(a~0(x0)-a0(x0)-h2μ2(1)a′′(x0)2)→dN(0,(1-ω′Ω-1ω)γ0(1)r22(x0)4),nh(a~(u0)-a(u0)-h2μ2(2)c(u0)2)→dNp(0,(Ω-ωω′)-1γ0(2)r12(u0)4).Remark 3.5.
Here, we have considered estimation method and asymptotic distributions for the case that two bandwidths are same. It is important to note that similar asymptotic theories can be obtained for the case that two bandwidths with the same order are different.Remark 3.6.
If we consider different bandwidthh1,h2 for kernel functions K1 and K2. Furthermore, suppose the Assumptions (2) and (7) are replaced by h15h2=o(1/n) and h1h25=o(1/n), for example, h1~n-1/6,h2~n-1/5, and max∥Zi∥=op(min{h2-2,nh1h2}), respectively. Similar results also will be obtained except that all the second-order derivatives in the results will disappear.Remark 3.7.
This paper restricts the study to one-dimensional variableX. The ideas used here can be adapted to higher dimensional variable X, for example, consider d-dimensional variable X for the case of same bandwidths. Similar asymptotic distribution results for nhd(a~(x0)-a(x0)) and nh(a~(u0)-a(u0)) can be obtained under the assumptions with Assumptions (2) and (7) being replaced by h~n-1/(5+d) and max∥Zi∥=op(n2/(5+d)), respectively.


## 4. Simulations

In this section, we carry out some simulations to illustrate the performance ofL1-method, and compare the performance of our L1-method with that of the L2-method. All the following simulations are conducted for sample size n=100.The following example is considered:(4.1)Y=a0(X)+a1(U)Z1+a2(U)Z2+ε,
where a0(x)=x+4exp(-8x2), a1(u)=2sin(1.5πu+u), a2(u)=4u(1-u), (X,U)′~U([-0.5,0.5]2), Z1, Z2 are normally distributed with correlation coefficient 1/2, the marginal distributions of Z1 and Z2 are standard normal, ε~N(0,0.22), and ε,(X,U) and Z1, Z2 are mutually independent.In each simulation, theL1 estimators of a0(x), a1(u), a2(u) were computed by solving the minimization problem (2.5) and using the integrated method described in (2.6) and (2.7). We use the Epanechnikov kernel, K(u)=(3/4)(1-u2)I|u|<1, for every Kl(·), l=1,2. All bandwidths in a model are selected by the method proposed in Section 2.To evaluate the asymptotic results given in Theorem3.3, the quantile-quantile plots of the estimators are constructed. Figure 1 presents the quantile-quantile plots for a~0(0),a~1(0),a~2(0) with sample size n=100 and 100 replications, respectively, and these plots reveal that the asymptotic approximation is reasonable.Normal QQ-plot fora~0(0),a~1(0) and a~2(0)(n=100).
(a)
(b)
(c)Figure2 displays the true function curves of a0(x),a1(u), and a2(u) and their estimated curves with sample size n=100 and one replication. We can see from the figure that the L1 estimates perform well.Simulation results, for example. Dotted curves areL1 estimators; solid curve are the true curves.
(a)
(b)
(c)In order to illustrate that theL1 method is a robust method. Figure 3 displays the estimated curves with four outliers, that is, (-0.4538,-0.4402,-10), (0.13377,0.1393,20), (0.2703,0.3085,20), (0.4174,0.4398,18) for the three element array (x,u,y). From Figure 3, we can see that the L1 estimate also has a good performance even in the presence of four large singular points of Y. The fact that outliers have little influence on the L1 estimates is displayed in Figure 3, so it is a robust method.Simulation results, for example, with outliers. Dotted curves areL1 estimators; solid curve are the true curves.
(a)
(b)
(c)By solving the following minimization problem:(4.2)min∑i=1n‍(Yi-a0-b0(Xi-x0)-a′Zi-b′(Ui-u0)Zi)2K1,h(Xi-x0)K2,h(Ui-u0),
we can obtain similarly the L2 estimators of the functions a0(x),a1(u),a2(u) by the equations (2.6) and (2.7). For comparing the L1-method with the L2 method. We simulated the function a0(x) by L2 method and display the fitted curves with (without) outliers for sample sizes n=100 and 1000 replications in Figure 4. We can see that L2 estimate cannot perform well for the data sparsity and singularity, and three small outlying data points make the estimated curve by the L2 method deviate from the true curve significantly. Combining Figures 2, 3, and 4, we conclude that the L1 method performs better than the L2 method, the L1 method is a robust method.Simulation results fora^0(x). Dotted curves are L2 estimators; solid curve are the true curves. (N) plot without outliers, (Y) plot with three outliers (x,y)=(-0.4538,-5), (0.1948,8), and (0.4649,5).
(a)
(b)Finally, for further comparing theL1-estimate with the L2-estimate method, we also assess their performance via the weighted average squared error (WASE), which is defined as
(4.3)WASE=1n∑i=1n‍((a0(xi)-a^0(xi))2(range(a0))2+(a1(ui)-a^1(ui))2(range(a1))2+(a2(ui)-a^2(ui))2(range(a2))2),
where range (ak) for k=0,1,2 are the ranges of the functions a0(x), a1(u), and a2(u). The weights are introduced to account for the different scales of the functions. We conducted 200 replications with sample size n=100. For the bandwidths and the Epanechnikov kernels used in the simulations, we obtain the mean and standard deviation of the WASE are 0.1201 and 0.0183 for the L1 method, and 1.6613 and 0.8936 for the L2 method. We can see that L1 method outperforms L2 method.


## 5. Application

A real data is analyzed by the proposedL1-method in this section. The classic gas furnace data was studied recently by Wong et al. [1]. The data set includes 296 samples (It,Ot), (t=1,…,296) measured at a fixed interval of 9 seconds, where It’s represent the input gas rate in cubic feet per minute, and Ot’s represent the concentration of carbon dioxide in the gas out of the furnace. Similar to the procedures of Wong et al. [1], the original data are transformed as xt=(It+2.716)/5.55, yt=(Ot-45.6)/14.9 for t=1,…,296 such that both x’s and y’s are limited in the interval [0,1], and the model
(5.1)yt=a0(xt-3)+∑j=14‍aj(yt-3)yt-j+ϵt
is used to fit the data. The first 250 samples are used to establish the model, and the remained 46 samples are used for prediction.In the proposed method, Epanechnikov kernel is used and all the bandwidthshj(j=0,…,4) are selected as 0.14 via cross validation for simplicity. The mean absolute error (MAE) and the mean squared error (MSE) for fitting and forecasting are listed as follows.Fitting:MAE=0.3189, MSE=0.1584; forecasting: MAE=0.3705, MSE=0.2127.Since the model is chosen based on theL2 errors, the results are not perfect as that in Wong et al. [1]. Moreover, the data set does  not contain obvious outliers, so the advantage of the L1 estimation method is not apparent. Compared to the results showed in Wong et al. [1], the difference between the fitting MAE/MSE and the forecasting MAE/MSE is small. The reason is that L1-methods are employed in our method. The fitted values {O^t}t=5250 and predictive values {O^t}t=251296 are shown in Figure 5. These results indicate that the estimated results are reasonable.(a) Fitted values ofOt(t=5,…,250) for gas data; (b) Predictive values Ot(t=251,…,296) for gas data. Dotted lines are fitted or predicted; solid lines are observed.
(a)
(b)


## 6. Proofs of the Main Results

Before completing the proofs of the main results, we give the following useful lemma first.Lemma 6.1.
Suppose Assumptions (1)–(7) hold, then for any fixedα0,β0,α,β,Rn converges to 0 in probability, that is, Rn=op(1) as n→∞.Proof of Lemma6.1.
Let(6.1)dn=max1≤i≤n{|a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij||[α0+β0(Xi-x0h)]1nh2+[α′Zi+β′(Ui-u0h)Zi]1nh2|dn5max1≤i≤n35+12|a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij|}.
Note that |Xi-x0|≤Ch, |Ui-u0|≤Ch and Assumptions (2), (5), and (7), for any fixed α0, β0, α and β, we have that
(6.2)dn=op(1)asn→∞.
Let(6.3)Tni=Lni+{-[12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi]|{∑j=1p‍aj′′(ηij)(Ui-u0)2Zij}[α0+β0(Xi-x0h)]1nh2+[α′Zi+β′(Ui-u0h)Zi]1nh2Tni=+1-[12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi]|Tni=+1-|12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi|[12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij)+εi]}×K1,h(Xi-x0)K2,h(Ui-u0).
It can be easily seen that Tni=0 if |εi|≥dn. Hence we have
(6.4)|Tni|≤2|{[α0+β0(Xi-x0h)]+[α′Zi+β′(Ui-u0h)Zi]}|1nh2×K1,h(Xi-x0)K2,h(Ui-u0)I{|εi|<dn},
where I(·) is the indicator function. For any ϵ>0, δ>0, we have
(6.5)P{|Rn|>ϵ}≤P{dn≥δ}+P{I{dn<δ}|Rn|>ϵ}.
Since E(sgn(εi))=0, we have
(6.6)E(I{dn<δ}Rn2)=E{I{dn<δ}(Sn-Fn+Ln)2}=E{I{dn<δ}(∑i=1n‍(Tni-E(Tni)))2},
here Ln=∑i=1n‍Lni. Combining (6.2) and (6.4), we have
(6.7)limn→∞E(I{dn<δ}Rn2)=limn→∞E{∑‍i≠jE((Tni-E(Tni))(Tnj-E(Tnj))∣X,U,Z)]I{dn<δ}[∑i=1n‍E((Tni-E(Tni))2∣X,U,Z)limn→∞I{dn<δ}555E+∑‍i≠jE((Tni-E(Tni))(Tnj-E(Tnj))∣X,U,Z)∑i=1n‍E((Tni-E(Tni))2∣X,U,Z)]}≤8limn→∞E{[(α0+β0(Xi-x0h))2+(α′Zi+β′(Ui-u0h)Zi)2][(α0+β0(Xi-x0h))2+(α′Zi+β′(Ui-u0h)Zi)2]8limn→∞55E×1h2K1,h2(Xi-x0)K2,h2(Ui-u0)P{|εi|≤δ∣X,U}[(α0+β0(Xi-x0h))2+(α′Zi+β′(Ui-u0h)Zi)2]}≤16δg(0∣x0,u0)f(x0,u0)[γ0(2)(α02γ0(1)+2α0β0γ1(1)+β02γ2(1))16δg(0∣x0,u0)ef(x0,u0)+γ0(1)(γ0(2)α′Ωα16δg(0∣x0,u0)ef(x0,u0)+11γ0(1)+2γ2(1)α′Ωβ+γ2(2)β′Ωβ)+o(1)(α02γ0(1)+2α0β0γ1(1)+β02γ2(1))]→0,
when δ→0. Combining (6.2), (6.5), and the argument E(I{dn<δ}Rn2)→0 as δ→0 and n→∞, the desired conclusion follows by using the Chebyshev’s inequality. This completes the proof of Lemma 6.1.Proof of Theorem3.1.
Set(6.8)A(Xi,Ui,Zi)=(α0+β0(Xi-x0h))1nh2+(α′Zi+β′(Ui-u0h)Zi)1nh2,B(Xi,Ui,Zi)=12(a0′′(ξi)(Xi-x0)2+∑j=1p‍aj′′(ηij)(Ui-u0)2Zij).
By Lemma 6.1, under Assumptions (1)–(7), A(Xi,Ui,Zi) and B(Xi,Ui,Zi) converge to zero as n→∞.
Start from the equality,(6.9)Fn=I{dn<δ}Fn+I{dn≥δ}Fn.
We first give the limit form of Fn, for fixed α0, β0, α, β, we have
(6.10)I{dn<δ}Fn=I{dn<δ}E{∑i=1n‍|A(Xi,Ui,Zi)-B(Xi,Ui,Zi)-εi|-|B(Xi,Ui,Zi)+εi|}×K1,h(Xi-x0)K2,h(Ui-u0)=I{dn<δ}nE{|A(X1,U1,Z1)-B(X1,U1,Z1)-ε1|-|B(X1,U1,Z1)+ε1|}×K1,h(X1-x0)K2,h(U1-u0)=I{dn<δ}nE{E{[|A(X1,U1,Z1)-B(X1,U1,Z1)-ε1|-|B(X1,U1,Z1)+ε1|]I{dn<δ}n{dn.}E×K1,h(X1-x0)K2,h(U1-u0)∣X,U,Z}}=I{dn<δ}nE{∫δ>w>A(X1,U1,Z1)-B(X1,U1,Z1)[w+B(X1,U1,Z1)-A(X1,U1,Z1)]g(w∣X1,U1)dwI{dn<δ}55E+∫-δ<w<A(X1,U1,Z1)-B(X1,U1,Z1)[A(X1,U1,Z1)-B(X1,U1,Z1)-w]g(w∣X1,U1)dwI{dn<δ}5+E-∫δ>w>-B(X1,U1,Z1)[w+B(X1,U1,Z1)]g(w∣X1,U1)dwI{dn<δ}5+E+∫-δ<w<-B(X1,U1,Z1)[B(X1,U1,Z1)+w]g(w∣X1,U1)dwI{dn<δ}5+E+∫w>δ-g(w∣X1,U1)A(X1,U1,Z1)dwI{dn<δ}5+E+∫w<-δg(w∣X1,U1)A(X1,U1,Z1)dw}K1,h(X1-x0)K2,h(U1-u0).
By the Integral Mean Value Theorem (refer to the Appendix), we have
(6.11)I{dn<δ}Fn=nI{dn<δ}E{[δ22+(A(X1,U1,Z1)-B(X1,U1,Z1))δ+(B(X1,U1,Z1)-A(X1,U1,Z1))22]g(τ1∣X1,U1)nI{dn<δ}E+E×[δ22+(B(X1,U1,Z1)-A(X1,U1,Z1))δ+(B(X1,U1,Z1)-A(X1,U1,Z1))22]nI{dn<δ}+EE+g(τ2∣X1,U1)nI{dn<δ}+EE×[δ22+(A(X1,U1,Z1)-B(X1,U1,Z1))δ+(B(X1,U1,Z1)-A(X1,U1,Z1))22]nI{dn<δ}+EE+A(X1,U1,Z1)[G(δ∣X1,U1)+G(-δ∣X1,U1)-1]nI{dn<δ}+EE-g(τ3∣X1,U1)[δ22+B(X1,U1,Z1)δ+B2(X1,U1,Z1)2]nI{dn<δ}+EE-g(τ4∣X1,U1)[δ22-B(X1,U1,Z1)δ+B2(X1,U1,Z1)2][δ22+(A(X1,U1,Z1)-B(X1,U1,Z1))δ+(B(X1,U1,Z1)-A(X1,U1,Z1))22]}×K1,h(X1-x0)K2,h(U1-u0),
where G(·|·) is the conditional probability distribution function of ε, and τ1,τ2,τ3,τ4 converge to zero as δ→0 and n→∞. Then by Assumptions (1)–(5), for any small enough δ>0, we have
(6.12)I{dn<δ}Fn=I{dn<δ}nE{(g(0∣X1,U1)+o(1))(A2(X1,U1,Z1)-2A(X1,U1,Z1)B(X1,U1,Z1))}×K1,h(X1-x0)K2,h(U1-u0)=I{dn<δ}(μ2(2)Ωα+μ3(2)Ωβ)(α02+2α0β0μ1(1)+β02μ2(1))}g(0∣x0,u0)f(x0,u0)α02+2α0β0×{(α02+2α0β0μ1(1)+β02μ2(1))×α02+2α0β0(1)+(α′Ωα+2μ1(2)α′Ωβ+μ2(2)β′Ωβ)×α02+2α0β0(1)+2α0(α′ω+μ1(2)β′ω)+2β0μ1(1)(α′ω+μ1(2)β′ω)×α02+2α0β0(1)-a′′(x0)(α0μ2(1)+β0μ3(1)+μ2(1)α′ω+μ2(1)μ1(2)β′ω)×α02+2α0β0(1)-c(u0)′(α0μ2(2)ω+β0μ1(1)μ2(2)ω×α02-c(u0)′1+2α0β0(1)+μ2(2)Ωα+μ3(2)Ωβ)(α02+2α0β0μ1(1)+β02μ2(1))}+o(1))=I{dn<δ}(F(α0,β0,α,β)+o(1)),Fn=F(α0,β0,α,β)+o(1)-[F(α0,β0,α,β)+o(1)-Fn]I{dn≥δ}.
Note that, for any fixed α0, β0, α, β,dn→0 as n→∞, we draw the desired conclusion. This completes the proof of Theorem 3.1.Proof of Theorem3.3.
Since the proofs of (3.7) and (3.8) are quite similar, we only give the proof of (3.8).
Start from the equality(6.13)Ln=∑i=1n‍{[α0+β0(Xi-x0h)]+[α′Zi+β′(Ui-u0h)Zi]}1nh2×K1,h(Xi-x0)K2,h(Ui-u0)sgn(εi).
By Lemma 6.1 and Theorem 3.1, for fixed α0, β0, α, β, we have
(6.14)Sn=F(α0,β0,α,β)-Ln+Rn+o(1)≜F(α0,β0,α,β)-Ln+R-n,
hence
(6.15)Sn+Ln=F(α0,β0,α,β)+R-n.
Note that
(6.16)ELn2=∑i=1n‍ELni2≤2nE{(α0+β0(Xi-x0h))21nh2K1,h2(Xi-x0)K2,h2(Ui-u0)2n5E+(α′Zi+β′(Ui-u0h)Zi)21nh2K1,h2(Xi-x0)K2,h2(Ui-u0)}=2[γ0(2)(α02γ0(1)+2α0β0γ1(1)+2β02γ2(1))+γ0(1)(α′Ωα+2γ1(1)α′Ωβ+γ2(2)β′Ωβ)].
We obtain that Ln is bounded in probability for any fixed α0, β0, α, β. Thus, for any fixed α0,β0,α,β, the random convex function Sn+Ln converges to F(α0,β0,α,β). According to the convexity lemma [17], we can deduce that, for any compact set 𝕂,
(6.17)sup(α0,β0,α′,β′)′∈𝕂|R-n|=op(1),
when δ→0. By the proof of Theorem 2 in Wang [18], we obtain that the “limit" here is not only in the sense of the limit of a sequence of random variables, but also in the sense of the limit of a sequence of stochastic process, and the minimizer of Sn converges to the minimizers of F(α0,β0,α,β)-Ln.
By the convexity of the functionF(α0,β0,α,β), we have
(6.18)β^0=∑i=1n[μ1(1)-((Xi-x0)/h)](1/nh2)K1,h(Xi-x0)K2,h(Ui-u0)sgn(εi)2g(0∣x0,u0)f(x0,u0)[μ12(1)-μ2(1)]+[μ2(1)μ1(1)-μ3(1)]a′′(x0)2[μ12(1)-μ2(1)],β^=Ω-1∑i=1n[μ1(2)-((Ui-u0)/h)](Zi/nh2)K1,h(Xi-x0)K2,h(Ui-u0)sgn(εi)2g(0∣x0,u0)f(x0,u0)[μ12(2)-μ2(2)]+[μ2(2)μ1(2)-μ3(2)]c(u0)2[μ12(2)-μ2(2)],α^=[μ22(2)-μ1(2)μ3(2)]c(u0)2[μ2(2)-μ12(2)]+12nh2∑i=1n‍(Ω-ωω′)-1(Zi-ω)K1,h(Xi-x0)K2,h(Ui-u0)sgn(εi)-μ1(2)Ω-1∑i=1n[μ1(2)-((Ui-u0)/h)]Zi(1/nh2)K1,h(Xi-x0)K2,h(Ui-u0)sgn(εi)2g(0∣x0,u0)f(x0,u0)[μ12(2)-μ2(2)]α^0=[μ22(1)-μ3(1)μ1(1)]a′′(x0)2[μ2(1)-μ12(1)]+12nh2∑i=1n‍(1-ω′Ω-1Zi)K1,h(Xi-x0)K2,h(Ui-u0)sgn(εi)-μ1(1)∑i=1n[μ1(1)-((Xi-x0)/h)](1/nh2)K1,h(Xi-x0)K2,h(Ui-u0)sgn(εi)2g(0∣x0,u0)f(x0,u0)[μ12(1)-μ2(1)].
Thus we obtain(6.19)α~(u0)h=1nh∑k=1n‍α^(Xk,u0)=12nhnh∑k=1n‍∑i=1n‍(Ω-ωω′)-1(Zi-ω)K1,h(Xi-Xk)K2,h(Ui-u0)sgn(εi)-μ1(2)Ω-1∑k=1n∑i=1n[μ1(2)-((Ui-u0)/h)]ZiK1,h(Xi-Xk)K2,h(Ui-u0)sgn(εi)2nhnhg(0∣x0,u0)f(x0,u0)[μ12(2)-μ2(2)]+[μ22(2)-μ1(2)μ3(2)]c(u0)2h[μ2(2)-μ12(2)]:=A1+[μ22(2)-μ1(2)μ3(2)]c(u0)2h[μ2(2)-μ12(2)].
By interchanging summation signs and noting that
(6.20)1n∑k=1n‍K1,h(Xi-Xk)=hf1(Xi){1+op(1)}A1 can be rewritten as
(6.21)A1=12nh∑i=1n‍(Ω-ωω′)-1(Zi-ω)K2,h(Ui-u0)f1(Xi)sgn(εi){1+op(1)}-μ1(2)Ω-1∑i=1n[μ1(2)-((Ui-u0)/h)]ZiK2,h(Ui-u0)f1(Xi)sgn(εi){1+op(1)}2nhg(0∣x0,u0)f(x0,u0)[μ12(2)-μ2(2)].
For all t∈Rp, by some straightforward computations, we get that the mean and variance of t′A1 are
(6.22)E(t′A1)=0E((t′A1)2)=t′{(Ω-ωω′)-1γ0(2)r22(u0)4t′1+μ12(2)[μ12(2)γ0(2)-2μ1(2)γ1(2)+γ2(2)]r12(u0)Ω-14g2(0∣x0,u0)f2(x0,u0)[μ12(2)-μ2(2)]2t′1-μ1(2)[μ1(2)γ1(2)-γ1(2)]r12(u0)2g(0∣x0,u0)f(x0,u0)[μ12(2)-μ2(2)]{∑j=1p‍aj′′(ηij)(Ui-u0)2Zij}}t{1+op(1)}.
By Assumptions (1)–(7), and using the methods used in the proof of Theorem 1 in Tang and Wang [19], we can easily check that the Lindeberg-Feller’s condition for t′A1 is held. So Cramer-Wold theorem tells us (3.8) follows. This completes the proof of Theorem 3.3.



---

*Source: 131085-2012-12-26.xml*