Model and Variable Selection Procedures for Semiparametric Time Series Regression

Abstract

Semiparametric regression models are very useful for time series analysis. They facilitate the detection of features resulting from external interventions. The complexity of semiparametric models poses new challenges for issues of nonparametric and parametric inference and model selection that frequently arise from time series data analysis. In this paper, we propose penalized least squares estimators which can simultaneously select significant variables and estimate unknown parameters. An innovative class of variable selection procedure is proposed to select significant variables and basis functions in a semiparametric model. The asymptotic normality of the resulting estimators is established. Information criteria for model selection are also proposed. We illustrate the effectiveness of the proposed procedures with numerical simulations.

1. Introduction

Non- and semiparametric regression has become a rapidly developing field of statistics in recent years. Various types of nonlinear model such as neural networks, kernel methods, as well as spline method, series estimation, local linear estimation have been applied in many fields. Non- and semiparametric methods, unlike parametric methods, make no or only mild assumptions about the trend or seasonal components and are, therefore, attractive when the data on hand does not meet the criteria for classical time series models. However, the price of this flexibility can be high; when multiple predictor variables are included in the regression equation, nonparametric regression faces the so-called curse of dimensionality.

A major problem associated with non- and semiparametric trend estimation involves the selection of a smoothing parameter and the number of basis functions. Most literature on nonparametric regression with dependent errors focuses on the kernel estimator of the trend function (see, e.g., Altman [1], Hart [2] and Herrmann et al. [3]). These results have been extended to the case with long-memory errors by Hall and Hart [4], Ray and Tsay [5], and Beran and Feng [6]. Kernel methods are affected by the so-called boundary effect. A well-known estimator with automatic boundary correction is the local polynomial approach which is asymptotically equivalent to some kernel estimates. For detailed discussions on local polynomial fitting see, for example, Fan and Gijbels [7] and Fan and Yao [8].

For semiparametric models with serially correlated errors, Gao [9] proposed the semiparametric least-square estimators (SLSEs) for the parametric component and studied its asymptotic properties. You and Chen [10] constructed a semiparametric generalized least-square estimator (SGLSE) with autoregressive errors. Aneiros-Pérez and Vilar-Fernández [11] constructed SLSE with correlated errors.

Like parametric regression models, variable selection of the smoothing parameter for the basis functions is important problem in non- and semiparametric models. It is common practice to include only important variables in the model to enhance predictability. The general approach to finding sensible parameters is to choose an optimal subset determined according to the model selection criterion. Several information criteria for evaluating models constructed by various estimation procedures have been proposed, see, for example, Konishi and Kitagawa [12]. The commonly used criteria are generalized cross-validation, the Akaike information criterion (AIC), and the Bayesian information criterion (BIC). Although best subset selection is practically useful, these selection procedures ignore stochastic errors inherited between the stages of variable selection. Furthermore, best subset selection lacks stability, see, for example, Breiman [13]. Nonconcave penalized likelihood approaches for selecting significant variables for parametric regression models have been proposed by Fan and Li [14]. This methodology can be extended to semiparametric generalized regression models with dependent errors. One of the advantages of this procedure is the simultaneous selection of variables and the estimation of unknown parameters.

The rest of this paper is organized as follows. In Section 2.1 we introduce our semiparametric regression models and explain classical partial ridge regression estimation. Rather than focus on the kernel estimator of the trend function, we use the basis functions to fit the trend component of time series. In Section 2.2, we propose a penalized weighted least-square approach with information criteria for estimation and variable selection. The estimation algorithms are explained in Section 2.3. In Section 2.4, the GIC proposed by Konishi and Kitagawa [15], the BIC𝑚 proposed by Hastie and Tibshirani [16], and the BIC𝑝 proposed by Konishi et al. [17] are applied to the evaluation of models estimated by penalized weighted least-square. Section 2.5 contains the asymptotic results of proposed estimators. In Section 3 the performance of these information criteria is evaluated by simulation studies. Section 4 contains the real data analysis. Section 5 concludes our results, and proofs of the theorems are given in the appendix.

2. Estimation Procedures

In this section, we present our semiparametric regression model and estimation procedures.

2.1. The Model and Penalized Estimation

We consider the semiparametric regression model:

𝑦_𝑖=𝛼(𝑡_𝑖)+𝜷^′𝐱_𝑖+𝜀_𝑖,𝑖=1,…,𝑛, (2.1) where 𝑦_𝑖 is the response variable and 𝐱_𝑖 is the 𝑑×1 covariate vector at time 𝑖, 𝛼(𝑡_𝑖) is an unspecified baseline function of 𝑡_𝑖 with 𝑡_𝑖=𝑖/𝑛, 𝜷 is a vector of unknown regression coefficients, and 𝜀_𝑖 is a Gaussian and zero mean covariance stationary process.

We assume the following properties for the error terms 𝜀_𝑖 and vectors of explanatory variables 𝐱_𝑖.

(A.1) It holds that {𝜀_𝑖} is a linear process given by 𝜀_𝑖=^∞_𝑗=0𝑏_𝑗𝑒_𝑖−𝑗, (2.2) where 𝑏₀=1 and {𝑒_𝑖} is an i.i.d. Gaussian random variable with 𝐸{𝑒_𝑖}=0 and 𝐸{𝑒_𝑖²}=𝜎²_𝑒.(A.2)The coefficients 𝑏_𝑗 satisfy the conditions that for all |𝑧|<1, ∑^∞_𝑗=0𝑏_𝑗𝑧^𝑗≠0 and ∑^∞_𝑗=0𝑗²|𝑏_𝑗|<∞.

We define 𝛾(𝑘)=cov(𝜀_𝑡,𝜀_𝑡+𝑘)=𝐸{𝜀_𝑡𝜀_𝑡+𝑘}.

The assumptions on covariate variables are as follows.

(B.1) Also 𝐱_𝑖=(𝑥_𝑖1,…,𝑥_𝑖𝑑)^′∈ℝ^𝑑 and {𝑥_𝑖𝑗},𝑗=1,…𝑑, have mean zero and variance 1.

The trend function 𝛼(𝑡_𝑖) is expressed as a linear combination of a set of 𝑚 underlying basis functions :

𝛼(𝑡_𝑖)=^𝑚_𝑘=1𝑤_𝑘𝜙_𝑘(𝑡_𝑖)=𝐰^′𝝓(𝑡), (2.3) where {𝝓(𝑡_𝑖)=(𝜙₁(𝑡_𝑖),…,𝜙_𝑚(𝑡_𝑖))^′} is an 𝑚-dimensional vector constructed from basis functions {𝜙_𝑘(𝑡_𝑖);𝑘=1,…,𝑚}, and 𝐰=(𝑤₁,…,𝑤_𝑚)^′ is an unknown parameter vector to be estimated. The examples of basis functions are B-spline, P-spline, and radial basis functions. A P-spline basis is given by

𝝓(𝑡_𝑖)=(𝑡_𝑖,…,𝑡^𝑝_𝑖,(𝑡_𝑖−𝜅₁)^𝑝₊,…,(𝑡_𝑖−𝜅_𝑘)^𝑝₊)^′, (2.4) where {𝜅_𝑘}_{𝑘=1,…,𝐾} are spline knots. This specification uses the so-called truncated power function basis. The choice of the number of knots 𝐾 and the knot locations are discussed by Yu and Ruppert [18].

Radial basis function (RBF) emerged as a variant of artificial renewal network in late 80s. Nonlinear specification of using RBF has been widely used in cognitive science, engineering, biology, linguistics, and so on. If we consider the RBF modeling, a basis function can take the form

𝜙_𝑘(𝑡_𝑖)=exp(−‖𝑡_𝑖−𝜇_𝑘‖²2𝑠²_𝑘), (2.5) where 𝜇_𝑘 determines the location and 𝑠²_𝑘 determines the width of the basis function.

Selecting appropriate basis functions, then the semiparametric regression model (2.1) can be expressed as a linear model

𝐲=𝐗𝜷+𝐁𝐰+𝜺, (2.6) where 𝐗=(𝐱₁,…,𝐱_𝑛)^′, 𝐲=(𝑦₁,…,𝑦_𝑛)^′, 𝐁=(𝝓₁,…,𝝓_𝑛)^′ with 𝝓_𝑖=(𝜙₁(𝑖/𝑛),…,𝜙_𝑚(𝑖/𝑛))^′. The penalized least-square estimator is then a minimizer of the function

12(𝐲−𝐗𝜷−𝐁𝐰)^′(𝐲−𝐗𝜷−𝐁𝐰)+𝑛𝜉𝐰^′𝐊𝐰, (2.7) where 𝜉 is the smoothing parameter controlling the tradeoff between the goodness-of-fit measured by weighted least-square and the roughness of the estimated function. Also 𝐊 is an appropriate positive semidefinite symmetric matrix. For example, if 𝐊 satisfies 𝐰^′𝐊𝐰=∫¹₀[𝛼^′′(𝑢)]²𝑑𝑢, we have the usual quadratic integral penalty (see, e.g., Green and Silverman [19]). By simple calculus, (2.7) is minimized when 𝜷 and 𝐰 satisfy the block matrix equation

𝐗^′𝐗𝐗^′𝐁𝐁^′𝐗𝐁^′𝜷𝐰=𝐗𝐁+𝑛𝜉𝐊𝐁^′′𝐲. (2.8) This equation can be solved without any iteration (see, e.g., Green [20]). First, we find 𝐁𝐰=𝐒(𝐲−𝐗𝜷), where 𝐒=𝐁(𝐁^′𝐁+𝛼𝐊)⁻¹𝐁^′ is usually called the smoothing matrix. Substituting 𝐁𝐰 into (2.6), we obtain

𝐲=𝐗𝜷+𝜺, (2.9) where 𝐲=(𝐈−𝐒)𝐲, 𝐗=(𝐈−𝐒)𝐗 ,and 𝐈 is the identity matrix of order 𝑛. Applying least-square to the linear model (2.9), we obtain the semiparametric ordinary least-square estimator (SOLSE) result:

𝜷_SOLSE=(𝐗^′𝐗)⁻¹𝐗^′𝐲𝐰_SOLSE,=(𝐁^′𝐁+𝑛𝜉𝐊)⁻¹𝐁^′(𝐲−𝐗𝜷_SOLSE). (2.10)(2.11)

Speckman [21] studied similar solutions for partial linear models with independent observations. Since the errors are serially correlated in model (2.1), 𝜷_SOLSE is not asymptotically efficient. To obtain an asymptotically efficient estimator for 𝜷, we use the prewhitening transformation. Note that the errors {𝜀_𝑖} in (2.6) are invertible. Let 𝑏(𝐿)=∑^∞_𝑗=1𝑏_𝑗𝑒_𝑖−𝑗, where 𝐿 is the lag operator and 𝑎(𝐿)=𝑏(𝐿)⁻¹=𝑎₀−∑^∞_𝑗=1𝑎_𝑗𝐿^𝑗 with 𝑎₀=1. Applying 𝑎(𝐿) to the model (2.6) and rewriting the corresponding equation, we obtain the new model:

