# Model and Variable Selection Procedures for Semiparametric Time Series Regression

Authors: Risa Kato; Takayuki Shiohama
Journal: Journal of Probability and Statistics
 (2009)

---

## 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.


---

## Body

## 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 Section2.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 BICm proposed by Hastie and Tibshirani [16], and the BICp 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:(2.1)yi=α(ti)+β′xi+εi,i=1,…,n,
where yi is the response variable and xi is the d×1 covariate vector at time i, α(ti) is an unspecified baseline function of ti with ti=i/n, β is a vector of unknown regression coefficients, and εi is a Gaussian and zero mean covariance stationary process.We assume the following properties for the error termsεi and vectors of explanatory variables xi.(A.1)
It holds that{εi} is a linear process given by
(2.2)εi=∑j=0∞bjei-j,
where b0=1 and {ei} is an i.i.d. Gaussian random variable with E{ei}=0 and E{ei2}=σe2.(A.2)
The coefficientsbj satisfy the conditions that for all |z|<1, ∑j=0∞bjzj≠0 and ∑j=0∞j2|bj|<∞.We defineγ(k)=cov(εt,εt+k)=E{εtεt+k}.The assumptions on covariate variables are as follows.(B.1)
Alsoxi=(xi1,…,xid)′∈ℝd and {xij},j=1,…d, have mean zero and variance 1.The trend functionα(ti) is expressed as a linear combination of a set of m underlying basis functions :(2.3)α(ti)=∑k=1mwkϕk(ti)=w′ϕ(t),
where {ϕ(ti)=(ϕ1(ti),…,ϕm(ti))′} is an m-dimensional vector constructed from basis functions {ϕk(ti);k=1,…,m}, and w=(w1,…,wm)′ 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)ϕ(ti)=(ti,…,tip,(ti-κ1)+p,…,(ti-κk)+p)′,
where {κk}k=1,…,K are spline knots. This specification uses the so-called truncated power function basis. The choice of the number of knots K 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(2.5)ϕk(ti)=exp(-∥ti-μk∥22sk2),
where μk determines the location and sk2 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)y=Xβ+Bw+ε,
where X=(x1,…,xn)′, y=(y1,…,yn)′, B=(ϕ1,…,ϕn)′ with ϕi=(ϕ1(i/n),…,ϕm(i/n))′. The penalized least-square estimator is then a minimizer of the function(2.7)12(y-Xβ-Bw)′(y-Xβ-Bw)+nξw′Kw,
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  K is an appropriate positive semidefinite symmetric matrix. For example, if K satisfies w′Kw=∫01[α′′(u)]2du, we have the usual quadratic integral penalty (see, e.g., Green and Silverman [19]). By simple calculus, (2.7) is minimized when β and w satisfy the block matrix equation(2.8)(X′XX′BB′XB′B+nξK)(βw)=(X′B′)y.
This equation can be solved without any iteration (see, e.g., Green [20]). First, we find Bw̃=S(y-Xβ), where S=B(B′B+αK)-1B′ is usually called the smoothing matrix. Substituting Bw̃ into (2.6), we obtain(2.9)ỹ=X̃β+ε,
where ỹ=(I-S)y, X̃=(I-S)X ,and I is the identity matrix of order n. Applying least-square to the linear model (2.9), we obtain the semiparametric ordinary least-square estimator (SOLSE) result:(2.10)β̂SOLSE=(X̃′X̃)-1X̃′ỹ,(2.11)ŵSOLSE=(B′B+nξK)-1B′(y-Xβ̂SOLSE).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 {εi} in (2.6) are invertible. Let b(L)=∑j=1∞bjei-j, where L is the lag operator and a(L)=b(L)-1=a0-∑j=1∞ajLj with a0=1. Applying a(L) to the model (2.6) and rewriting the corresponding equation, we obtain the new model:(2.12)y̲=X̲β+B̲w+e,
where y̲=(y̲1,…,y̲n)′,X̲=(x̲1,…,x̲n)′,B̲=(ϕ1̲,…,ϕn̲)′ and e=(e1,…,en)′. Here
(2.13)y̲i=yi-∑j=1∞ajyi-j,ϕ̲i=ϕi-∑j=1∞ajϕi-j,x̲i=xi-∑j=1∞ajxi-j.
The regression errors in (2.12) are i.i.d. Because, in practice, the response variable yi is unknown, we use a reasonable approximation by y̲i based on the work by Xiao et al. [22] and Aneiros-Pérez and Vilar-Fernández [11].Under the usual regularity conditions the coefficientsaj decrease geometrically so, letting τ=τ(n) denote a truncation parameter, we may consider the truncated autoregression on εi:(2.14)ei=εi-∑j=1∞ajεi-j,
where ei are i.i.d. random variables with E(ei)=0. We make the following assumption about the truncation parameter.(C.1)
The truncation parameterτ satisfies τ(n)=clogn for some c>0.The expansion rate of the truncation parameter given in (C.1) is also for convenience. LetTτ be the n×n transformation matrix such that eτ=Tτε. Then the model (2.12) can be expressed as(2.15)Tτy=TτXβ+TτBw+Tτε,
where(2.16)Tτ=(δ110⋯0δ21-δ220⋯0⋮δτ1⋯-δττ-aτ⋯-a110-aτ⋯-a11⋮0⋯0-aτ⋯-a11)
with δ11=σe/γ(0),δ22=σe/(1-ρ2(1))γ(0),δ21=ρ(1)(σe/(1-ρ2(1))γ(0)),…. Here ρ(h)=γ(h)/ρ(0) denotes the lag h autocorrelation function of {εi}.Now our estimation problem for the semiparametric time series regression model can be expressed as the minimization of the function(2.17)ℒ(β,w)=12(y-Xβ-Bw)′V-1(y-Xβ-Bw)+αw′Kw,
where V-1=σe-2Tτ'Tτ and σe2=n-1∥Tτε∥2. Based on the work by Aneiros-Pérez and Vilar-Fernández [11], an estimator for Tτ is constructed as follows. We use the residuals ε̂=y-Xβ̂SOLSE-BŵSOLSE to construct an estimate of Tτ using the ordinary least square method applied to the model(2.18)ε̂i=a1ε̂i-1+⋯+aτε̂i-τ+residuali.
Define the estimate âτ=(â1,â2,…,âτ)′ of aτ=(a1,a2,…,aτ)′, where
(2.19)âτ=(Êτ′Êτ)-1Êτ′ε̂,
where ε̂=(ε̂τ+1,…,ε̂n) and Êτ is the (n-τ)×τ matrix of regressors with the typical element ε̂i-j. Then T̂τ is obtained from Tτ by replacing aj with âj, σe2 with σ̂e2, and so forth. Applying least-square to the linear model, we obtain
(2.20)T̂τy=T̂τXβ+T̂τBw+T̂τε.
Then(2.21)β̂SGLSE=(X̃τ̂′X̃τ̂)-1X̃τ̂′ỹτ̂,ŵSGLSE=(Bτ̂'Bτ̂+nξK)-1Bτ̂'(yτ̂-Xτ̂β̂SGLSE),
where X̃τ̂=(I-S)Xτ̂ and ỹτ̂=(I-S)yτ̂, with yτ̂=T̂τy and Xτ̂=T̂τX. 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𝚺1=limn→∞n-1X̃′V-1X̃ is nonsingular. Let β0 denote the true value of β, then
(2.22)n(β̂-β0)=n(β̂SGLSE-β0)+Op((τn)1/2),(2.23)n(β̂SGLSE-β0)→DN(0,𝚺1-1),
where →D denotes convergence in distribution and β̂=(Xτ'Xτ)-1Xτ'yτ. Assume that 𝚺2=limn→∞n-1B′V-1B is nonsingular and let w0 denote the true value of w, then one has
(2.24)n(ŵ-w0)=n(ŵSGLSE-w0)+Op((τn)1/2),(2.25)n(ŵSGLSE-w0)→DN(0,𝚺2-1),
where ŵ=(Bτ'Bτ+nξK)-1Bτ'(yτ-Xτβ̂).


### 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 θ=(β′,w′)′ in Section 2.3. In Section 2.4 we present the information criteria for the model selection.From now on, we assume that the matricesXτ and Bτ 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 w, which we will denote by l(β,w). Then expression (2.7) can be written as(2.26)ℒ(β,w)=l(β,w)+nξw′Kw.
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(2.27)𝒮(β,w)=l(β,w)+n(∑i=1dpλ1(|βi|)+∑j=1mpλ2(|wj|))+nξw′Kw.
where  pλi(·) are penalty functions and  λi 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θ=(β′,w′)′ and z=(z1,…,zd+m)′ the true parameters and the estimates, respectively. By taking the hard thresholding penalty function(2.28)pλ(|θ|)=λ2+(|θ|-λ)2I(|θ|<λ),
we obtain the hard thresholding rule(2.29)θ̂=zI(|z|>λ).
The L2 penalty Pλ(|θ|)=λ|θ|2 results in a ridge regression and the L1 penalty Pλ(|θ|)=λ|θ| yields a soft thresholding rule(2.30)θ̂=sgn(z)I(|z|>λ)+.
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 L1 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θ=(β′,w′)′. 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 Tτ by using ε obtained in Step 1.  Then θ̂SGLSEHT is obtained using T̂τ(Step 3). Here the penalty parameters λ, and ξ, and the number of basis functions m are chosen using information criteria that will be discussed in Section 2.4.Step 1.
First we obtainβ̂SOLSE and ŵSOLSE by (2.10) and (2.11), respectively. Then we have the model
(2.31)ŷ=BŵSOLSE+Xβ̂SOLSE+ε.Step 2.
An estimator forTτ is constructed followings the work of Aneiros-Pérez and Vilar-Fernández [4]. We use the residuals ε̂=y-BŵSOLSE-Xβ̂SOLSE to construct an estimate of Tτ using the ordinary least square method applied to the model
(2.32)ε̂i=a1ε̂i-1+⋯+aτε̂i-τ+residuali.
The estimator T̂τ is obtained from Tτ by replacing parameters with their estimates.Step 3.
Our SGLSE ofθ is obtained by using the model
(2.33)yτ̂=Bτ̂w+Xτ̂β+ετ̂,
where yτ̂=T̂τy, Bτ̂=T̂τB,  Xτ̂=T̂τX, and ετ=T̂τε. Finding the solution of the penalized least-square of (2.27) needs the local quadratic approximation, because the L1 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 θ(0) that is close to the minimizer of (2.27). If θj(0) is very close to 0, then set θ̂j(0)=0. Otherwise they can be locally approximated by a quadratic function as
(2.34)[pλj(θj)]′=pλj′(|θj|)sgn(θj)≈{pλj′(|θj(0)|)|θj(0)|}θj,
when θ̂j(0)≠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
(2.35)l(β0,w0)+∇lβ(β,w)′(β-β0)+∇lw(β0,w0)′(w-w0)+12(β-β0)′∇ββ2(β0,w0)(β-β0)+12(w-w0)′∇ww2(β0,w0)(w-w0)12(β0-β)′∇βw2(w-w0)+nβ′𝚺λ1(β0)β+nw′𝚺λ2(w0)w,,
where
(2.36)∇lβ(β0,w0)=∂l(β0,w0)∂β,∇lw(β0,w0)=∂l(β0,w0)∂w,∇2lββ(β0,w0)=∂2l(β0,w0)∂β∂β′,∇2lww(β0,w0)=∂2l(β0,w0)∂w∂w′,∇2lβ,w(β0,w0)=∂2l(β0,w0)∂β∂w,𝚺λ1(β0)=diag{pλ1′(|β1(0)|)|β1(0)|,…,pλ1′(|βd(0)|)|βd(0)|},𝚺λ2(w0)=diag{pλ2′(|w1(0)|)|w1(0)|,…,pλ2′(|wm(0)|)|wm(0)|}.
The solution can be found by iteratively computing the block matrix equation:
(2.37)(Xτ̂'Xτ̂+n𝚺λ1(β(0))Xτ̂'Bτ̂Bτ̂'Xτ̂Bτ̂'Bτ̂+αK+n𝚺λ2(w(0)))(βw)=(Xτ̂'Bτ̂')y.
This gives the estimators
(2.38)β̂SGLSEHT=(X̃τ̂′X̃τ̂+n𝚺λ1(β(0)))-1X̃τ̂′ỹτ̂,ŵSGLSEHT=(Bτ̂'Bτ̂+nξK+n𝚺λ2(w(0)))(-1)Bτ̂'(yτ̂-Xτ̂β̂SGLSEHT),
where ỹτ̂=(I-Sτ̂)yτ̂, X̃τ̂=(I-Sτ̂)Xτ̂, and Sτ̂=Bτ̂(Bτ̂'Bτ̂+nξK+n𝚺λ2(w(0)))-1Bτ̂'.


### 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)yτ=Aτθ+e,
where Aτ=(Xτ,Bτ) and θ=(β′,w′)′. In many applications, the number of basis functions m needs to be large to adequately capture the trend. To determine the number of basis functions, all models with m≤mmax are fitted and the preferred model minimizes some model selection criteria.The Schwarz BIC is given by(2.40)BIC=nlog(2πσ̂e2)+logn(thenumberofparameters),
where σ̂e2 is the least-square estimate of σe2 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 trSβ, we formally obtain information criteria for the basis function Gaussian regression model in the form(2.41)BICm=nlog(2πσ̂e2)+(trSθ)logn,
where σ̂e2=n-1∥y-Sθy∥2 and(2.42)trSθ=Aτ(Aτ'Aτ+nξK̃+n𝚺λ(θ))-1Aτ'.
Here 𝚺λ(θ) is defined by (2.44) in what follows.We also consider the use of the BICp criterion to choose appropriate values for these unknown parameters. Denote
(2.43)𝚺λ1(β)=diag{pλ1′′(|β10|),…,pλ1′′(|βd0|)},𝚺λ2(w)=diag{pλ2′′(|w10|),…,pλ2′′(|wm0|)},(2.44)𝚺λ(θ)=(𝚺λ1(β),𝚺λ2(w)).
Let N1 and N2 be the number of zero components in β0 and w0, respectively. Then the BICp criterion is(2.45)BICp=nlog(2πσ̂e2)+nθ̂′𝚺λ(θ)θ̂+nξθ̂K̃θ̂+log|JG(θ̂)|-log|K̃|+-log|Σλ(θ)|+-(m-N2)logξ+Const,
where JG(θ̂) is the (d+m+1)×(d+m+1) matrix of second derivatives of the penalized likelihood defined by(2.46)JG=1nσ̂e2(Aτ'Aτ+nΣλ(θ)+nξK̃Aτ'Λ1nσe21n'ΛAτσ̂e2n2σ̂e2).
Here Λ is a diagonal matrix with ith element Λi=diag[e1,…,en] and 1n=(1,…,1)′. The n-dimensional vector q has ith element (Tijyj-Aτ,ijθj)2/2σ̂e4-1/2σ̂e2 where Tij is the element in the ith row and jth column of Tτ. Also K̃ is the (d+m)×(d+m) matrix defined by(2.47)K̃=(KOd,mOm,dOm,m),
and |K̃|+ and |Σλ(θ)|+ are the product of the (m-N1) and (d+m-N1-N2) nonzero eigenvalues of K̃ 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(2.48)GIC=nlog(2πσ̂e2)+n+2tr{IGJG-1},
where IG is a (m+d+1)×(m+d+1) matrix.  Also IG is basically the product of the empirical influence function and the score function. It is defined by
(2.49)IG=1nσ̂e2(Aτ'Λσe2-𝚺λ(θ)θ̂1n'-ξK̃θ̂1n'q′)(ΛAτ,σ̂e2q).The number of basis functionsm, penalty parameters ξ,λ1,λ2 are determined by minimizing BICm, BICp 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 functionspλ1j′(·) and pλ2j′(·) are negative and nondecreasing with pλ1j′(0)=pλ2j′(0)=0. Let β0 and w0 denote the true values of β and w, respectively. Also let(2.50)a1n=maxj{|pλ1j′(|βj0|)|:βj0≠0},a2n=maxj{|pλ2j′(|wj0|)|:wj0≠0},b1n=maxj{|pλ1j′′(|βj0|)|:βj0≠0},b2n=maxj{|pλ2j′′(|wj0|)|:wj0≠0}.Theorem 2.2.
Under the conditions of Theorem2.1, if a1n,b1n,a2n, and b2n tend to 0 as n→∞, then with probability tending to 1, there exist local minimizers β̂ and ŵ of ℒ(β,w) such that ∥β̂SGLSEHT-β0∥=Op(n-1/2+a1n) and ∥ŵSLOSEHT-w0∥=Op(n-1/2+a2n).Theorem2.2 demonstrates how the rate of convergence of the penalized least-square estimator θ̂SGLSEHT=(β̂SGLSE'HT,ŵSGLSE'HT)′ of ℒ(θ) depends on λij for i=1,2. To achieve the root n convergence rate, we have to take λij small enough so that an=Op(n-1/2).Next we establish the oracle property for the penalized least-square estimator. LetβS1 consist of all nonzero components of β0 and let βN1 consist of all zero components. Let wS2 consist of all nonzero components of w0 and let wN2 consist of all zero components. Let(2.51)x̂(t)′β0=x̂S1′βS1+x̂N1′(t)βN1=x̂S1(t)′βS1,ϕ(t)′w0=ϕS2′wS2+ϕN2′(t)wN2=ϕS2(t)′wS2.
Write(2.52)bβ=(pλ1n′(|β10|)sgn(β10),…,pλS1n′(|βS10|)sgn(βS10))′,bw=(pλ2n′(|w10|)sgn(w10),…,pλS2n′(|wS20|)sgn(wS20))′.
Further, let β̂1 consist of the first S1 components of β̂ and let β̂2 consist of the last d-S1 components of β̂SGLSEHT. Let ŵ1 consist of the first S2 components of ŵ and let ŵ2 consist of the last m-S2 components of ŵSGLSEHT.Theorem 2.3.
Assume that forj=1,…,d and k=1,…,m, one has λ1→0, nλ1→∞, λ2→0 and nλ2→∞. Assume that the penalty functions pλ1′(|βj|) and pλ2′(|wk|) satisfy
(2.53)lim infn→∞lim infβj→0+pλ1′(βj)λ1>0,lim infn→∞lim infωk→0+pλ2′(ωk)λ2>0.
If a1n=a2n=Op(n-1/2) then, under the conditions of Theorem 2.1, with probability tending to 1, the root n consistent local minimizers β̂SGLSEHT=(β̂1′,β̂2′)′ and ŵSGLSEHT=(ŵ1′,ŵ2′)′ in Theorem 2.2 must satisfy the following: (1)
(sparsity)β̂2=ŵ2=0;(2)
(asymptotic normality)(2.54)n(IS1+𝚺λ1(β))(β̂1-β10+(IS1+𝚺λ1(β))-1bβ)→NS1(0,𝚺1(1)-1),n(IS2+𝚺λ1(w))(ŵ1-w10+(IS2+𝚺λ2(w)+ξK)-1bw)→NS2(0,𝚺2(1)-1)Here 𝚺1(1)-1 and 𝚺2(1)-1 consist of the first S1 and S2 rows and columns of 𝚺1 and 𝚺2 defined in Theorem 2.1, respectively.


## 2.1. The Model and Penalized Estimation

We consider the semiparametric regression model:(2.1)yi=α(ti)+β′xi+εi,i=1,…,n,
where yi is the response variable and xi is the d×1 covariate vector at time i, α(ti) is an unspecified baseline function of ti with ti=i/n, β is a vector of unknown regression coefficients, and εi is a Gaussian and zero mean covariance stationary process.We assume the following properties for the error termsεi and vectors of explanatory variables xi.(A.1)
It holds that{εi} is a linear process given by
(2.2)εi=∑j=0∞bjei-j,
where b0=1 and {ei} is an i.i.d. Gaussian random variable with E{ei}=0 and E{ei2}=σe2.(A.2)
The coefficientsbj satisfy the conditions that for all |z|<1, ∑j=0∞bjzj≠0 and ∑j=0∞j2|bj|<∞.We defineγ(k)=cov(εt,εt+k)=E{εtεt+k}.The assumptions on covariate variables are as follows.(B.1)
Alsoxi=(xi1,…,xid)′∈ℝd and {xij},j=1,…d, have mean zero and variance 1.The trend functionα(ti) is expressed as a linear combination of a set of m underlying basis functions :(2.3)α(ti)=∑k=1mwkϕk(ti)=w′ϕ(t),
where {ϕ(ti)=(ϕ1(ti),…,ϕm(ti))′} is an m-dimensional vector constructed from basis functions {ϕk(ti);k=1,…,m}, and w=(w1,…,wm)′ 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)ϕ(ti)=(ti,…,tip,(ti-κ1)+p,…,(ti-κk)+p)′,
where {κk}k=1,…,K are spline knots. This specification uses the so-called truncated power function basis. The choice of the number of knots K 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(2.5)ϕk(ti)=exp(-∥ti-μk∥22sk2),
where μk determines the location and sk2 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)y=Xβ+Bw+ε,
where X=(x1,…,xn)′, y=(y1,…,yn)′, B=(ϕ1,…,ϕn)′ with ϕi=(ϕ1(i/n),…,ϕm(i/n))′. The penalized least-square estimator is then a minimizer of the function(2.7)12(y-Xβ-Bw)′(y-Xβ-Bw)+nξw′Kw,
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  K is an appropriate positive semidefinite symmetric matrix. For example, if K satisfies w′Kw=∫01[α′′(u)]2du, we have the usual quadratic integral penalty (see, e.g., Green and Silverman [19]). By simple calculus, (2.7) is minimized when β and w satisfy the block matrix equation(2.8)(X′XX′BB′XB′B+nξK)(βw)=(X′B′)y.
This equation can be solved without any iteration (see, e.g., Green [20]). First, we find Bw̃=S(y-Xβ), where S=B(B′B+αK)-1B′ is usually called the smoothing matrix. Substituting Bw̃ into (2.6), we obtain(2.9)ỹ=X̃β+ε,
where ỹ=(I-S)y, X̃=(I-S)X ,and I is the identity matrix of order n. Applying least-square to the linear model (2.9), we obtain the semiparametric ordinary least-square estimator (SOLSE) result:(2.10)β̂SOLSE=(X̃′X̃)-1X̃′ỹ,(2.11)ŵSOLSE=(B′B+nξK)-1B′(y-Xβ̂SOLSE).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 {εi} in (2.6) are invertible. Let b(L)=∑j=1∞bjei-j, where L is the lag operator and a(L)=b(L)-1=a0-∑j=1∞ajLj with a0=1. Applying a(L) to the model (2.6) and rewriting the corresponding equation, we obtain the new model:(2.12)y̲=X̲β+B̲w+e,
where y̲=(y̲1,…,y̲n)′,X̲=(x̲1,…,x̲n)′,B̲=(ϕ1̲,…,ϕn̲)′ and e=(e1,…,en)′. Here
(2.13)y̲i=yi-∑j=1∞ajyi-j,ϕ̲i=ϕi-∑j=1∞ajϕi-j,x̲i=xi-∑j=1∞ajxi-j.
The regression errors in (2.12) are i.i.d. Because, in practice, the response variable yi is unknown, we use a reasonable approximation by y̲i based on the work by Xiao et al. [22] and Aneiros-Pérez and Vilar-Fernández [11].Under the usual regularity conditions the coefficientsaj decrease geometrically so, letting τ=τ(n) denote a truncation parameter, we may consider the truncated autoregression on εi:(2.14)ei=εi-∑j=1∞ajεi-j,
where ei are i.i.d. random variables with E(ei)=0. We make the following assumption about the truncation parameter.(C.1)
The truncation parameterτ satisfies τ(n)=clogn for some c>0.The expansion rate of the truncation parameter given in (C.1) is also for convenience. LetTτ be the n×n transformation matrix such that eτ=Tτε. Then the model (2.12) can be expressed as(2.15)Tτy=TτXβ+TτBw+Tτε,
where(2.16)Tτ=(δ110⋯0δ21-δ220⋯0⋮δτ1⋯-δττ-aτ⋯-a110-aτ⋯-a11⋮0⋯0-aτ⋯-a11)
with δ11=σe/γ(0),δ22=σe/(1-ρ2(1))γ(0),δ21=ρ(1)(σe/(1-ρ2(1))γ(0)),…. Here ρ(h)=γ(h)/ρ(0) denotes the lag h autocorrelation function of {εi}.Now our estimation problem for the semiparametric time series regression model can be expressed as the minimization of the function(2.17)ℒ(β,w)=12(y-Xβ-Bw)′V-1(y-Xβ-Bw)+αw′Kw,
where V-1=σe-2Tτ'Tτ and σe2=n-1∥Tτε∥2. Based on the work by Aneiros-Pérez and Vilar-Fernández [11], an estimator for Tτ is constructed as follows. We use the residuals ε̂=y-Xβ̂SOLSE-BŵSOLSE to construct an estimate of Tτ using the ordinary least square method applied to the model(2.18)ε̂i=a1ε̂i-1+⋯+aτε̂i-τ+residuali.
Define the estimate âτ=(â1,â2,…,âτ)′ of aτ=(a1,a2,…,aτ)′, where
(2.19)âτ=(Êτ′Êτ)-1Êτ′ε̂,
where ε̂=(ε̂τ+1,…,ε̂n) and Êτ is the (n-τ)×τ matrix of regressors with the typical element ε̂i-j. Then T̂τ is obtained from Tτ by replacing aj with âj, σe2 with σ̂e2, and so forth. Applying least-square to the linear model, we obtain
(2.20)T̂τy=T̂τXβ+T̂τBw+T̂τε.
Then(2.21)β̂SGLSE=(X̃τ̂′X̃τ̂)-1X̃τ̂′ỹτ̂,ŵSGLSE=(Bτ̂'Bτ̂+nξK)-1Bτ̂'(yτ̂-Xτ̂β̂SGLSE),
where X̃τ̂=(I-S)Xτ̂ and ỹτ̂=(I-S)yτ̂, with yτ̂=T̂τy and Xτ̂=T̂τX. 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𝚺1=limn→∞n-1X̃′V-1X̃ is nonsingular. Let β0 denote the true value of β, then
(2.22)n(β̂-β0)=n(β̂SGLSE-β0)+Op((τn)1/2),(2.23)n(β̂SGLSE-β0)→DN(0,𝚺1-1),
where →D denotes convergence in distribution and β̂=(Xτ'Xτ)-1Xτ'yτ. Assume that 𝚺2=limn→∞n-1B′V-1B is nonsingular and let w0 denote the true value of w, then one has
(2.24)n(ŵ-w0)=n(ŵSGLSE-w0)+Op((τn)1/2),(2.25)n(ŵSGLSE-w0)→DN(0,𝚺2-1),
where ŵ=(Bτ'Bτ+nξK)-1Bτ'(yτ-Xτβ̂).


## 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 θ=(β′,w′)′ in Section 2.3. In Section 2.4 we present the information criteria for the model selection.From now on, we assume that the matricesXτ and Bτ 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 w, which we will denote by l(β,w). Then expression (2.7) can be written as(2.26)ℒ(β,w)=l(β,w)+nξw′Kw.
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(2.27)𝒮(β,w)=l(β,w)+n(∑i=1dpλ1(|βi|)+∑j=1mpλ2(|wj|))+nξw′Kw.
where  pλi(·) are penalty functions and  λi 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θ=(β′,w′)′ and z=(z1,…,zd+m)′ the true parameters and the estimates, respectively. By taking the hard thresholding penalty function(2.28)pλ(|θ|)=λ2+(|θ|-λ)2I(|θ|<λ),
we obtain the hard thresholding rule(2.29)θ̂=zI(|z|>λ).
The L2 penalty Pλ(|θ|)=λ|θ|2 results in a ridge regression and the L1 penalty Pλ(|θ|)=λ|θ| yields a soft thresholding rule(2.30)θ̂=sgn(z)I(|z|>λ)+.
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 L1 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θ=(β′,w′)′. 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 Tτ by using ε obtained in Step 1.  Then θ̂SGLSEHT is obtained using T̂τ(Step 3). Here the penalty parameters λ, and ξ, and the number of basis functions m are chosen using information criteria that will be discussed in Section 2.4.Step 1.
First we obtainβ̂SOLSE and ŵSOLSE by (2.10) and (2.11), respectively. Then we have the model
(2.31)ŷ=BŵSOLSE+Xβ̂SOLSE+ε.Step 2.
An estimator forTτ is constructed followings the work of Aneiros-Pérez and Vilar-Fernández [4]. We use the residuals ε̂=y-BŵSOLSE-Xβ̂SOLSE to construct an estimate of Tτ using the ordinary least square method applied to the model
(2.32)ε̂i=a1ε̂i-1+⋯+aτε̂i-τ+residuali.
The estimator T̂τ is obtained from Tτ by replacing parameters with their estimates.Step 3.
Our SGLSE ofθ is obtained by using the model
(2.33)yτ̂=Bτ̂w+Xτ̂β+ετ̂,
where yτ̂=T̂τy, Bτ̂=T̂τB,  Xτ̂=T̂τX, and ετ=T̂τε. Finding the solution of the penalized least-square of (2.27) needs the local quadratic approximation, because the L1 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 θ(0) that is close to the minimizer of (2.27). If θj(0) is very close to 0, then set θ̂j(0)=0. Otherwise they can be locally approximated by a quadratic function as
(2.34)[pλj(θj)]′=pλj′(|θj|)sgn(θj)≈{pλj′(|θj(0)|)|θj(0)|}θj,
when θ̂j(0)≠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
(2.35)l(β0,w0)+∇lβ(β,w)′(β-β0)+∇lw(β0,w0)′(w-w0)+12(β-β0)′∇ββ2(β0,w0)(β-β0)+12(w-w0)′∇ww2(β0,w0)(w-w0)12(β0-β)′∇βw2(w-w0)+nβ′𝚺λ1(β0)β+nw′𝚺λ2(w0)w,,
where
(2.36)∇lβ(β0,w0)=∂l(β0,w0)∂β,∇lw(β0,w0)=∂l(β0,w0)∂w,∇2lββ(β0,w0)=∂2l(β0,w0)∂β∂β′,∇2lww(β0,w0)=∂2l(β0,w0)∂w∂w′,∇2lβ,w(β0,w0)=∂2l(β0,w0)∂β∂w,𝚺λ1(β0)=diag{pλ1′(|β1(0)|)|β1(0)|,…,pλ1′(|βd(0)|)|βd(0)|},𝚺λ2(w0)=diag{pλ2′(|w1(0)|)|w1(0)|,…,pλ2′(|wm(0)|)|wm(0)|}.
The solution can be found by iteratively computing the block matrix equation:
(2.37)(Xτ̂'Xτ̂+n𝚺λ1(β(0))Xτ̂'Bτ̂Bτ̂'Xτ̂Bτ̂'Bτ̂+αK+n𝚺λ2(w(0)))(βw)=(Xτ̂'Bτ̂')y.
This gives the estimators
(2.38)β̂SGLSEHT=(X̃τ̂′X̃τ̂+n𝚺λ1(β(0)))-1X̃τ̂′ỹτ̂,ŵSGLSEHT=(Bτ̂'Bτ̂+nξK+n𝚺λ2(w(0)))(-1)Bτ̂'(yτ̂-Xτ̂β̂SGLSEHT),
where ỹτ̂=(I-Sτ̂)yτ̂, X̃τ̂=(I-Sτ̂)Xτ̂, and Sτ̂=Bτ̂(Bτ̂'Bτ̂+nξK+n𝚺λ2(w(0)))-1Bτ̂'.


## 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)yτ=Aτθ+e,
where Aτ=(Xτ,Bτ) and θ=(β′,w′)′. In many applications, the number of basis functions m needs to be large to adequately capture the trend. To determine the number of basis functions, all models with m≤mmax are fitted and the preferred model minimizes some model selection criteria.The Schwarz BIC is given by(2.40)BIC=nlog(2πσ̂e2)+logn(thenumberofparameters),
where σ̂e2 is the least-square estimate of σe2 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 trSβ, we formally obtain information criteria for the basis function Gaussian regression model in the form(2.41)BICm=nlog(2πσ̂e2)+(trSθ)logn,
where σ̂e2=n-1∥y-Sθy∥2 and(2.42)trSθ=Aτ(Aτ'Aτ+nξK̃+n𝚺λ(θ))-1Aτ'.
Here 𝚺λ(θ) is defined by (2.44) in what follows.We also consider the use of the BICp criterion to choose appropriate values for these unknown parameters. Denote
(2.43)𝚺λ1(β)=diag{pλ1′′(|β10|),…,pλ1′′(|βd0|)},𝚺λ2(w)=diag{pλ2′′(|w10|),…,pλ2′′(|wm0|)},(2.44)𝚺λ(θ)=(𝚺λ1(β),𝚺λ2(w)).
Let N1 and N2 be the number of zero components in β0 and w0, respectively. Then the BICp criterion is(2.45)BICp=nlog(2πσ̂e2)+nθ̂′𝚺λ(θ)θ̂+nξθ̂K̃θ̂+log|JG(θ̂)|-log|K̃|+-log|Σλ(θ)|+-(m-N2)logξ+Const,
where JG(θ̂) is the (d+m+1)×(d+m+1) matrix of second derivatives of the penalized likelihood defined by(2.46)JG=1nσ̂e2(Aτ'Aτ+nΣλ(θ)+nξK̃Aτ'Λ1nσe21n'ΛAτσ̂e2n2σ̂e2).
Here Λ is a diagonal matrix with ith element Λi=diag[e1,…,en] and 1n=(1,…,1)′. The n-dimensional vector q has ith element (Tijyj-Aτ,ijθj)2/2σ̂e4-1/2σ̂e2 where Tij is the element in the ith row and jth column of Tτ. Also K̃ is the (d+m)×(d+m) matrix defined by(2.47)K̃=(KOd,mOm,dOm,m),
and |K̃|+ and |Σλ(θ)|+ are the product of the (m-N1) and (d+m-N1-N2) nonzero eigenvalues of K̃ 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(2.48)GIC=nlog(2πσ̂e2)+n+2tr{IGJG-1},
where IG is a (m+d+1)×(m+d+1) matrix.  Also IG is basically the product of the empirical influence function and the score function. It is defined by
(2.49)IG=1nσ̂e2(Aτ'Λσe2-𝚺λ(θ)θ̂1n'-ξK̃θ̂1n'q′)(ΛAτ,σ̂e2q).The number of basis functionsm, penalty parameters ξ,λ1,λ2 are determined by minimizing BICm, BICp 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 functionspλ1j′(·) and pλ2j′(·) are negative and nondecreasing with pλ1j′(0)=pλ2j′(0)=0. Let β0 and w0 denote the true values of β and w, respectively. Also let(2.50)a1n=maxj{|pλ1j′(|βj0|)|:βj0≠0},a2n=maxj{|pλ2j′(|wj0|)|:wj0≠0},b1n=maxj{|pλ1j′′(|βj0|)|:βj0≠0},b2n=maxj{|pλ2j′′(|wj0|)|:wj0≠0}.Theorem 2.2.
Under the conditions of Theorem2.1, if a1n,b1n,a2n, and b2n tend to 0 as n→∞, then with probability tending to 1, there exist local minimizers β̂ and ŵ of ℒ(β,w) such that ∥β̂SGLSEHT-β0∥=Op(n-1/2+a1n) and ∥ŵSLOSEHT-w0∥=Op(n-1/2+a2n).Theorem2.2 demonstrates how the rate of convergence of the penalized least-square estimator θ̂SGLSEHT=(β̂SGLSE'HT,ŵSGLSE'HT)′ of ℒ(θ) depends on λij for i=1,2. To achieve the root n convergence rate, we have to take λij small enough so that an=Op(n-1/2).Next we establish the oracle property for the penalized least-square estimator. LetβS1 consist of all nonzero components of β0 and let βN1 consist of all zero components. Let wS2 consist of all nonzero components of w0 and let wN2 consist of all zero components. Let(2.51)x̂(t)′β0=x̂S1′βS1+x̂N1′(t)βN1=x̂S1(t)′βS1,ϕ(t)′w0=ϕS2′wS2+ϕN2′(t)wN2=ϕS2(t)′wS2.
Write(2.52)bβ=(pλ1n′(|β10|)sgn(β10),…,pλS1n′(|βS10|)sgn(βS10))′,bw=(pλ2n′(|w10|)sgn(w10),…,pλS2n′(|wS20|)sgn(wS20))′.
Further, let β̂1 consist of the first S1 components of β̂ and let β̂2 consist of the last d-S1 components of β̂SGLSEHT. Let ŵ1 consist of the first S2 components of ŵ and let ŵ2 consist of the last m-S2 components of ŵSGLSEHT.Theorem 2.3.
Assume that forj=1,…,d and k=1,…,m, one has λ1→0, nλ1→∞, λ2→0 and nλ2→∞. Assume that the penalty functions pλ1′(|βj|) and pλ2′(|wk|) satisfy
(2.53)lim infn→∞lim infβj→0+pλ1′(βj)λ1>0,lim infn→∞lim infωk→0+pλ2′(ωk)λ2>0.
If a1n=a2n=Op(n-1/2) then, under the conditions of Theorem 2.1, with probability tending to 1, the root n consistent local minimizers β̂SGLSEHT=(β̂1′,β̂2′)′ and ŵSGLSEHT=(ŵ1′,ŵ2′)′ in Theorem 2.2 must satisfy the following: (1)
(sparsity)β̂2=ŵ2=0;(2)
(asymptotic normality)(2.54)n(IS1+𝚺λ1(β))(β̂1-β10+(IS1+𝚺λ1(β))-1bβ)→NS1(0,𝚺1(1)-1),n(IS2+𝚺λ1(w))(ŵ1-w10+(IS2+𝚺λ2(w)+ξK)-1bw)→NS2(0,𝚺2(1)-1)Here 𝚺1(1)-1 and 𝚺2(1)-1 consist of the first S1 and S2 rows and columns of 𝚺1 and 𝚺2 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)yi=α(ti)+β′xi+εi,
where α(ti)=exp(-3(i/n))sin(3πi/n), β=(3,1.5,0,0.2,0,0,0)′  and ε(t) 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 x from a normal distribution with mean 0 and cov(xi,xj)=0.5|i-j|. In each case, the autoregressive coefficient is set to 0, 0.25, 0.5 or 0.75 and the sample size n is set to 50, 100 or 200. Figure 1 depicts some examples of simulation data.Simulation data with (a)n=50 and ρ=0.5, (b) n=100 and ρ=0.5, (c) n=200 and ρ=0.5. The dotted lines represent α(t); the solid lines α(t)+ε(t).
(a)(b)(c)We compare the effectiveness of our proposed procedure (PLS + HT) with an existing procedure (PLS). We also compare the performance of the information criteria BICm, GIC and BICp 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α̂(t) by the square root of average squared errors (RASEα):(3.2)RASEα=ngrid-1∑k=1ngrid{α̂(tk)-α(tk)}2,
where {tk,k=1,…,ngrid} are the grid points at which the baseline function α(·) is estimated. In our simulation, we use ngrid=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 BICp work better than those based on BICm or GIC.Table 1
Means (standard deviations) ofRASEα.
ρ=0.0ρ=0.25ρ=0.50ρ=0.75n=50PLSBICm0.069 (0.013)0.081(0.017)0.114 (0.037)0.232 (0.115)GIC0.047 (0.013)0.062 (0.015)0.106 (0.037)0.229 (0.120)BICp0.042 (0.010)0.060 (0.019)0.103 (0.039)0.226 (0.124)PLS + HTBICm0.061 (0.040)0.070 (0.021)0.101 (0.038)0.226 (0.103)GIC0.053 (0.017)0.068 (0.020)0.101 (0.034)0.218 (0.097)BICp0.046 (0.015)0.060 (0.019)0.093 (0.034)0.214 (0.101)n=100PLSBICm0.041(0.008)0.052 (0.012)0.080 (0.025)0.172 (0.080)GIC0.034 (0.008)0.044 (0.011)0.074 (0.026)0.170 (0.080)BICp0.036 (0.010)0.044 (0.010)0.070 (0.024)0.163 (0.079)PLS + HTBICm0.042 (0.008)0.051 (0.016)0.080 (0.024)0.172 (0.079)GIC0.040 (0.015)0.048 (0.016)0.073 (0.024)0.168 (0.078)BICp0.037 (0.011)0.041 (0.011)0.068 (0.023)0.158 (0.075)n=200PLSBICm0.029 (0.005)0.040 (0.016)0.058 (0.018)0.129 (0.056)GIC0.025 (0.008)0.033 (0.010)0.056 (0.018)0.125 (0.057)BICp0.029 (0.006)0.031 (0.007)0.050 (0.015)0.114 (0.052)PLS + HTBICm0.030 (0.005)0.040 (0.016)0.058 (0.019)0.127 (0.053)GIC0.027 (0.009)0.033 (0.011)0.054 (0.015)0.123 (0.054)BICp0.019 (0.008)0.028 (0.009)0.047 (0.018)0.109 (0.048)Then the performance ofβ̂ is assessed by the square root of average squared errors (RASEβ):(3.3)RASEβ=1d∑i=1d(β̂i-βi)2.
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, BICp 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) ofRASEβ.
ρ=0.0ρ=0.25ρ=0.50ρ=0.75n=50PLSBICm0.022 (0.007)0.023 (0.007)0.022 (0.007)0.020 (0.007)GIC0.021 (0.006)0.023 (0.007)0.023 (0.010)0.021 (0.007)BICp0.021(0.006)0.022 (0.007)0.022 (0.009)0.020 (0.007)PLS + HTBICm0.011(0.005)0.013 (0.007)0.012 (0.007)0.010 (0.005)GIC0.010 (0.004)0.013 (0.007)0.013 (0.009)0.011(0.006)BICp0.010 (0.004)0.011 (0.005)0.011 (0.006)0.010 (0.005)n=100PLSBICm0.014 (0.004)0.014 (0.004)0.014 (0.005)0.012 (0.004)GIC0.013 (0.004)0.014 (0.004)0.013 (0.004)0.012 (0.004)BICp0.014 (0.004)0.014 (0.004)0.013 (0.004)0.011 (0.004)PLS + HTBICm0.007 (0.003)0.008 (0.004)0.007 (0.004)0.006 (0.003)GIC0.007 (0.003)0.008 (0.004)0.007 (0.003)0.006 (0.003)BICp0.007 (0.003)0.007 (0.003)0.006 (0.003)0.006 (0.003)n=200PLSBICm0.009 (0.003)0.009 (0.003)0.009 (0.003)0.007 (0.002)GIC0.009 (0.003)0.009 (0.003)0.008 (0.003)0.007 (0.002)BICp0.009 (0.003)0.009 (0.003)0.008 (0.002)0.007 (0.002)PLS + HTBICm0.004 (0.002)0.005 (0.002)0.005 (0.002)0.005 (0.002)GIC0.005 (0.002)0.005 (0.002)0.005 (0.002)0.004 (0.002)BICp0.005 (0.002)0.005 (0.002)0.004 (0.002)0.005 (0.002)The first step ahead prediction error (PE), which is defined as(3.4)PE=(ŷn+1-yn+1∣n)2,
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 BICm perform well for large sample sizes.Table 3
Means (standard deviations) of first step ahead prediction errors.
ρ=0.0ρ=0.25ρ=0.50ρ=0.75n=50PLSBICm0.136 (0.115)0.150 (0.116)0.140 (0.120)0.158 (0.117)GIC0.111 (0.088)0.127 (0.097)0.134 (0.098)0.149 (0.122)BICp0.111 (0.088)0.127 (0.097)0.131 (0.095)0.149 (0.122)PLS + HTBICm0.121 (0.096)0.106 (0.086)0.119 (0.092)0.139 (0.112)GIC0.094 (0.071)0.118 (0.093)0.126 (0.094)0.139 (0.112)BICp0.095 (0.071)0.116 (0.092)0.124 (0.093)0.139 (0.112)n=100PLSBICm0.101 (0.086)0.105 (0.082)0.130 (0.112)0.145 (0.124)GIC0.090 (0.070)0.101 (0.078)0.105 (0.082)0.137 (0.109)BICp0.091 (0.070)0.096 (0.072)0.105 (0.092)0.137 (0.109)PLS + HTBICm0.097 (0.082)0.096 (0.078)0.098 (0.088)0.140 (0.162)GIC0.084 (0.063)0.091 (0.071)0.103 (0.081)0.130 (0.111)BICp0.084 (0.063)0.091 (0.071)0.103 (0.081)0.130 (0.111)n=200PLSBICm0.091 (0.070)0.105 (0.081)0.114 (0.087)0.174 (0.129)GIC0.087 (0.068)0.095 (0.072)0.102 (0.077)0.139 (0.114)BICp0.086 (0.068)0.095 (0.072)0.102 (0.077)0.139 (0.114)PLS + HTBICm0.084 (0.066)0.090 (0.069)0.091 (0.068)0.123 (0.096)GIC0.083 (0.063)0.090 (0.069)0.098 (0.076)0.126 (0.100)BICp0.082 (0.063)0.092 (0.070)0.098 (0.076)0.126 (0.100)The means and standard deviations of the number and deviation of basis functions are shown in Tables4 and 5. The BICp gives a smaller number of basis functions than the other information criteria. The models evaluated by BICp also give smaller standard deviations of the number of basis functions. The models determined by BICp tend to choose larger deviations of basis functions than those based on BICm 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 n. From Table 5, it also appears that the deviations of basis functions do not depend on the sample size n and ρ.Table 4
Means (standard deviations) of the number of basis functions.
ρ=0.0ρ=0.25ρ=0.50ρ=0.75n=50BICm7.87 (1.38)8.85 (1.13)8.76 (1.24)8.82 (1.17)GIC8.06 (1.44)8.75 (1.27)8.84 (1.20)8.84 (1.24)BICp6.02 (0.14)6.15 (0.53)6.17 (0.37)6.21 (0.48)n=100BICm7.98 (1.31)8.83 (1.17)8.71 (1.30)8.71 (1.30)GIC8.01 (1.37)8.91 (1.18)8.67 (1.29)8.95 (1.20)BICp6.20 (0.50)6.22 (0.44)6.31 (0.60)6.35 (0.66)n=200BICm7.93 (1.33)8.18 (1.44)8.25 (1.48)8.20 (1.39)GIC8.11 (1.35)8.11 (1.52)8.39 (1.41)8.55 (1.37)BICp6.15 (0.66)6.22 (0.73)6.46 (1.03)6.93 (1.43)Table 5
Means (standard deviations) of the deviations of basis functions.
ρ=0.0ρ=0.25ρ=0.50ρ=0.75n=50BICm0.10 (0.02)0.10 (0.02)0.10 (0.02)0.10 (0.02)GIC0.11 (0.03)0.10 (0.03)0.10 (0.03)0.10 (0.03)BICp0.14 (0.02)0.18 (0.03)0.16 (0.03)0.16 (0.03)n=100BICm0.10 (0.02)0.09 (0.02)0.09 (0.02)0.09 (0.03)GIC0.11 (0.03)0.09 (0.02)0.10 (0.03)0.09 (0.02)BICp0.15 (0.02)0.15 (0.04)0.15 (0.03)0.13 (0.03)n=200BICm0.10 (0.02)0.11 (0.03)0.11 (0.03)0.10 (0.03)GIC0.11 (0.03)0.12 (0.04)0.11 (0.04)0.10 (0.03)BICp0.15 (0.03)0.17 (0.02)0.16 (0.03)0.14 (0.04)We now compare the performance of our procedure with existing procedures in terms of the reduction of model complexity. Table6 shows simulation results of the means and standard deviations of the number of parameters excluded (β=0 or w=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 BICp 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 BICp give a larger number of excluded parameters as the sample size increases. On the other hand, the models evaluated by BICm 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.
ρ=0.0ρ=0.25ρ=0.50ρ=0.75n=50PLS + HTBICm7.715 (0.915)6.910 (1.087)7.300 (1.364)6.888 (1.343)GIC8.345 (1.568)7.404 (1.850)7.620 (1.715)7.337 (1.598)BICp4.950 (0.419)5.020 (0.502)5.070 (0.492)5.092 (0.421)n=100PLS + HTBICm7.506 (0.784)7.334 (1.251)5.698 (0.772)5.460 (0.700)GIC7.916 (1.239)7.718 (1.435)5.906 (0.919)5.740 (0.866)BICp4.990 (0.184)5.076 (0.332)5.092 (0.316)5.086 (0.327)n=200PLS + HTBICm7.062 (0.723)5.594 (0.744)5.544 (0.736)5.460 (0.702)GIC7.450 (1.116)5.764 (0.847)5.656 (0.864)5.586 (0.802)BICp5.008 (0.109)5.152 (0.359)5.162 (0.385)5.086 (0.356)Table7 shows the means and standard deviations of the number of basis functions excluded as w=0 by the proposed procedure. From Table 7 it appears that the models evaluated by BICp 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 BICp 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.
ρ=0.0ρ=0.25ρ=0.50ρ=0.75n=50BICm3.52 (2.29)4.21 (2.23)3.98 (1.60)3.96 (1.49)GIC3.74 (2.15)4.40 (1.90)4.18 (1.51)4.26 (1.46)BICp1.03 (0.22)1.20 (0.60)1.28 (0.54)1.24 (0.49)n=100BICm3.35 (2.19)4.49 (2.04)3.78 (1.58)3.95 (1.60)GIC3.67 (2.15)4.62 (1.84)3.91 (1.53)4.30 (1.60)BICp1.06 (0.31)1.78 (0.96)1.31 (0.60)1.36 (0.66)n=200BICm3.64 (2.13)3.26 (1.71)3.26 (1.71)3.61 (1.60)GIC3.86 (2.02)3.43 (1.81)3.65 (1.69)3.89 (1.76)BICp1.12 (0.34)1.23 (0.75)1.46 (1.03)1.93 (1.44)Table8 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 BICp 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.
ρ=0.0ρ=0.25ρ=0.50ρ=0.75n=50BICm4.14 (1.60)4.15 (1.63)4.69 (0.92)4.79 (0.74)GIC4.28 (1.47)4.35 (1.41)4.70 (0.89)4.72 (0.87)BICp4.97 (0.21)4.95 (0.26)4.97 (0.23)4.99 (0.14)n=100BICm4.15 (1.59)4.17 (1.55)4.72 (0.92)4.77 (0.87)GIC4.22 (1.51)4.29 (1.47)4.77 (0.84)4.65 (1.03)BICp4.98 (0.14)4.95 (0.26)5.00 (0.04)5.00 (0.06)n=200BICm4.14 (1.59)4.78 (0.82)4.78 (0.82)4.72 (0.86)GIC4.16 (1.55)4.68 (1.01)4.75 (0.88)4.66 (1.04)BICp4.99 (0.11)4.99 (0.15)5.00 (0.00)5.00 (0.04)Table9 shows the percentage of times that various βi were estimated as being zero. As for the parameters βj≠0,j=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 BICp is superior to the other criteria.Table 9
Percentage of timesβi is estimated as zero.
β3=0β4=0β6=0β7=0β8=0n=50ρ=0BICm0.840.830.820.830.82GIC0.870.850.870.860.83BICp1.000.990.991.001.00ρ=0.25BICm0.830.830.840.830.83GIC0.860.860.860.890.87BICp0.990.990.990.990.98ρ=0.50BICm0.950.930.940.940.93GIC0.940.940.930.940.95BICp0.990.991.001.000.99ρ=0.75BICm0.960.960.950.950.97GIC0.940.930.950.940.96BICp1.001.001.001.001.00n=100ρ=0BICm0.830.830.840.820.82GIC0.850.840.850.840.84BICp1.000.991.000.991.00ρ=0.25BICm0.830.840.830.820.85GIC0.870.850.880.850.85BICp0.990.990.990.991.00ρ=0.50BICm0.950.930.950.950.94GIC0.960.950.940.960.95BICp1.001.001.001.001.00ρ=0.75BICm0.960.950.950.950.95GIC0.930.940.940.920.92BICp1.001.001.001.001.00n=200ρ=0BICm0.920.930.920.910.94GIC0.940.940.940.950.95BICp1.001.001.001.000.99ρ=0.25BICm0.950.940.940.950.93GIC0.940.940.940.940.93BICp1.001.001.001.001.00ρ=0.50BICm0.970.950.950.960.95GIC0.960.950.950.940.95BICp1.001.001.001.001.00ρ=0.75BICm0.960.940.950.950.93GIC0.930.930.930.940.94BICp1.001.001.001.001.00


## 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 yi is the logarithm of the annual per capita consumption of spirits. The explanatory variables xi1 and xi2 are the logarithms of per capita income and the price of spirits, respectively, and xi3=xi1xi2. Figure 2 shows that there is a change-point at the start of the First World War (1914). Therefore, we prepare a variable z: z=0 from 1870 to 1914 and z=1 from 1915 to 1933. From this we derive another three explanatory variables: xi4=xi1z, xi5=xi2z, and xi6=xi3z. We consider the semiparametric model:Application data-set: (a):yi=log(the annual per capita consumption of spirits); (b): xi1=log(per capita income); (c): xi2=log(price of spirits).
(a)(b)(c)(4.1)yi=α(ti)+β′xi+εi,i=1,…,69.We computed the basis function estimate forα using the existing procedure (PLS) and the proposed procedure (PLS + HT) with BICp. 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 s is estimated as 0.12. Therefore, we obtain the modelTable 10
Estimated Coefficients for Model 4.1.
PLS estimatorsSEPLS + HT estimatorsSEβ1-0.6533.0800β2-1.1215.9620β31.8429.1640β43.5703.7612.3950.421β5-2.5534.4550β6-1.257.763-2.4110.524(4.2)ŷi=α̂(ti)+2.395xi4-2.411xi5,i=1,…,69.
The results indicate that the proposed procedure excludes some variable and reduces model complexity. Table 10 shows that PLS + HT selects only β4 and β6. 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-4.Figure 3
Plots of estimated curves. The solid line representsy; the dotted lines are the estimates of y; 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:(4.3)yi=α(ti)+β1xi1+β2xi2+εi,i=1,…,69.
The semiparametric least-square (SLS) regression gives ŷi=α̂(ti)+0.65xi1-0.95xi2. The residual mean square is 2.2×10-4, 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 ŷi=α̂(ti)-0.71xi2. The result indicates that xi1 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 oflog(TFRi) on log(FLPi). As the law of the Equal Employment Act came into force in 1985, we use the interaction variables “dummy for 1968–1984 ×log(TFR)” (xi2) and for 1985–2007 (xi3). We also use dummy variables for 1990–1999 and 2000–2007 (xi4, xi5) and consider the semiparametric model(4.4)log(TFR)i=α(ti)+β′xi+εi,i=1,…,40.We applied the existing procedure (PLS) and proposed procedure (PLS + HT) withBICp. The resulting estimates and standard errors (SE) of β are given in Table 11.  Therefore, we obtain the modelTable 11
Estimated coefficients for Model 4.4.
log(FLP)PLS estimatorsSEPLS + HT estimatorsSEfor 1968–1984-0.32-1.99-0.36-2.16for 1984–2007-0.28-2.00-0.31-2.18for 1968–19790.020.170for 1980–198901.370for 1990–1999-0.04-0.510for 2000–20070.040.170(4.5)ŷi=α̂(ti)-0.27xi1-0.20xi2,i=1,…,40.
The residual mean square of PLS + HT is 2.24×10-2 and that of PLS is 2.47×10-2. The selected number of basis functions is six with one excluded basis function and the spread parameter s is estimated as 0.30. Table 11 shows that PLS + HT selects only log(FLPi) 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 representsy; the dotted lines are the estimated curves of y; the dashed lines are the estimated curves of α.Figure 8
Plot of residuals.Figure 9
ACF plot of residuals.


## 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 yi is the logarithm of the annual per capita consumption of spirits. The explanatory variables xi1 and xi2 are the logarithms of per capita income and the price of spirits, respectively, and xi3=xi1xi2. Figure 2 shows that there is a change-point at the start of the First World War (1914). Therefore, we prepare a variable z: z=0 from 1870 to 1914 and z=1 from 1915 to 1933. From this we derive another three explanatory variables: xi4=xi1z, xi5=xi2z, and xi6=xi3z. We consider the semiparametric model:Application data-set: (a):yi=log(the annual per capita consumption of spirits); (b): xi1=log(per capita income); (c): xi2=log(price of spirits).
(a)(b)(c)(4.1)yi=α(ti)+β′xi+εi,i=1,…,69.We computed the basis function estimate forα using the existing procedure (PLS) and the proposed procedure (PLS + HT) with BICp. 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 s is estimated as 0.12. Therefore, we obtain the modelTable 10
Estimated Coefficients for Model 4.1.
PLS estimatorsSEPLS + HT estimatorsSEβ1-0.6533.0800β2-1.1215.9620β31.8429.1640β43.5703.7612.3950.421β5-2.5534.4550β6-1.257.763-2.4110.524(4.2)ŷi=α̂(ti)+2.395xi4-2.411xi5,i=1,…,69.
The results indicate that the proposed procedure excludes some variable and reduces model complexity. Table 10 shows that PLS + HT selects only β4 and β6. 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-4.Figure 3
Plots of estimated curves. The solid line representsy; the dotted lines are the estimates of y; 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:(4.3)yi=α(ti)+β1xi1+β2xi2+εi,i=1,…,69.
The semiparametric least-square (SLS) regression gives ŷi=α̂(ti)+0.65xi1-0.95xi2. The residual mean square is 2.2×10-4, 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 ŷi=α̂(ti)-0.71xi2. The result indicates that xi1 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 oflog(TFRi) on log(FLPi). As the law of the Equal Employment Act came into force in 1985, we use the interaction variables “dummy for 1968–1984 ×log(TFR)” (xi2) and for 1985–2007 (xi3). We also use dummy variables for 1990–1999 and 2000–2007 (xi4, xi5) and consider the semiparametric model(4.4)log(TFR)i=α(ti)+β′xi+εi,i=1,…,40.We applied the existing procedure (PLS) and proposed procedure (PLS + HT) withBICp. The resulting estimates and standard errors (SE) of β are given in Table 11.  Therefore, we obtain the modelTable 11
Estimated coefficients for Model 4.4.
log(FLP)PLS estimatorsSEPLS + HT estimatorsSEfor 1968–1984-0.32-1.99-0.36-2.16for 1984–2007-0.28-2.00-0.31-2.18for 1968–19790.020.170for 1980–198901.370for 1990–1999-0.04-0.510for 2000–20070.040.170(4.5)ŷi=α̂(ti)-0.27xi1-0.20xi2,i=1,…,40.
The residual mean square of PLS + HT is 2.24×10-2 and that of PLS is 2.47×10-2. The selected number of basis functions is six with one excluded basis function and the spread parameter s is estimated as 0.30. Table 11 shows that PLS + HT selects only log(FLPi) 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 representsy; the dotted lines are the estimated curves of y; 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.



---

*Source: 487194-2009-09-22.xml*