𝐲=𝐗𝜷+𝐁𝐰+𝐞, (2.12) where 𝐲=(𝑦₁,…,𝑦_𝑛)^′,𝐗=(𝐱₁,…,𝐱_𝑛)^′,𝐁=(𝝓₁,…,𝝓_𝑛)^′ and 𝐞=(𝑒₁,…,𝑒_𝑛)^′. Here 𝑦_𝑖=𝑦_𝑖−^∞_𝑗=1𝑎_𝑗𝑦_𝑖−𝑗,𝝓_𝑖=𝝓_𝑖−^∞_𝑗=1𝑎_𝑗𝝓_𝑖−𝑗,𝐱_𝑖=𝐱_𝑖−^∞_𝑗=1𝑎_𝑗𝐱_𝑖−𝑗. (2.13) The regression errors in (2.12) are i.i.d. Because, in practice, the response variable 𝑦_𝑖 is unknown, we use a reasonable approximation by 𝑦_𝑖 based on the work by Xiao et al. [22] and Aneiros-Pérez and Vilar-Fernández [11].

Under the usual regularity conditions the coefficients 𝑎_𝑗 decrease geometrically so, letting 𝜏=𝜏(𝑛) denote a truncation parameter, we may consider the truncated autoregression on 𝜀_𝑖:

𝑒_𝑖=𝜀_𝑖−^∞_𝑗=1𝑎_𝑗𝜀_𝑖−𝑗, (2.14) where 𝑒_𝑖 are i.i.d. random variables with 𝐸(𝑒_𝑖)=0. We make the following assumption about the truncation parameter.

(C.1)The truncation parameter 𝜏 satisfies 𝜏(𝑛)=𝑐log𝑛 for some 𝑐>0.

The expansion rate of the truncation parameter given in (C.1) is also for convenience. Let 𝐓_𝜏 be the 𝑛×𝑛 transformation matrix such that 𝐞_𝜏=𝐓_𝜏𝜺. Then the model (2.12) can be expressed as

𝐓_𝜏𝐲=𝐓_𝜏𝐗𝜷+𝐓_𝜏𝐁𝐰+𝐓_𝜏𝜺, (2.15) where

𝐓_𝜏=𝛿₁₁𝛿0⋯0₂₁−𝛿₂₂⋮𝛿0⋯0_𝜏1⋯−𝛿_𝜏𝜏−𝑎_𝜏⋯−𝑎₁10−𝑎_𝜏⋯−𝑎₁1⋮0⋯0−𝑎_𝜏⋯−𝑎₁1 (2.16) with 𝛿₁₁=𝜎_𝑒/𝛾(0),𝛿₂₂=𝜎_𝑒/(1−𝜌²(1))𝛾(0),𝛿₂₁=𝜌(1)(𝜎_𝑒/(1−𝜌²(1))𝛾(0)),…. Here 𝜌(ℎ)=𝛾(ℎ)/𝜌(0) denotes the lag ℎ autocorrelation function of {𝜀_𝑖}.

Now our estimation problem for the semiparametric time series regression model can be expressed as the minimization of the function

1ℒ(𝜷,𝐰)=2(𝐲−𝐗𝜷−𝐁𝐰)^′𝐕⁻¹(𝐲−𝐗𝜷−𝐁𝐰)+𝛼𝐰^′𝐊𝐰, (2.17) where 𝐕⁻¹=𝜎_𝑒⁻²𝐓^′_𝜏𝐓_𝜏 and 𝜎_𝑒²=𝑛⁻¹‖𝐓_𝜏𝜺‖². Based on the work by Aneiros-Pérez and Vilar-Fernández [11], an estimator for 𝐓_𝝉 is constructed as follows. We use the residuals 𝜺=𝐲−𝐗𝜷_SOLSE−𝐁𝐰_SOLSE to construct an estimate of 𝐓_𝜏 using the ordinary least square method applied to the model

𝜀_𝑖̂=𝑎₁𝜀̂_𝑖−1+⋯+𝑎_𝜏𝜀̂_𝑖−𝜏+residual_𝑖. (2.18) Define the estimate 𝐚_𝜏=(𝑎₁,𝑎₂,…,𝑎_𝜏)^′ of 𝐚_𝜏=(𝑎₁,𝑎₂,…,𝑎_𝜏)^′, where 𝐚_𝜏=(𝐄^′_𝜏𝐄_𝜏)⁻¹𝐄^′_𝜏𝜺, (2.19) where 𝜺=(𝜀̂_𝜏+1,…,𝜀_𝑛̂) and 𝐄_𝜏 is the (𝑛−𝜏)×𝜏 matrix of regressors with the typical element 𝜀̂_𝑖−𝑗. Then 𝐓_𝜏 is obtained from 𝐓_𝜏 by replacing 𝑎_𝑗 with 𝑎_𝑗, 𝜎²_𝑒 with 𝜎²_𝑒, and so forth. Applying least-square to the linear model, we obtain 𝐓_𝜏𝐲=𝐓_𝜏𝐗𝜷+𝐓_𝜏𝐁𝐰+𝐓_𝜏𝜺. (2.20) Then

𝜷_SGLSE=(𝐗^′𝜏𝐗𝜏)⁻¹𝐗𝜏^′𝐲𝜏,𝐰_SGLSE=(𝐁^′𝜏𝐁𝜏+𝑛𝜉𝐊)⁻¹𝐁^′𝜏(𝐲𝜏−𝐗𝜏𝜷_SGLSE), (2.21) where 𝐗𝜏=(𝐈−𝐒)𝐗𝜏 and 𝐲𝜏=(𝐈−𝐒)𝐲𝜏, with 𝐲𝜏=𝐓_𝜏𝐲 and 𝐗𝜏=𝐓_𝜏𝐗. The following theorem shows that the loss in efficiency associated with the estimation of the autocorrelation structure is modest in large samples.

Theorem 2.1. Let the conditions of (A.1), (A.2), (B.1), and (C.1) hold, and assume that 𝚺₁=lim_𝑛→∞𝑛⁻¹𝐗^′𝐕⁻¹𝐗 is nonsingular. Let 𝜷₀ denote the true value of 𝜷, then √𝑛(𝜷−𝜷₀)=√𝑛(𝜷_SGLSE−𝜷₀)+𝑂_𝑝𝜏((𝑛)^1/2)√,𝑛(𝜷_SGLSE−𝜷₀)→_𝐷𝑁(0,𝚺₁⁻¹), (2.22)(2.23) where →_𝐷 denotes convergence in distribution and 𝜷=(𝐗^′_𝜏𝐗_𝜏)⁻¹𝐗^′_𝜏𝐲_𝜏. Assume that 𝚺₂=lim_𝑛→∞𝑛⁻¹𝐁^′𝐕⁻¹𝐁 is nonsingular and let 𝐰_𝟎 denote the true value of 𝐰, then one has √𝑛(𝐰−𝐰₀)=√𝑛(𝐰_SGLSE−𝐰₀)+𝑂_𝑝𝜏((𝑛)^1/2)√,𝑛(𝐰_SGLSE−𝐰₀)→_𝐷𝑁(0,𝚺₂⁻¹), (2.24)(2.25) where 𝐰=(𝐁^′_𝜏𝐁_𝜏+𝑛𝜉𝐊)⁻¹𝐁^′_𝜏(𝐲_𝜏−𝐗_𝜏𝜷).

2.2. Variable Selection and Penalized Least Squares

Variable and model selection are an indispensable tool for statistical data analysis. However, it has rarely been studied in the semiparametric context. Fan and Li [23] studied penalized weighted least-square estimation with variable selection in semiparametric models for longitudinal data. In this section, we introduce the penalized weighted least-square approach. We propose an algorithm for calculating the penalized weighted least-square estimator of 𝜽=(𝜷^′,𝐰^′)^′ in Section 2.3. In Section 2.4 we present the information criteria for the model selection.

From now on, we assume that the matrices 𝐗_𝜏 and 𝐁_𝜏 are standardized so that each column has mean 0 and variance 1. The first term in (2.7) can be regarded as a loss function of 𝜷 and 𝐰, which we will denote by 𝑙(𝜷,𝐰). Then expression (2.7) can be written as

ℒ(𝜷,𝐰)=𝑙(𝜷,𝐰)+𝑛𝜉𝐰^′𝐊𝐰. (2.26) The methodology in the previous section can be applied to the variable selection via penalized least-square. A form of penalized weighted least-square is

𝒮(𝜷,𝐰)=𝑙(𝜷,𝐰)+𝑛(^𝑑_𝑖=1𝑝_𝜆1(|𝛽_𝑖|)+^𝑚_𝑗=1𝑝_𝜆2(|𝑤_𝑗|))+𝑛𝜉𝐰^′𝐊𝐰. (2.27) where 𝑝_𝜆𝑖(⋅) are penalty functions and 𝜆_𝑖 are regularization parameters, which control the model complexity. By minimizing (2.27) with a special construction of the penalty function given in what following some coefficients are estimated as 0, which deletes the corresponding variables, whereas others are not. Thus, the procedure selects variables and estimates coefficients simultaneously. The resulting estimate is called a penalized weighted least-square estimate.

Many penalty functions have been used for penalized least-square and penalized likelihood in various non- and semiparametric models. There are strong connections between the penalized weighted least-square and the variable selection. Denote by 𝜽=(𝜷^′,𝐰^′)^′ and 𝐳=(𝑧₁,…,𝑧_𝑑+𝑚)^′ the true parameters and the estimates, respectively. By taking the hard thresholding penalty function

𝑝_𝜆(|𝜃|)=𝜆²+(|𝜃|−𝜆)²𝐼(|𝜃|<𝜆), (2.28) we obtain the hard thresholding rule

𝜃=𝑧𝐼(|𝑧|>𝜆). (2.29) The 𝐿₂ penalty 𝑃_𝜆(|𝜃|)=𝜆|𝜃|² results in a ridge regression and the 𝐿₁ penalty 𝑃_𝜆(|𝜃|)=𝜆|𝜃| yields a soft thresholding rule

𝜃=sgn(𝑧)𝐼(|𝑧|>𝜆)₊. (2.30) This solution gives the best subset selection via stepwise deletion and addition. Tibshirani [24, 25] has proposed LASSO, which is the penalized least-square estimate with the 𝐿₁ penalty, in the general least-square and likelihood settings.

2.3. An Estimation Algorithm

In this section we describe an algorithm for calculating the penalized least-square estimator of 𝜽=(𝜷^′,𝐰^′)^′. The estimate of 𝜽 minimizes the penalized sum of squares ℒ(𝜽) given by (2.17). First we obtain 𝜽_SOLSE in Step 1. In Step 2, we estimate 𝐓_𝜏 by using 𝜀 obtained in Step 1. Then 𝜽^HT_SGLSE is obtained using 𝐓_𝜏(Step 3). Here the penalty parameters 𝜆, and 𝜉, and the number of basis functions 𝑚 are chosen using information criteria that will be discussed in Section 2.4.

Step 1. First we obtain 𝜷_SOLSE and 𝐰_SOLSE by (2.10) and (2.11), respectively. Then we have the model 𝐲=𝐁𝐰_SOLSE+𝐗𝜷_SOLSE+𝜺. (2.31)

Step 2. An estimator for 𝐓_𝜏 is constructed followings the work of Aneiros-Pérez and Vilar-Fernández [4]. We use the residuals 𝜺=𝐲−𝐁𝐰_SOLSE−𝐗𝜷_SOLSE to construct an estimate of 𝐓_𝜏 using the ordinary least square method applied to the model 𝜀_𝑖̂=𝑎₁𝜀̂_𝑖−1+⋯+𝑎_𝜏𝜀̂_𝑖−𝜏+residual_𝑖. (2.32) The estimator 𝐓_𝜏 is obtained from 𝐓_𝜏 by replacing parameters with their estimates.

Step 3. Our SGLSE of 𝜽 is obtained by using the model 𝐲𝜏=𝐁𝜏𝐰+𝐗𝜏𝜷+𝜺𝜏, (2.33) where 𝐲𝜏=𝐓_𝜏𝐲, 𝐁𝜏=𝐓_𝜏𝐁, 𝐗𝜏=𝐓_𝜏𝐗, and 𝜺_𝜏=𝐓_𝜏𝜺. Finding the solution of the penalized least-square of (2.27) needs the local quadratic approximation, because the 𝐿₁ and hard thresholding penalty are irregular at the origin and may not have second derivatives at some points. We follow the methodology of Fan and Li [14]. Suppose that we are given an initial value 𝜽⁽⁰⁾ that is close to the minimizer of (2.27). If 𝜃_𝑗⁽⁰⁾ is very close to 0, then set 𝜃_𝑗⁽⁰⁾=0. Otherwise they can be locally approximated by a quadratic function as [𝑝_𝜆𝑗(𝜃_𝑗)]^′=𝑝^′_𝜆𝑗(|𝜃_𝑗|)sgn(𝜃_𝑗𝑝)≈{^′𝜆𝑗(|𝜃_𝑗⁽⁰⁾|)|𝜃_𝑗⁽⁰⁾|}𝜃_𝑗, (2.34) when 𝜃_𝑗⁽⁰⁾≠0. Therefore, the minimization problem (2.27) can be reduced to a quadratic minimization problem and the Newton-Raphson algorithm can be used. The right-hand side of equation (2.27) can be locally approximated by 𝑙(𝜷₀,𝐰₀)+∇𝑙_𝛽(𝜷,𝐰)^′(𝜷−𝜷₀)+∇𝑙_𝑤(𝜷₀,𝐰₀)^′(𝐰−𝐰₀)+12(𝜷−𝜷₀)^′∇²_𝛽𝛽(𝜷₀,𝐰₀)(𝜷−𝜷₀1)+2(𝐰−𝐰₀)^′∇²_𝑤𝑤(𝜷₀,𝐰₀)(𝐰−𝐰₀)12(𝜷₀−𝜷)^′∇²_𝛽𝑤(𝐰−𝐰₀)+𝑛𝜷^′𝚺_𝜆1(𝜷₀)𝜷+𝑛𝐰^′𝚺_𝜆2(𝐰_𝟎)𝐰,, (2.35) where ∇𝑙_𝛽(𝜷₀,𝐰₀)=𝜕𝑙(𝜷₀,𝐰₀)𝜕𝜷,∇𝑙_𝑤(𝜷₀,𝐰₀)=𝜕𝑙(𝜷₀,𝐰₀),∇𝜕𝐰₂𝑙_𝛽𝛽(𝜷₀,𝐰₀𝜕)=²𝑙(𝜷₀,𝐰₀)𝜕𝜷𝜕𝜷^′,∇²𝑙_𝑤𝑤(𝜷₀,𝐰₀𝜕)=²𝑙(𝜷₀,𝐰₀)𝜕𝐰𝜕𝐰^′,∇²𝑙_𝛽,𝑤(𝜷₀,𝐰₀𝜕)=²𝑙(𝜷₀,𝐰₀),𝚺𝜕𝜷𝜕𝐰_𝜆1(𝜷₀𝑝)=diag{^′𝜆1(|𝛽₁⁽⁰⁾|)|𝛽₁⁽⁰⁾|𝑝,…,^′𝜆1(|𝛽_𝑑⁽⁰⁾|)|𝛽_𝑑⁽⁰⁾|𝚺},_𝜆2(𝐰₀𝑝)=diag{^′𝜆2(|𝑤₁⁽⁰⁾|)|𝑤₁⁽⁰⁾|𝑝,…,^′𝜆2(|𝑤_𝑚⁽⁰⁾|)|𝑤_𝑚⁽⁰⁾|}. (2.36) The solution can be found by iteratively computing the block matrix equation: 𝐗𝜏^′𝐗𝜏+𝑛𝚺_𝜆1(𝜷⁽⁰⁾)𝐗𝜏^′𝐁𝜏𝐁𝜏^′𝐗𝜏𝐁𝜏^′𝐁𝜏+𝛼𝐊+𝑛𝚺_𝜆2(𝐰⁽⁰⁾)𝜷𝐰=𝐗𝜏^′𝐁𝜏^′𝐲. (2.37) This gives the estimators 𝜷^HT_SGLSE=(𝐗𝜏^′𝐗𝜏+𝑛𝚺_𝜆1(𝜷⁽⁰⁾))⁻¹𝐗𝜏^′𝐲𝜏,𝐰^HT_SGLSE=(𝐁𝜏^′𝐁𝜏+𝑛𝜉𝐊+𝑛𝚺_𝜆2(𝐰⁽⁰⁾))⁽⁻¹⁾𝐁𝜏^′(𝐲𝜏−𝐗𝜏𝜷^HT_SGLSE), (2.38) where 𝐲𝜏=(𝐈−𝐒𝜏)𝐲𝜏, 𝐗𝜏=(𝐈−𝐒𝜏)𝐗𝜏, and 𝐒𝜏=𝐁𝜏(𝐁𝜏^′𝐁𝜏+𝑛𝜉𝐊+𝑛𝚺_𝜆2(𝐰⁽⁰⁾))⁻¹𝐁𝜏^′.

2.4. Information Criteria

Selecting suitable values for the penalty parameters and number of basis functions is crucial to obtaining good curve fitting and variable selection. The estimate of 𝜽 minimizes the penalized sum of squares ℒ(𝜽) given by (2.17). In this section, we express the model (2.15) as

𝐲_𝜏=𝐀_𝜏𝜽+𝐞, (2.39) where 𝐀_𝜏=(𝐗_𝜏,𝐁_𝜏) and 𝜽=(𝜷^′,𝐰^′)^′. In many applications, the number of basis functions 𝑚 needs to be large to adequately capture the trend. To determine the number of basis functions, all models with 𝑚≤𝑚_max are fitted and the preferred model minimizes some model selection criteria.

The Schwarz BIC is given by

BIC=𝑛log(2𝜋𝜎²_𝑒)+log𝑛(thenumberofparameters), (2.40) where 𝜎²_𝑒 is the least-square estimate of 𝜎²_𝑒 without a degree of freedom correction. Hastie and Tibshirani [16] used the trace of the smoother matrix as an approximation to the effective number of parameters. By replacing the number of parameters in BIC by tr𝐒_𝛽, we formally obtain information criteria for the basis function Gaussian regression model in the form

BIC𝑚=𝑛log(2𝜋𝜎²_𝑒)+(tr𝐒_𝜽)log𝑛, (2.41) where 𝜎²_𝑒=𝑛⁻¹‖𝑦−𝐒_𝜽𝑦‖² and

tr𝐒_𝜽=𝐀_𝜏(𝐀^′_𝜏𝐀_𝜏+𝑛𝜉𝐊+𝑛𝚺_𝜆(𝜽))⁻¹𝐀^′_𝜏. (2.42) Here 𝚺_𝜆(𝜽) is defined by (2.44) in what follows.

We also consider the use of the BIC𝑝 criterion to choose appropriate values for these unknown parameters. Denote 𝚺_𝜆1(𝜷)=diag{𝑝_𝜆^′′₁(|𝛽₁₀|),…,𝑝_𝜆^′′₁(|𝛽_𝑑0𝚺|)},_𝜆2(𝐰)=diag{𝑝_𝜆^′′₂(|𝑤₁₀|),…,𝑝_𝜆^′′₂(|𝑤_𝑚0𝚺|)},_𝜆(𝜽)=(𝚺_𝜆1(𝜷),𝚺_𝜆2(𝐰)). (2.43)(2.44) Let 𝑁₁ and 𝑁₂ be the number of zero components in 𝜷₀ and 𝐰₀, respectively. Then the BIC𝑝 criterion is

BIC𝑝=𝑛log(2𝜋𝜎²_𝑒)+𝑛𝜽^′𝚺_𝜆(𝜽)𝜽+𝑛𝜉𝜽𝐊𝜽+log|𝐉_𝐺(𝜽)|−log|𝐊|₊−log|Σ_𝜆(𝜽)|₊−(𝑚−𝑁₂)log𝜉+Const, (2.45) where 𝐉_𝐺(𝜽) is the (𝑑+𝑚+1)×(𝑑+𝑚+1) matrix of second derivatives of the penalized likelihood defined by

𝐉_𝐺=1𝑛𝜎²_𝑒𝐀^′_𝜏𝐀_𝜏+𝑛Σ_𝜆𝐊𝐀(𝜽)+𝑛𝜉^′𝜏Λ𝟏_𝑛𝜎²_𝑒𝟏^′_𝑛Λ𝐀_𝜏𝜎²_𝑒𝑛2𝜎²_𝑒. (2.46) Here Λ is a diagonal matrix with 𝑖th element Λ_𝑖=diag[𝑒₁,…,𝑒_𝑛] and 𝟏_𝑛=(1,…,1)^′. The 𝑛-dimensional vector 𝐪 has 𝑖th element (𝐓_𝑖𝑗𝐲_𝑗−𝐀_{𝜏,𝑖𝑗}𝜽_𝑗)²/2𝜎⁴_𝑒−1/2𝜎²_𝑒 where 𝐓_𝑖𝑗 is the element in the 𝑖th row and 𝑗th column of 𝐓_𝜏. Also 𝐊 is the (𝑑+𝑚)×(𝑑+𝑚) matrix defined by

𝐊=𝐊𝐎_𝑑,𝑚𝐎_𝑚,𝑑𝐎_𝑚,𝑚, (2.47) and |𝐊|₊ and |Σ_𝜆(𝜽)|₊ are the product of the (𝑚−𝑁₁) and (𝑑+𝑚−𝑁₁−𝑁₂) nonzero eigenvalues of 𝐊 and Σ_𝜆(𝜽), respectively.

Konishi and Kitagawa [15] proposed a framework of Generalized Information Criteria (GIC) to the case where the models are not estimated by maximum likelihood. Hence, we also consider the use of GIC for the model evaluations. The GIC for the hard thresholding penalty function is given by

GIC=𝑛log(2𝜋𝜎²_𝑒)+𝑛+2tr{𝐈_𝐺𝐉_𝐺⁻¹}, (2.48) where 𝐈_𝐺 is a (𝑚+𝑑+1)×(𝑚+𝑑+1) matrix. Also 𝐈_𝐺 is basically the product of the empirical influence function and the score function. It is defined by 𝐈_𝐺=1𝑛𝜎²_𝑒𝐀^′_𝜏Λ𝜎²_𝑒−𝚺_𝜆(𝜽)𝜽𝟏^′_𝐧−𝜉𝐊𝜽𝟏^′_𝐧𝐪^′(Λ𝐀_𝜏,𝜎²_𝑒𝐪). (2.49)

The number of basis functions 𝑚, penalty parameters 𝜉,𝜆₁,𝜆₂ are determined by minimizing BIC𝑚, BIC𝑝 or GIC.

2.5. Sampling Properties

We now study the asymptotic properties of the estimate resulting from the penalized least-square function (2.27).

First we establish the convergence rate of the penalized profile least-square estimator. Assume that penalty functions 𝑝^′_𝜆1𝑗(⋅) and 𝑝^′_𝜆2𝑗(⋅) are negative and nondecreasing with 𝑝^′_𝜆1𝑗(0)=𝑝^′_𝜆2𝑗(0)=0. Let 𝜷₀ and 𝐰₀ denote the true values of 𝜷 and 𝐰, respectively. Also let

𝑎_1𝑛=max_𝑗{|𝑝^′_𝜆1𝑗(|𝛽_𝑗0|)|∶𝛽_𝑗0≠0},𝑎_2𝑛=max_𝑗{|𝑝^′_𝜆2𝑗(|𝑤_𝑗0|)|∶𝑤_𝑗0𝑏≠0},_1𝑛=max_𝑗{|𝑝_𝜆^′′_1𝑗(|𝛽_𝑗0|)|∶𝛽_𝑗0≠0},𝑏_2𝑛=max_𝑗{|𝑝_𝜆^′′_2𝑗(|𝑤_𝑗0|)|∶𝑤_𝑗0≠0}. (2.50)

Theorem 2.2. Under the conditions of Theorem 2.1, if 𝑎_1𝑛,𝑏_1𝑛,𝑎_2𝑛, and 𝑏_2𝑛 tend to 0 as 𝑛→∞, then with probability tending to 1, there exist local minimizers 𝜷 and 𝐰 of ℒ(𝜷,𝐰) such that ‖𝜷^HT_SGLSE−𝜷₀‖=𝑂_𝑝(𝑛^−1/2+𝑎_1𝑛) and ‖𝐰^HT_SLOSE−𝐰₀‖=𝑂_𝑝(𝑛^−1/2+𝑎_2𝑛).

Theorem 2.2 demonstrates how the rate of convergence of the penalized least-square estimator 𝜽^HT_SGLSE=(𝜷^'HT_SGLSE,𝐰^'HT_SGLSE)^′ of ℒ(𝜽) depends on 𝜆_𝑖𝑗 for 𝑖=1,2. To achieve the root 𝑛 convergence rate, we have to take 𝜆_𝑖𝑗 small enough so that 𝑎_𝑛=𝑂_𝑝(𝑛^−1/2).

Next we establish the oracle property for the penalized least-square estimator. Let 𝜷_𝑆1 consist of all nonzero components of 𝜷₀ and let 𝜷_𝑁1 consist of all zero components. Let 𝐰_𝑆2 consist of all nonzero components of 𝐰₀ and let 𝐰_𝑁2 consist of all zero components. Let

𝐱(𝑡)^′𝜷₀=𝐱^′_𝑆1𝜷_𝑆1+𝐱^′_𝑁1(𝑡)𝜷_𝑁1=𝐱_𝑆1(𝑡)^′𝜷_𝑆1,𝝓(𝑡)^′𝐰₀=𝜙^′_𝑆2𝐰_𝑆2+𝜙^′_𝑁2(𝑡)𝐰_𝑁2=𝜙_𝑆2(𝑡)^′𝐰_𝑆2. (2.51) Write

𝐛_𝛽=(𝑝^′_𝜆1𝑛(|𝛽₁₀|)sgn(𝛽₁₀),…,𝑝^′_{𝜆𝑆1𝑛}(|𝛽_𝑆10|)sgn(𝛽_𝑆10))^′,𝐛_𝑤=(𝑝^′_𝜆2𝑛(|𝑤₁₀|)sgn(𝑤₁₀),…,𝑝^′_{𝜆𝑆2𝑛}(|𝑤_𝑆20|)sgn(𝑤_𝑆20))^′. (2.52) Further, let 𝜷₁ consist of the first 𝑆₁ components of 𝜷 and let 𝜷₂ consist of the last 𝑑−𝑆₁ components of 𝜷^HT_SGLSE. Let 𝐰_𝟏 consist of the first 𝑆₂ components of 𝐰 and let 𝐰₂ consist of the last 𝑚−𝑆₂ components of 𝐰^HT_SGLSE.

Theorem 2.3. Assume that for 𝑗=1,…,𝑑 and 𝑘=1,…,𝑚, one has 𝜆₁→0, √𝑛𝜆₁→∞, 𝜆₂→0 and √𝑛𝜆₂→∞. Assume that the penalty functions 𝑝^′_𝜆1(|𝛽_𝑗|) and 𝑝^′_𝜆2(|𝑤_𝑘|) satisfy liminf_𝑛→∞liminf_{𝛽𝑗→0+}𝑝^′_𝜆1(𝛽_𝑗)𝜆₁>0,liminf_𝑛→∞liminf_{𝜔𝑘→0+}𝑝^′_𝜆2(𝜔_𝑘)𝜆₂>0. (2.53) If 𝑎_1𝑛=𝑎_2𝑛=𝑂_𝑝(𝑛^−1/2) then, under the conditions of Theorem 2.1, with probability tending to 1, the root 𝑛 consistent local minimizers 𝜷^HT_SGLSE=(𝜷^′₁,𝜷^′₂)^′ and 𝐰^HT_SGLSE=(𝐰^′₁,𝐰^′₂)^′ in Theorem 2.2 must satisfy the following: (1)(sparsity) 𝜷₂=𝐰₂=𝟎;(2)(asymptotic normality) √𝑛(𝐈_𝑆1+𝚺_𝜆1(𝜷))(𝜷₁−𝜷₁₀+(𝐈_𝑆1+𝚺_𝜆1(𝜷))⁻¹𝐛_𝛽)→𝑁_𝑆1(0,𝚺⁻¹₁₍₁₎√),𝑛(𝐈_𝑆2+𝚺_𝜆1(𝐰))(𝐰₁−𝐰₁₀+(𝐈_𝑆2+𝚺_𝜆2(𝐰)+𝜉𝐊)⁻¹𝐛_𝑤)→𝑁_𝑆2(0,𝚺⁻¹₂₍₁₎) (2.54)Here 𝚺⁻¹₁₍₁₎ and 𝚺⁻¹₂₍₁₎ consist of the first 𝑆₁ and 𝑆₂ rows and columns of 𝚺₁ and 𝚺₂ defined in Theorem 2.1, respectively.

3. Numerical Simulations

We now assess the performance of semiparametric estimators of the proposed in previous section via simulations. We generate simulation data from the model

𝑦_𝑖=𝛼(𝑡_𝑖)+𝜷^′𝐱_𝑖+𝜀_𝑖, (3.1) where 𝛼(𝑡_𝑖)=exp(−3(𝑖/𝑛))sin(3𝜋𝑖/𝑛), 𝜷=(3,1.5,0,0.2,0,0,0)^′ and 𝜀(𝑡) is a Gaussian AR(1) process with autoregressive coefficient 𝜌. We used the radial basis function network modeling to fit the trend component. We simulate the covariate vector 𝑥 from a normal distribution with mean 0 and cov(𝑥_𝑖,𝑥_𝑗)=0.5^{|𝑖−𝑗|}. In each case, the autoregressive coefficient is set to 0, 0.25, 0.5 or 0.75 and the sample size 𝑛 is set to 50, 100 or 200. Figure 1 depicts some examples of simulation data.

(a)

(b)

(c)

Figure 1 

Simulation data with (a) 𝑛=50 and 𝜌=0.5, (b) 𝑛=100 and 𝜌=0.5, (c) 𝑛=200 and 𝜌=0.5. The dotted lines represent 𝛼(𝑡); the solid lines 𝛼(𝑡)+𝜀(𝑡).

We compare the effectiveness of our proposed procedure (PLS + HT) with an existing procedure (PLS). We also compare the performance of the information criteria BIC𝑚, GIC and BIC𝑝 for evaluating the models. As discussed in Section 3, the proposed procedure (PLS + HT) excludes basis functions as well as explanatory variables.

First we assess the performance of 𝛼(𝑡) by the square root of average squared errors (RASE_𝛼):

RASE_𝛼=𝑛⁻¹_grid^𝑛grid_𝑘=1{𝛼(𝑡_𝑘)−𝛼(𝑡_𝑘)}², (3.2) where {𝑡_𝑘,𝑘=1,…,𝑛_grid} are the grid points at which the baseline function 𝛼(⋅) is estimated. In our simulation, we use 𝑛_grid=200. Table 1 shows the means and standard deviations of RASE_𝛼 for 𝜌=0,0.25,0.5,0.75 based on 500 simulations. RASE_𝛼 increases as the autoregressive coefficient increases but decreases as the sample size increases. From Table 1, we see that the proposed procedure (PLS + HT) works better than PLS and that models evaluated by BIC𝑝 work better than those based on BIC𝑚 or GIC.

Table 1 

Means (standard deviations) of RASE_𝛼.

Then the performance of 𝜷 is assessed by the square root of average squared errors (RASE_𝛽):

RASE_𝛽=1𝑑^𝑑_𝑖=1(𝛽_𝑖−𝛽_𝑖)². (3.3) The means and standard deviations of RASE_𝛽 for 𝜌=0,0.25,0.5,0.75 based on 500 simulations are shown in Table 2. We can see that the proposed procedure (PLS + HT) works better than the existing procedure. There is almost no change in RASE_𝛽 as the autoregressive coefficient changes (unlike the procedure of You and Chen [10]), whereas RASE_𝛽 depends strongly on the information, BIC𝑝 works the best among the criteria. We can also confirm the consistency of the estimator, that is RASE_𝛽 decreases as the sample size increases.

Table 2 

Means (standard deviations) of RASE_𝛽.

The first step ahead prediction error (PE), which is defined as

PE=(𝑦_𝑛+1−𝑦_{𝑛+1∣𝑛})², (3.4) is also investigated. Table 3 shows the means and standard errors of PE for 𝜌=0,0.25,0.5,0.75 based on 500 simulations. The PE increases as the autoregressive coefficient increases, but the PE decreases as the sample size increases. From Table 3, we see that PLS + HT works better than the existing procedures and there is almost no difference in the PE depending on the information criteria. The models evaluated by BIC𝑚 perform well for large sample sizes.

Table 3 

Means (standard deviations) of first step ahead prediction errors.

The means and standard deviations of the number and deviation of basis functions are shown in Tables 4 and 5. The BIC𝑝 gives a smaller number of basis functions than the other information criteria. The models evaluated by BIC𝑝 also give smaller standard deviations of the number of basis functions. The models determined by BIC𝑝 tend to choose larger deviations of basis functions than those based on BIC𝑚 and GIC. The number of basis functions increases gradually as the sample size or 𝜌 increase. From Table 4, it appears that the number of basis functions does not depend on the sample size 𝑛. From Table 5, it also appears that the deviations of basis functions do not depend on the sample size 𝑛 and 𝜌.

Table 4 

Means (standard deviations) of the number of basis functions.

Table 5 

Means (standard deviations) of the deviations of basis functions.

We now compare the performance of our procedure with existing procedures in terms of the reduction of model complexity. Table 6 shows simulation results of the means and standard deviations of the number of parameters excluded (𝛽=0 or 𝑤=0) by the proposed procedure. The results indicate that the proposed procedure reduces model complexity. From Table 6, It appears that the models determined by BIC𝑝 tend to exclude fewer parameters and give smaller standard deviations for the number of parameters excluded. This is due to the selection of a smaller number of basis functions compared to the selection based on the other criteria (see Table 4). There is almost no dependence of the number of excluded parameters on 𝜌. The models evaluated by BIC𝑝 give a larger number of excluded parameters as the sample size increases. On the other hand, the models evaluated by BIC𝑚 or GIC give a smaller number of excluded parameters as the sample size increases.

Table 6 

Means (standard deviations) of the number of parameters excluded.

Table 7 shows the means and standard deviations of the number of basis functions excluded as 𝑤=0 by the proposed procedure. From Table 7 it appears that the models evaluated by BIC𝑝 tend to exclude fewer basis functions than those based on GIC and BIC. Again this is due to the selection of a smaller number of basis functions (see Table 4). The models determined by BIC𝑝 also give smaller standard deviations of the number of basis functions than the other criteria. There is almost no dependence of the number of basis functions on 𝜌.

Table 7 

Means (stnadard deviations) of the number of basis functions excluded.

Table 8 shows the means and standard deviations of the number of basis functions excluded as 𝛽=0 by the proposed procedure. The number of 𝛽 which values really 0 was five. From Table 8 we see that the proposed procedure gives nearly five. The models determined by BIC𝑝 give results more close to five and smaller standard deviations of the number of basis functions than the other criteria. The number of basis functions approaches five as the sample size increases. The standard deviations of the number of basis functions excluded decrease as 𝜌 increases. These results indicate that the proposed procedure reduces model complexity.

Table 8 

Means (standard deviations) of the number of explanatory variables excluded.

Table 9 shows the percentage of times that various 𝛽_𝑖 were estimated as being zero. As for the parameters 𝛽_𝑗≠0,𝑗=1,2,5, these parameters were not estimated zero for every simulations, we omit to show the corresponding results on Table 9. The results indicate that the proposed procedure excludes insignificant variables and selects significant variables. It can be seen that the proposed procedure gives a better performance as the sample size increases and that BIC𝑝 is superior to the other criteria.

Table 9 

Percentage of times 𝛽_𝑖 is estimated as zero.

4. Real Data Analysis

In this section we present the consequence of analyzing the real-time series data using proposed procedure. We use two data in this study; the data about the spirit consumption data in United Kingdom and the association between fertility and female employment in Japan.

4.1. The Spirit Consumption Data in the United Kingdom

We now illustrate our theory through an application to spirit consumption data for the United Kingdom from 1870 to 1938. The data-set can be found in Fuller [26, page 523]. In this data-set, the dependent variable 𝑦_𝑖 is the logarithm of the annual per capita consumption of spirits. The explanatory variables 𝑥_𝑖1 and 𝑥_𝑖2 are the logarithms of per capita income and the price of spirits, respectively, and 𝑥_𝑖3=𝑥_𝑖1𝑥_𝑖2. Figure 2 shows that there is a change-point at the start of the First World War (1914). Therefore, we prepare a variable 𝑧: 𝑧=0 from 1870 to 1914 and 𝑧=1 from 1915 to 1933. From this we derive another three explanatory variables: 𝑥_𝑖4=𝑥_𝑖1𝑧, 𝑥_𝑖5=𝑥_𝑖2𝑧, and 𝑥_𝑖6=𝑥_𝑖3𝑧. We consider the semiparametric model:

(a)

(b)

(c)

Figure 2 

Application data-set: (a): 𝑦𝑖=log(the annual per capita consumption of spirits); (b): 𝑥𝑖1=log(per capita income); (c): 𝑥𝑖2=log(price of spirits).

𝑦_𝑖=𝛼(𝑡_𝑖)+𝜷^′𝐱_𝑖+𝜀_𝑖,𝑖=1,…,69. (4.1)

We computed the basis function estimate for 𝛼 using the existing procedure (PLS) and the proposed procedure (PLS + HT) with BIC𝑝. The resulting estimates and standard errors (SEs) of 𝜷 are given in Table 10. The selected number of basis function is seven with one excluded basis function and the spread parameter 𝑠 is estimated as 0.12. Therefore, we obtain the model

Table 10 

Estimated Coefficients for Model 4.1.

𝑦_𝑖=𝛼(𝑡_𝑖)+2.395𝑥_𝑖4−2.411𝑥_𝑖5,𝑖=1,…,69. (4.2) The results indicate that the proposed procedure excludes some variable and reduces model complexity. Table 10 shows that PLS + HT selects only 𝛽₄ and 𝛽₆. That indicates possible interactions between consumption and income and between consumption and income×price after 1915. Consumption increases as income increases; however, as income increases and the price also increases, consumption decreases. We plot the estimated trend curve, residuals and autocorrelations functions in Figures 3 to 5. The residual mean square is 1.7×10⁻⁴.

Figure 3 

Plots of estimated curves. The solid line represents 𝑦; the dotted lines are the estimates of 𝑦; the dashed lines are the estimated curve 𝛼.

Figure 4 

Plot of residuals.

Figure 5 

ACF plot of residuals.

Figure 6 

Plots of standardized total fertility rate and female labor force participation rate for women aged 15 to 49 in Japan, 1968–2007. The solid line represents standardized TRF; the dotted lines are standardized FLP.

You and Chen [10] used the following semiparametric partially linear regression model:

𝑦_𝑖=𝛼(𝑡_𝑖)+𝛽₁𝑥_𝑖1+𝛽₂𝑥_𝑖2+𝜀_𝑖,𝑖=1,…,69. (4.3) The semiparametric least-square (SLS) regression gives 𝑦_𝑖=𝛼(𝑡_𝑖)+0.65𝑥_𝑖1−0.95𝑥_𝑖2. The residual mean square is 2.2×10⁻⁴, which is more than that of our SGLSE fit. For a fair comparison, we use model (4.3) to revise You and Chen's estimation. Our semiparametric generalized least-square gives 𝑦_𝑖=𝛼(𝑡_𝑖)−0.71𝑥_𝑖2. The result indicates that 𝑥_𝑖1 is insignificant in model (4.3).

4.2. The Association between Fertility and Female Employment in Japan

Recent literature finds that in OECD countries the cross-country correlation between the total fertility rate and the female labor force participation rate, which until the beginning of the 1980s had a negative value, has since acquired a positive value. See for example, Brewster and Rindfuss [27], Ahn and Mira [28], and Engelhardt et al. [29]. This result is often interpreted as evidence for a changing sign in the time series association between fertility and female employment within OECD countries.

However, OECD countries, including Japan, have different cultural backgrounds. We investigate whether or not the relation between the total fertility rate (TFR) and the female labor force participation rate (FLP) has changed in Japan from a negative value to a positive value. This application challenges previous findings and could be good news for policy makers, as a positive relationship implies that a rising FLP is associated with an increasing TFR.

Usually, FLP contains all women aged 15 to 64. However, TFR is a combination of fertility rates for ages 15–49, so we use the FLP of women aged 15 to 49 instead of women aged 15 to 64. We take the TFR from 1968 to 2007 in Japan. The estimation is a semiparametric regression of log(TFR_𝑖) on log(FLP_𝑖). As the law of the Equal Employment Act came into force in 1985, we use the interaction variables “dummy for 1968–1984 ×log(TFR)” (𝑥_𝑖2) and for 1985–2007 (𝑥_𝑖3). We also use dummy variables for 1990–1999 and 2000–2007 (𝑥_𝑖4, 𝑥_𝑖5) and consider the semiparametric model

log(TFR)_𝑖=𝛼(𝑡_𝑖)+𝜷^′𝐱_𝑖+𝜀_𝑖,𝑖=1,…,40. (4.4)

We applied the existing procedure (PLS) and proposed procedure (PLS + HT) with BIC_𝑝. The resulting estimates and standard errors (SE) of 𝜷 are given in Table 11. Therefore, we obtain the model

Table 11 

Estimated coefficients for Model 4.4.

𝑦_𝑖=𝛼(𝑡_𝑖)−0.27𝑥_𝑖1−0.20𝑥_𝑖2,𝑖=1,…,40. (4.5) The residual mean square of PLS + HT is 2.24×10⁻² and that of PLS is 2.47×10⁻². The selected number of basis functions is six with one excluded basis function and the spread parameter 𝑠 is estimated as 0.30. Table 11 shows that PLS + HT selects only log(FLP_𝑖) 1968–1984 and 1985–2007. That indicates a negative correlation between TFR and FLP for 1968–2007, especially for 1968–1984, which means TFR decreases as FLP increases. We could not see the positive association in 80s which has been reported in recent studies, see, for example, Brewster and Rindfuss [27], Ahn and Mira [28], and Engelhardt et al. [29]. We plot the estimated trend curve, residuals and autocorrelation functions in Figures 7 to 9.

Figure 7 

Plots of estimated curves, the solid line represents 𝑦; the dotted lines are the estimated curves of 𝑦; the dashed lines are the estimated curves of 𝛼.

Figure 8 

Plot of residuals.

Figure 9 

ACF plot of residuals.

5. Concluding Remarks

In this article we have proposed variable and model selection procedures for the semiparametric time series regression. We used the basis functions to fit the trend component. An algorithm of estimation procedures is provided and the asymptotic properties are investigated. From the numerical simulations, we have confirmed that the models determined by the proposed procedure are superior to those based on the existing procedure. They reduce the complexity of models and give good fitting by excluding basis functions and nuisance explanatory variables.

The development here is limited to the case with Gaussian stationary errors, but it seems likely our approach can be extended to the case with non-Gaussian long-range dependent errors, along with the lines explored in recent work by Aneiros-Pérez et al. [30]. However, the efficient estimation for regression parameter is an open question in case of long-range dependence. This is a question we hope to address in future work. We also plan to explore the question of whether the proposed techniques can be extended to the cointegrating regression models with an autoregressive distributed lag framework.

Appendix

Proofs

In this appendix we give the proofs of the theorems in Section 2. We use ‖𝑥‖ to denote the Euclidian norm of 𝑥.

Let 𝐚_𝜏,𝑛=(𝑎_1,𝑛,…,𝑎_𝜏,𝑛)^′ be the infeasible estimator for 𝐚_𝜏=(𝑎₁,…,𝑎_𝜏)^′ constructed using OLS methods. That is 𝐚_𝜏,𝑛=(𝑎_1,𝑛,𝑎_2,𝑛,…,𝑎_𝜏,𝑛)^′=(𝐄^′_𝜏𝐄_𝜏)⁻¹𝐄^′_𝜏𝜺, where 𝜺=(𝜀_𝜏+1,…,𝜀_𝑛)^′ and 𝐄_𝜏=[𝜀_𝑖,𝑗],𝑖=1,…,𝑛,𝑗=1,…,𝜏 with 𝜀_𝑖,𝑗=𝜀_{𝑖−𝑗−𝜏}. For ease of notation, we set 𝑎_𝑗,𝑛=𝑎_𝑗,𝑛=0 for 𝑗>𝜏, and 𝑎_0,𝑛=𝑎_0,𝑛=1. We write Γ(𝑘) for cov(𝜀₀,𝜀_𝑘). Then we can construct the infeasible estimate 𝐕 using 𝐚_𝜏,𝑛 and Γ(𝑘),𝑘=0,…,𝜏. The following lemma states that the estimators 𝜷 and 𝐰 given in Theorem 2.1 have asymptotically normal distributions.

Lemma .1. Under the assumptions of Theorem 2.1, one has √𝑛(𝜷−𝜷)→_𝐷𝑁(0,𝚺₁⁻¹√),𝑛(𝐰−𝐰)→_𝐷𝑁(0,𝚺₂⁻¹), (A.1)(A.2) where 𝚺₁⁻¹ and 𝚺₂⁻¹ are defined in Theorem 2.1.

Proof of Lemma .1. From model (2.6), 𝐲−𝐗𝜷−𝐁𝐰 can be written as 𝐲−𝐗𝜷−𝐁𝐰=𝐲−𝐁𝐰+(𝐲−𝐗𝜷)−(𝐲−𝐗𝜷)−𝐗𝜷=(𝐲−𝐗𝜷)+(𝐗−𝐗)𝜷−(𝐲−𝐲)−𝐁𝐰=(𝐲−𝐗𝜷)+𝐒(𝐲−𝐗𝜷)−𝐁𝐰=(𝐲−𝐗𝜷)−𝐁(𝐰−𝐰), (A.3) where 𝐲, 𝐗, and 𝐰 are given by 𝐲=(𝐈−𝐒)𝐲, 𝐗=(𝐈−𝐒)𝐗, and 𝐰=(𝐁^′𝐁+𝑛𝜉𝐊)⁻¹𝐁^′𝜺, respectively. Hence 𝐰 can be expressed without using 𝜷. Then the minimization function ℒ(𝜷,𝐰) in (2.17) can be written as 1ℒ(𝜷,𝐰)=2{(𝐲−𝐗𝜷)^′𝐕⁻¹(𝐲−𝐗𝜷)−2(𝐰−𝐰)^′𝐁^′𝐕⁻¹(𝐲−𝐗𝜷)+(𝐰−𝐰)^′𝐁^′𝐕⁻¹𝐁(𝐰−𝐰)}+𝛼𝐰^′𝐊𝐰≡𝐼₁(𝜷)+𝐼₂(𝜷,𝐰)+𝐼₃(𝐰)+𝐼₄(𝐰). (A.4) First we consider asymptotic normality for 𝐰, using the model 𝐲=𝐗𝜷₀+𝐁𝐰₀+𝜺. (A.5) The estimators 𝐰 minimize the function ℒ(𝜷,𝐰), which yields 𝜕ℒ(𝜷,𝐰)𝜕(𝐰)=𝐼^′₂(𝜷,𝐰)+𝐼^′₃(𝐰)+𝐼^′₄(𝐰)=−𝐁^′𝐕⁻¹(𝐲−𝐗𝜷)+𝐁^′𝐕⁻¹𝐁(𝐰−𝐰₀)+𝐁^′𝐕⁻¹𝐁(𝐰₀−𝐰)+2𝑛𝜉𝐊(𝐰−𝐰₀)+𝑛𝜉𝐊𝐰₀. (A.6) Then the minimization of this quadratic function is given by 𝐰−𝐰₀=(𝐁^′𝐕⁻¹𝐁+𝑛𝜉𝐊)⁻¹𝐁^′𝐕⁻¹{(𝐲−𝐗𝜷)−𝐁(𝐰−𝐰_𝟎)−𝑛𝜉𝐕𝐁⁻¹𝐊𝐰₀}=(𝐁^′𝐕⁻¹𝐁+𝑛𝜉𝐊)⁻¹𝐁^′𝐕⁻¹(𝐲−𝐗𝜷)+(𝐁^′𝐕⁻¹𝐁+𝑛𝜉𝐊)⁻¹𝐁^′𝐕⁻¹𝐁(𝐰₀−𝐰)−𝑛𝜉(𝐁^′𝐕⁻¹𝐁+𝑛𝜉𝐊)⁻¹𝐁^′𝐕⁻¹𝐕𝐁⁻¹𝐊𝐰₀≡𝐴₁+𝐴₂+𝐴₃. (A.7) We now deal with 𝐴₁, 𝐴₂, and 𝐴₃. First we evaluate 𝐴₁. From the expansion (𝐀+𝑎𝐁)⁻¹=𝐀⁻¹−𝑎𝐀⁻¹𝐁𝐀⁻¹+𝑂(𝑎²), we can see that 𝐴₁=(𝐁^′𝐕⁻¹𝐁+𝑛𝜉𝐊)⁻¹𝐁^′𝐕⁻¹𝜺𝐁=(^′𝐕⁻¹𝐁𝑛+𝜉𝐊)⁻¹𝐁^′𝐕⁻¹𝜺𝑛𝐁={(^′𝐕⁻¹𝐁𝑛)⁻¹𝐊𝐰₀𝐁+𝜉(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁𝐊(^′𝐕⁻¹𝐁𝑛)⁻¹𝐊𝐰₀+𝑂(𝜉²𝐁)}^′𝐕⁻¹𝑛𝜺𝐁=(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁^′𝐕⁻¹𝑛𝐁𝜺−𝜉(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁𝐊(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁^′𝐕⁻¹𝑛𝜺+𝑂(𝜉²)𝐁′𝐕⁻¹𝑛𝐁𝜺,=(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁^′𝐕⁻¹𝑛𝜺+𝑂(𝜉)+𝑂(𝜉²). (A.8) Similarly, we obtain 𝐴₂=(𝐁^′𝐕⁻¹𝐁+𝑛𝜉𝐊)⁻¹𝐁^′𝐕⁻¹𝐁(𝐰₀−𝐁𝐰)=(^′𝐕⁻¹𝐁𝑛+𝜉𝐊)⁻¹𝐁^′𝐕⁻¹𝐁𝑛(𝐰₀−𝐁𝐰)={(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁−𝜉(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁𝐊(^′𝐕⁻¹𝐁𝑛)⁻¹+𝑂(𝜉²𝐁)}×(^′𝐕⁻¹𝐁𝑛)(𝐰₀−𝐰)=(𝐰₀−𝐁𝐰)−𝜉(^′𝐕⁻¹𝐁𝑛)⁻¹𝐊(𝐰₀−𝐰)+𝑂(𝜉²𝐁)(^′𝐕⁻¹𝐁𝑛)(𝐰₀−𝐰). (A.9) Finally, we can evaluate 𝐴₃ as follows: 𝐴₃=−(𝐁^′𝐕⁻¹𝐁+𝑛𝜉𝐊)⁻¹𝐁^′𝐕⁻¹𝐁⁻¹𝑛𝜉𝐊𝐰₀𝐁=(^′𝐕⁻¹𝐁𝑛+𝜉𝐊)⁻¹𝜉𝐊𝐰₀𝐁={(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁−𝜉(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁𝐊(^′𝐕⁻¹𝐁𝑛)⁻¹+𝑂(𝜉²)}𝜉𝐊𝐰₀𝐁=−𝜉(^′𝐕⁻¹𝐁𝑛)⁻¹𝐊𝐰₀+𝜉²(𝐁^′𝐕⁻¹𝐁𝑛𝐁)𝐊(^′𝐕⁻¹𝐁𝑛)𝐊𝐰₀+𝑂(𝜉³)𝐊𝐰₀. (A.10)(1)(2)(3)(4) We can also observe that the weighted least-square estimates 𝐰 have a normal distribution. Hence 𝐰−𝐰₀=𝑂_𝑝(𝑛^−1/2). (A.11) If 𝜉=𝑂(𝑛^𝜂) and 𝜂<−1/2, then 𝐴₁, 𝐴₂, and 𝐴₃ become 𝐴₁𝐁=(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁^′𝐕⁻¹𝑛𝜺+𝑂(𝜉)+𝑂(𝜉²𝐴),₂=(𝐰₀−𝐰)+𝑂(𝜉)×𝑂_𝑝(𝐰₀−𝐰)+𝑂(𝜉²)×𝑂_𝑝(𝐰₀−𝐰)=𝑂_𝑝(𝑛^−1/2), (A.12) and 𝐴₃=𝑂(𝜉)+𝑂(𝜉²)+𝑂(𝜉³)=𝑂(𝜉). Therefore, (A.7) can be written as 𝐰−𝐰₀𝐁=(^′𝐕⁻¹𝐁𝑛)⁻¹𝐁^′𝐕⁻¹𝑛𝜺+𝑂_𝑝(𝑛^−1/2). (A.13) By the law of large numbers and the central limit theorem, √𝑛(𝐰−𝐰₀)→_𝐷𝑁(0,𝚺₂⁻¹). (A.14) Next we deal with the estimators 𝜷. These minimize the function ℒ(𝜷,𝐰), which yields 𝜕ℒ(𝜷,𝐰)𝜕𝜷=𝐼^′₁(𝜷)+𝐼^′₂(𝜷,𝐰)=𝐗^′𝐕⁻¹𝜺+𝐗^′𝐕⁻¹𝐗(𝜷−𝜷₀)+𝐗^′𝐕⁻¹𝐁(𝐰−𝐰). (A.15) The minimization of this quadratic function is given by 𝜷=𝜷₀+(𝐗^′𝐕⁻¹𝐗)⁻¹{𝐗^′𝐕⁻¹𝜺+𝐗^′𝐕⁻¹𝐁(𝐰−𝐰)}=𝜷₀+(𝐗^′𝐕⁻¹𝐗)⁻¹𝐗^′𝐕⁻¹{𝜺+𝐁(𝐰−𝐰)}. (A.16) If we substitute 𝐰 for its estimator 𝐰₀, from (A.14) and (A.11), we have 𝜷=𝜷₀+(𝐗^′𝐕⁻¹𝐗)⁻¹{𝐗^′𝐕⁻¹𝜺+𝐗^′𝐕⁻¹𝐁(𝐰₀−𝐰)}=𝜷₀+(𝐗^′𝐕⁻¹𝐗)⁻¹𝐗^′𝐕⁻¹𝜺+𝑂_𝑝(𝑛^−1/2). (A.17) Similarly, by the law of large numbers and the central limit theorem, √𝑛(𝜷−𝜷₀)→_𝐷𝑁(0,𝚺₁⁻¹). (A.18) This completes the proof of the lemma.

Proof of Theorem 2.1. Let the estimator 𝐚_𝜏,𝑛=(𝑎_1,𝑛,…𝑎_𝜏,𝑛)^′ be the ordinary least-square estimate applied to model (2.18). For the ease of notation, we set 𝑎_𝑗,𝑛=𝑎_𝑗,𝑛=0 for 𝑗>𝜏 and 𝑎_0,𝑛=𝑎_0,𝑛=1. Then we write 𝑒̂_𝑖,𝑛=𝑒_𝑖−𝑆_𝑖,𝑛+𝑅_𝑖,𝑛+𝑄_𝑖,𝑛, (A.19) where 𝑆_𝑖,𝑛=^∞_𝑗=0𝑎_𝑗{(𝜷−𝜷)^′𝐱_𝑖−𝑗+(𝐰−𝐰)^′𝜙_𝑖−𝑗𝑅},_𝑖,𝑛=^𝜏_𝑗=0(𝑎_𝑗,𝑛−𝑎_𝑗,𝑛)(𝑦_𝑖−𝑗−𝜷^′𝐱_𝑖−𝑗−𝐰^′𝜙_𝑖−𝑗𝑄),_𝑖,𝑛=^∞_𝑗=0(𝑎_𝑗,𝑛−𝑎_𝑗)(𝑦_𝑖−𝑗−𝜷^′𝐱_𝑖−𝑗−𝐰^′𝜙_𝑖−𝑗). (A.20) From assumptions (A.1), (A.2), and Lemma .1 we can see that under the assumptions about 𝜏 and by the Caucy-Schwarz inequality |𝑆_𝑖,𝑛|≤^∞_𝑗=0|𝑎_𝑗||(𝜷−𝜷)^′𝐱_𝑖−𝑗+(𝐰−𝐰)^′𝜙_𝑖−𝑗|≤‖𝜷−𝜷‖‖^∞_𝑗=0𝑎_𝑗𝐱_𝑖−𝑗‖+‖𝐰−𝐰‖‖^∞_𝑗=0𝑎_𝑗𝜙_𝑖−𝑗‖=𝑂_𝑝(𝑛^−1/2). (A.21) Next we evaluate 𝑅_𝑖,𝑛. In An et al. [31, proof of Theorem  5]: it is shown that, under the assumptions about 𝜏(𝑛), ^𝜏_𝑗=0(𝑎_𝑗,𝑛−𝑎_𝑗,𝑛)²=𝑜((log(𝑛)𝑛)^1/2). (A.22) Thus, by the Cauchy-Schwarz inequality |𝑅_𝑖,𝑛|≤(^𝜏_𝑗=0(𝑎_𝑗,𝑛−𝑎_𝑗,𝑛)²)^1/2(^𝜏_𝑗=0(𝑦_𝑖−𝑗−𝜷^′𝐱_𝑖−𝑗−𝐰^′𝜙_𝑖−𝑗)²)^1/2, (A.23) which yields 𝑅_𝑖,𝑛=𝑜((log(𝑛)/𝑛)^1/4))𝑂_𝑝(𝜏^1/2)=𝑜_𝑝(1). Finally, we evaluate 𝑄_𝑖,𝑛. By the extended Baxter inequality from Bühlmann [32, proof of Theorem  3.1], we have ^∞_𝑗=0|𝑎_𝑗,𝑛−𝑎_𝑗|≤𝐶^∞_𝑗=𝜏+1|𝑎_𝑗|. (A.24) Notice that 𝑦_𝑖−𝑗−𝜷^′𝐱_𝑖−𝑗−𝐰^′𝝓_𝑖−𝑗=𝑒_𝑖,𝑛. Since 𝑒_𝑖,𝑛 is a stationary and invertible process whose linear process coefficients satisfy the given summability assumptions, we have for some 𝑀>0, |𝑄_𝑖,𝑛|≤𝑀^∞_𝑗=0|𝑎_𝑗,𝑛−𝑎_𝑗|≤𝑀^∞_𝑗=𝜏+1|𝑎_𝑗|=𝑜_𝑝(1). (A.25) From the above decomposition and evaluation, we can see that 𝐲−𝐗𝜷−𝐁𝐰=𝐲𝜏−𝐗𝜏𝜷−𝐁𝜏𝐰+𝑜_𝑝(1). (A.26) Therefore, in order to prove the second equation in Theorem 2.1 we just need to show 1𝑛(𝐗^′𝐕_𝜏⁻¹𝐗−𝐗^′𝐕_𝜏⁻¹𝐗)=𝑂_𝑝𝜏((𝑛)^1/21),𝑛(𝐁^′𝐕_𝜏⁻¹𝐁−𝐁^′𝐕_𝜏⁻¹𝐁)=𝑂_𝑝𝜏((𝑛)^1/21),√𝑛𝐗^′(𝐕_𝜏⁻¹−𝐕_𝜏⁻¹)𝜺=𝑂_𝑝𝜏((𝑛)^1/21),√𝑛𝐁^′(𝐕_𝜏⁻¹−𝐕_𝜏⁻¹)𝜺=𝑂_𝑝𝜏((𝑛)^1/2). (A.27) To see the above results are true, let 𝑦_𝜏,𝑖 be the 𝑖th element 𝐲_𝜏 of model (2.20). We have for 𝐓_𝜏,𝑖 (the 𝑖th row of 𝐓_𝜏), 𝐗_𝜏,𝑖 (the 𝑖th column of 𝐗_𝜏), and 𝐁_𝜏,𝑖 (the 𝑖th column of 𝐁_𝜏) 𝑒_𝑖̂=𝐓_𝜏,𝑖𝜺=𝐞_𝑖+^𝜏_𝑗=1(𝑎_𝑗−𝑎_𝑗)𝜀_𝑖−𝑗,𝑋𝜏_,𝑖𝑗=𝐓_𝜏,𝑗⋅𝐗_𝜏,⋅𝑖=𝑋_{𝜏,𝑖𝑗}+^𝜏_𝑗=1𝑎_𝑗𝑋_{𝑖−𝑗,𝑖}+^𝜏_𝑗=1(𝑎_𝑗−𝑎_𝑗)𝑋_{𝑖−𝑗,𝑖}≡𝑋_{𝜏,𝑖𝑗}+^𝜏_𝑗=1(𝑎_𝑗−𝑎_𝑗)𝑋_{𝑖−𝑗,𝑖},𝐵𝜏_,𝑖𝑗=𝐓_𝜏,𝑗⋅𝐁_𝜏,⋅𝑖=𝐵_{𝜏,𝑖𝑗}+^𝜏_𝑗=1𝑎_𝑗𝐵_{𝑖−𝑗,𝑖}+^𝜏_𝑗=1(𝑎_𝑗−𝑎_𝑗)𝐵_{𝑖−𝑗,𝑖}≡𝐵_{𝜏,𝑖𝑗}+^𝜏_𝑗=1(𝑎_𝑗−𝑎_𝑗)𝐵_{𝑖−𝑗,𝑖} (A.28) for 𝑖=𝜏+1,𝜏+2,…,𝑛, with similar expressions holding for 𝑖=1,2,…,𝜏. By (A.26) and the fact that ‖𝐚_𝜏−𝐚‖=𝑂_𝑝((𝜏/𝑛)^1/2) (see Xiao et al. [22]), it follows that 1√𝑛^𝑛_𝑖=1𝑒_𝑖̂𝑋𝜏_,𝑖𝑗=1√𝑛^𝑛_𝑖=1𝑒_𝑖𝑋_{𝜏,𝑖𝑗}+𝑂_𝑝𝜏((𝑛)^1/21),𝑛^𝑛_𝑖=1𝑋𝜏_,𝑖𝑗𝑋𝜏_,𝑖𝑘=1𝑛^𝑛_𝑖=1𝑋_{𝜏,𝑖𝑗}𝑋_{𝜏,𝑖𝑘}+𝑂_𝑝𝜏((𝑛)^1/21),√𝑛^𝑛_𝑖=1𝑒_𝑖̂𝐵_{𝜏,𝑖𝑗}=1√𝑛^𝑛_𝑖=1𝑒_𝑖𝐵_{𝜏,𝑖𝑗}+𝑂_𝑝𝜏((𝑛)^1/21),𝑛^𝑛_𝑖=1𝐵_{𝜏,𝑖𝑗}𝐵_{𝜏,𝑖𝑘}=1𝑛^𝑛_𝑖=1𝐵_{𝜏,𝑖𝑗}𝐵_{𝜏,𝑖𝑘}+𝑂_𝑝𝜏((𝑛)^1/2). (A.29) Therefore, using the expansion (𝐀+𝑎𝐁)⁻¹=𝐀⁻¹−𝑎𝐀⁻¹𝐁𝐀⁻¹+𝑂(𝑎²) and from (A.17), (A.13) and (A.27), we have √𝑛(𝜷_SGLSE−𝜷₀𝐗)=(^′𝐕⁻¹𝐗𝑛)⁻¹(1√𝑛{𝐗^′𝐕⁻¹𝜺+𝑂_𝑝(𝑛^−1/2𝐗)})=(^′𝐕⁻¹𝐗𝑛+𝑂_𝑝𝜏((𝑛)^1/2))⁻¹𝐗×(^′𝐕⁻¹√𝑛𝜺+𝑂_𝑝𝜏((𝑛)^1/2)+𝑂_𝑝(𝑛⁻¹=√))𝑛(𝜷−𝜷₀)+𝑂_𝑝𝜏((𝑛)^1/2√),𝑛(𝐰_SGLSE−𝐰₀𝐁)=(^′𝐕⁻¹𝐁𝑛)⁻¹{1√𝑛(𝐁^′𝐕⁻¹𝜺+𝑂_𝑝(𝑛^−1/2𝐁))}=(^′𝐕⁻¹𝐁𝑛+𝑂_𝑝𝜏((𝑛)^1/2))⁻¹(𝐁^′𝐕⁻¹√𝑛𝜺+𝑂_𝑝𝜏((𝑛)^−1/2)+𝑂_𝑝(𝑛⁻¹=√))𝑛(𝜷−𝜷₀)+𝑂_𝑝𝜏((𝑛)^1/2). (A.30) This completes the proof of Theorem 2.1.

Proof of Theorem 2.2. We write 𝛼_𝑛=𝑛^−1/2+𝑎_𝑛. It is sufficient to show that, for any given 𝜁>0, there exist large constants 𝐶 such that 𝑃{inf_{‖𝐮‖=𝐶}𝒮(𝜽₀+𝛼_𝑛𝐮)≥𝒮(𝜽₀)}≤1−𝜁. (A.31) This implies, with probability at least 1−𝜁, that there exists a local minimizer in the balls {𝜽₀+𝛼_𝑛𝐮∶‖𝐮‖≤𝐶}. Define 𝐷_𝑛(𝐮)=𝒮(𝜽₀+𝛼_𝑛𝐮)−𝒮(𝜽₀). (A.32) Note that 𝑝_𝜆𝑗𝑛(0)=0 and that 𝑝_𝜆𝑗𝑛(|𝜃_𝑗|) is nonnegative, so 𝐷_𝑛(𝐮)≥𝑛⁻¹{𝑙(𝜽₀+𝛼_𝑛𝐮)−𝑙(𝜽₀+)}_{𝑑+𝑚𝑗=1}{𝑝_𝜆𝑗𝑛(|𝜃_𝑗0+𝛼_𝑛𝑢_𝑗|)−𝑝_𝜆𝑗𝑛(|𝜃_𝑗0|)}+𝜉(𝜽₀+𝛼_𝑛𝐮)^′𝐊(𝜽₀+𝛼_𝑛𝐮)−𝜉𝜽^′₀𝐊𝜽₀, (A.33) where 𝑙(𝜽) is the first term of (2.7) and 𝐊 is defined in (2.47). Now 12𝑛⁻¹{𝑙(𝜽₀+𝛼_𝑛=𝛼𝐮)−𝑙(𝜽)}_21𝑛2𝐮^′₁{𝐁^′𝐕𝐁𝑛+𝑜_𝑝(1)}𝐮₁−𝛼_1𝑛𝐮^′₁{𝐁^′𝐕𝑛𝜺+𝑜_𝑝(𝑛^−1/2+𝛼)}_22𝑛2𝐮^′₂{𝐗^′𝐕⁻¹𝐗𝑛+𝑜_𝑝(1)}𝐮₂−𝛼_2𝑛𝐮^′₂{𝐗^′𝐕⁻¹𝑛𝜺+𝑜_𝑝(𝑛^−1/2)}. (A.34) Note that 𝐁^′𝐕⁻¹𝐁/𝑛→𝚺₁, 𝐁^′𝐕⁻¹𝜺/𝑛→𝝃₁, 𝐗^′𝐕⁻¹𝐗/𝑛→Σ₂, and 𝐗^′𝐕⁻¹𝜺/𝑛→𝝃₂ are finite positive definite matrices in probability. So the first term in the right side of (A.34) is of order 𝑂_𝑝(𝐶²₁𝛼²_1𝑛), and the second term is of order 𝑂_𝑝(𝐶₁𝑛^−1/2𝛼_1𝑛)=𝑂_𝑝(𝐶𝛼²_1𝑛). Similarly, the third term of (A.34) is of order 𝑂_𝑝(𝐶²₂𝛼²_2𝑛) and the fourth term is order 𝑂_𝑝(𝐶₂𝑛^−1/2𝛼²_2𝑛). Furthermore, ^𝑚_𝑗=1{𝑝_{𝜆𝑗1𝑛}(|𝑤_𝑗0+𝛼_1𝑛𝑢_𝑗|)−𝑝_{𝜆𝑗1𝑛}(|𝑤_𝑗0|)},_𝑑𝑗=1{𝑝_{𝜆𝑗2𝑛}(|𝛽_𝑗0+𝛼_2𝑛𝑢_𝑗|)−𝑝_{𝜆𝑗2𝑛}(|𝛽_𝑗0|)}, (A.35)(A.36) are bounded by √𝑚𝛼_1𝑛𝑎_1𝑛‖𝐮‖+𝛼²_1𝑛𝑏_1𝑛‖𝐮‖²=𝐶𝛼²_𝑛(√𝑚+𝑏_1𝑛√𝐶),𝑑𝛼_2𝑛𝑎_2𝑛‖𝐮‖+𝛼²_𝑛𝑏_2𝑛‖𝐮‖²=𝐶𝛼²_2𝑛(√𝑑+𝑏_2𝑛𝐶) (A.37) by the Taylor expansion and the Cauchy-Schwarz inequality. As 𝑏_𝑛→0, the first term on the right side of (A.34) will dominate (A.35) and (A.36) as well as the second term on the right side of (A.34), by taking 𝐶 sufficiently large. Hence (A.31) holds for sufficiently large 𝐶. This completes the proof of the theorem.

Lemma .2. Under the conditions of Theorem 2.3, with probability tending 1, for any given 𝜷 and 𝐰, satisfying ‖𝜷₁−𝜷₁₀‖=‖𝐰₁−𝐰₁₀‖=𝑂_𝑝(𝑛^−1/2) and any constant 𝐶, 𝒮{(𝜷^′_𝟏,𝟎^′)^′,(𝐰^′_𝟏,𝟎^′)}=min_{‖𝜷2‖≤𝐶1𝑛−1/2,‖𝐰2‖≤𝐶2𝑛−1/2,}𝒮{(𝜷^′_𝟏,𝜷^′_𝟐)^′,(𝐰^′_𝟏,𝐰^′_𝟐)^′}. (A.38)

Proof. We show that with probability tending to 1, as 𝑛→∞, for any 𝜷₁ and 𝐰₁ satisfying ‖𝜷₁−𝜷₁₀‖=‖𝐰₁−𝐰₁₀‖=𝑂_𝑝(𝑛^−1/2), ‖𝜷₂‖≤𝐶₁𝑛^−1/2, and ‖𝐰₂‖≤𝐶₂𝑛^−1/2, 𝜕𝑙(𝜷,𝐰)/𝜕𝛽_𝑗 and 𝛽_𝑗 have the same signs for 𝛽_𝑗∈(−𝐶₁𝑛^−1/2,𝐶₁𝑛^−1/2), for 𝑗=𝑆₁+1,…,𝑑. Also 𝜕𝑙(𝜷,𝐰)/𝜕𝑤_𝑗 and 𝑤_𝑗 have the same signs for 𝑤_𝑗∈(−𝐶₂𝑛^−1/2,𝐶₂𝑛^−1/2), for 𝑗=𝑆₂+1,…,𝑚. Thus the minimization is attained at 𝜷₂=𝐰₂=0.
For 𝛽_𝑗≠0 and 𝑗=𝑆₁+1,…,𝑑,
𝜕𝒮(𝜷)𝜕𝛽_𝑗=𝑙^′_𝑗(𝛽)+𝑛𝑝_{𝜆𝑗2𝑛}(|𝛽_𝑗|)sgn(𝛽_𝑗), (A.39) where 𝑙^′_𝑗(𝜷)=𝜕𝑙(𝜷)/𝜕𝛽_𝑗. By the proof of Theorem 2.1, 𝑙^′_𝑗(𝜷)=−𝑛{𝜉̂_2𝑗−(𝜷−𝜷₀)^′𝚺_1𝑗+𝑂_𝑝(𝑛^−1/2)}, (A.40) where 𝜉_2𝑗 is the 𝑗th component of 𝝃_2𝑛 and 𝚺_1𝑗 is the 𝑗th column of 𝚺₁. From ‖𝜷−𝜷₀‖=𝑂_𝑝(𝑛^−1/2), 𝑛⁻¹𝑙^′_𝑗(𝜷) is of order 𝑂_𝑝(𝑛^−1/2). Therefore, 𝜕𝒮(𝜷)𝜕𝛽_𝑗=𝑛𝜆_𝑗2𝑛{𝜆_𝑗1𝑛𝑝^′_{𝜆𝑗1𝑛}(|𝛽_𝑗|)sgn(𝛽_𝑗)+𝑂_𝑝(𝑛^−1/2/𝜆_1𝑛)}. (A.41) Because liminf_𝑛→∞liminf_{𝛽𝑗→0+}𝜆⁻¹_𝑗1𝑛𝑝^′_{𝜆𝑗1𝑛}(|𝛽_𝑗|)>0 and 𝑛^−1/2𝜆_𝑗1𝑛→0, the sign of the derivative is completely determined by that of 𝛽_𝑗.
For 𝑤_𝑗≠0 and 𝑗=𝑆₁+1,…,𝑚,
𝜕𝒮(𝐰)𝜕𝑤_𝑗=𝑙^′_𝑗(𝑤)+𝑛𝑝_{𝜆𝑗1𝑛}(|𝑤_𝑗|)sgn(𝑤_𝑗)+2𝑛𝜉𝐊𝐰, (A.42) where 𝑙^′_𝑗(𝐰)=𝜕𝑙(𝐰)/𝜕𝑤_𝑗. By the proof of Theorem 2.1, 𝑙^′_𝑗(𝐰)=−𝑛{𝜉̂_1𝑗−(𝐰−𝐰₀)^′𝚺_2𝑗+𝑂_𝑝(𝑛^−1/2)}, (A.43) where 𝜉_1𝑗 is the 𝑗th component of 𝝃_1𝑛 and 𝚺_2𝑗 is the 𝑗th column of 𝚺₂. From ‖𝐰−𝐰₀‖=𝑂_𝑝(𝑛^−1/2), 𝑛⁻¹𝑙^′_𝑗(𝐰) is of order 𝑂_𝑝(𝑛^−1/2). Therefore, 𝜕𝒮(𝐰)𝜕𝑤_𝑗=𝑛𝜆_𝑗2𝑛{𝜆_𝑗2𝑛𝑝^′_{𝜆𝑗2𝑛}(|𝑤_𝑗|)sgn(𝑤_𝑗)+𝑂_𝑝(𝑛^−1/2/𝜆_2𝑛)}. (A.44) Because liminf_𝑛→∞liminf_{𝑤𝑗→0+}𝜆⁻¹_𝑗2𝑛𝑝^′_{𝜆𝑗2𝑛}(|𝑤_𝑗|)>0 and 𝑛^−1/2𝜆_𝑗2𝑛→0, 𝑛^−1/2𝜆_𝑗2𝑛→0, the sign of the derivative is completely determined by that of 𝑤_𝑗. This completes the proof.

Proof of Theorem 2.3. Part (a) follows directly from follows by Lemma .2. Now we prove part (b). Using an argument similar to the proof of Theorem 2.1, it can be shown that there exist a 𝐰₁ and 𝜷₁ in Theorem 2.3 that are a root-𝑛 consistent local minimizer of 𝒮{(𝐰^′_𝟏,𝟎^′)^′} and 𝒮{(𝜷^′_𝟏,𝟎^′)^′}, satisfying the penalized least-square equations: 𝜕𝒮(𝐰^′_𝟏,𝟎^′)^′𝜕𝐰₁=𝟎,𝜕𝒮(𝜷^′_𝟏,𝟎^′)^′𝜕𝜷₁=𝟎. (A.45) Following the proof of Theorem 2.1, we have 𝜕𝓢(𝐰^′_𝟏,𝟎^′)^′𝜕𝑤_𝑗=𝑛[𝝃₂₍₁₎+𝑂_𝑝(𝑛^−1/2)+𝑛𝜉𝐊𝐰+{𝚺₂₍₁₎+𝑂_𝑝(1)}(𝐰₁−𝐰₁₀)]+𝑛[𝐛_𝑛+𝚺_𝜆1(𝐰){1+𝑂_𝑝(1)}(𝐰₁−𝐰₁₀)],𝜕𝓢(𝜷^′_𝟏,𝟎^′)^′𝜕𝛽_𝑗=𝑛[𝝃₁₍₁₎+𝑂_𝑝(𝑛^−1/2)+{𝚺₁₍₁₎+𝑂_𝑝(1)}(𝜷₁−𝜷₁₀)]+𝑛[𝐛_𝑛+𝚺_𝟐+𝚺_𝜆2(𝜷){1+𝑂_𝑝(1)}(𝜷₁−𝜷₁₀)] (A.46) where 𝜉̂₁₍₁₎ and 𝜉̂₂₍₁₎ consist of the first 𝑆_𝑗, 𝑗=1,2,…,𝑆₁ and 𝑗=1,2,…,𝑆₂ components of 𝝃₁ and 𝝃₂ respectively. Also 𝚺₁₍₁₎ and 𝚺₂₍₁₎ consist of the first 𝑆_𝑗, 𝑗=1,2,…,𝑆₁ and 𝑗=1,2,…,𝑆₂ rows and columns of 𝚺₁ and 𝚺₂, respectively.
Therefore, similar to the proof of Theorem 2.1 and by Slutsky's theorem, it follows that
√𝑛(𝐈_𝑆1+𝚺_𝜆1(𝜷))(𝜷₁−𝜷₁₀+(𝐈_𝑆1+𝚺_𝜆1(𝜷))⁻¹𝐛_𝛽)→𝑁_𝑆1(0,𝚺⁻¹₁₍₁₎√),𝑛(𝐈_𝑆2+𝚺_𝜆2(𝐰))(𝐰₁−𝐰₁₀+(𝐈_𝑆2+𝚺_𝜆2(𝐰)+𝜉𝐊)⁻¹𝐛_𝑤)→𝑁_𝑆2(0,𝚺⁻¹₂₍₁₎) (A.47) This completes the proof of Theorem 2.3.

Acknowledgments

The authors are grateful to two anonymous referees whose probing questions have led to a substantial improvement of the paper. This research was supported by the Norinchukin Bank and Nochu Information System endowed chair of Financial Engineering in the Department of Management Science, Tokyo University of Science.