2.2.1 Penalized splines in the cross-sectional context
We start our considerations by discussing penalized spline specifications for cross-sectional data. Consider the additive regression model
yi=β0+
H
X
h=1
fh(xhi) +ui, ui∼N 0, σu2
, i= 1, . . . , n, (2.1)
where yi is the response variable of interest, β0 is an overall intercept term, f1(x1i), . . . , fp(xHi) are smooth functions representing the potentially nonlinear effects of H deterministic covariates andui are independent and identically distributed normal error terms with variance σu2.2 To approximate the potentially nonlinear effects fh, we use the weighted sum ofdh B-spline basis functions, Bh1, . . . , Bhdh, such that
fh(xhi)≈
dh
X
j=1
Bhj(xhi)βhj=zTh(xhi)βh, (2.2) whereβhis adh-dimensonal column vector of basis coefficients andzh(xhi) is thedh-dimensonal column vector containing the evaluations of the basis functions at the observed covariate value xhi. Thereby, the amount of basis functions and coefficients dh is steered by the number of knots which divide the domain of the covariate. The bias introduced by the spline representation of a smooth function converges to zero with growing number of knots, see Claeskens et al. (2009) for details. We assume this bias to be negligible by using sufficiently many knots and subsequently postulate equality between an arbitrary smooth function and its spline representation, which leads to the expression
fh(xh) =Zhβh
in compact matrix notation, where
Zh=
Bh1(xh1) . . . Bhdh(xh1) ... . .. ... Bh1(xhn) . . . Bhdh(xhn)
is a design matrix of dimension n×dh assumed to be of full rank. In order to avoid an overfit to the data, a matrix Kh penalizing to much variability of adjacent coefficients in the coefficient vectorβh is assigned to each smooth function resulting in the penalized least squares criterion
y−β01n−
H
X
h=1
Zhβh
!T
y−β01n−
H
X
h=1
Zhβh
! +
H
X
h=1
λhβThKhβh, (2.3)
whereydenotes then-dimensional response vector,1nis ann-dimensonal column vector of ones andλhis a smoothing parameter determining the impact of the penalty on the minimization criterion. Thedh×dh -dimensional matrix Kh of first-order differences, that is penalizing differences of directly contiguous
2For notational simplicity, we refrain from adding stochastic covariates and covariates with strictly parametric effects.
However, as can be seen in Section 2.5, semiparametric partially linear models can also be handled easily within our framework.
coefficients, has the form
Kh=
1 −1
−1 2 −1 . .. . .. . ..
−1 2 −1
−1 1
.
Difference matrices of higher orders can be easily constructed. Details on such penalties for B-spline functions can be found, for example, in Eilers and Marx (1996).
Let nowxhbe an arbitrary value on the domain ofxh.Defining the smoothing matrixLh(xh) as Lh(xh) = (In−S−h)Zh
ZTh(In−S−h)Zh+λhKh−1zTh(xh) (2.4) with In denoting the identity matrix of dimension n×n, S−h = Z−h(ZT−hZ−h +λ−hK−h)−1ZT−h with λ−h = (λ1, . . . , λh−1, λh+1, . . . , λH), Z−h = (Z1, . . . ,Zh−1,Zh+1, . . . ,ZH), K−h= (K1, . . . ,Kh−1,Kh+1, . . . ,KH), andzh(xh) defined as in (2.2), the estimator of eachfh(xh) can be written as
fˆh(xh) =LTh(xh)y.
It follows that Varh
fˆh(xh)i
= Var
LTh(xh)y
=LTh(xh)Var(y)Lh(xh) =LTh(xh)σu2InLh(xh) (2.5) holds for homoscedastic and independent errors.
Naturally, the smoothing parametersλhare unknown. One way to estimate them is to exploit the mixed model representation of penalized splines. In particular, it is always possible to rewrite
Zhβh=Zh(Fhfαhf+Fhrαhr) =Xhfαhf+Xhrαhr,
where (Fhf,Fhr) is of full rank,FThfFhr =FThrFhf =FThfKhFhf = 0 and FThuKhFhu =Idh−q with the difference penalty orderq.3 It follows thatXhfis of dimensionn×qandXhris of dimensionn×(dh−q).
Then, αhf containsq fixed coefficients andαhr is a vector of (dh−q) virtually penalized random coef-ficients which are assumed to be mutually independent and normally distributed with constant variance σ2hr and independent from the errorsui. In this mixed model formulation, we obtain estimates both for the coefficients (fixed and random) and smoothing parameters by a single (restricted) maximum likelihood estimation. The smoothing parameter estimators ˆλh = σˆσˆ22u
hr
are then given as ratios of two (estimated) variances. For details we refer the reader to Ruppert and Wand (2003) or Fahrmeir et al. (2013). In Section 2.3 we will make use of the the mixed model formulation to construct confidence bands.
Asymptotic properties of the penalized spline estimator have been studied, among others, by Claeskens et al. (2009), Kauermann et al. (2009), Wang et al. (2011), Antoniadis et al. (2012), Yoshida and Naito (2012) and Yoshida and Naito (2014). Under mild conditions, consistency of the estimator is shown by Claeskens et al. (2009) for a univariate model with i.i.d. errors. Antoniadis et al. (2012), Yoshida and Naito (2012) and Yoshida and Naito (2014) discuss the asymptotic properties for additive models and derive consistency within different frameworks, always including the case of i.i.d. errors. As we will show, our models can be transformed in such a way that they fit into the class of additive models with i.i.d.
errors.
3One way to obtain such a decomposition is described in Wood (2006, pp. 316-317).
It should be noted that each row in the initial design matrix Zh (i.e., before applying the mixed model reformulation) for each covariate sums up to one, that is,Pd
j=1Bhj(xhi) = 1∀i= 1, . . . , n.Obviously, this leads to an identification problem in an additive model with an intercept or multiple smooth components.
This issue can be solved by imposing a centering constraint on each functionfhsuch that
n
X
i=1
fh(xhi) =
n
X
i=1
zTh(xhi)βh=0
holds for all h= 1, . . . , H. Following the ideas of Wood (2006, pp. 167-168), this can be achieved by constructing appropriate matrices Wh of dimension dh×(dh−1) with orthogonal columns, leading to a reparameterized model with design matrices ˜Zh = ZhWh and penalty matrices ˜Kh = WThKhWh. If the mixed model framework is used to determine the smoothing parameters as described above, the reparameterizing procedure to ensure identifiability is done before the mixed model reformulation of the model.
2.2.2 Penalized splines for panel data: A first-difference estimator
In comparison to cross-sectional data leading to model (2.1) introduced in the previous section, we now consider individuals (e.g., persons) observed at T consecutive points of time.4 We therefore consider an additive panel data model
yit=γi+
H
X
h=1
fh(xhit) +uit, i= 1, . . . , N, t= 1, . . . , T, (2.6)
where uit are assumed to be independent and normally distributed errors with constant variance and γi are individual-specific, time-invariant fixed effects allowed to be correlated with other covariates. As model (2.6) holds for each point of time, we obtain
yi,t−1=γi+
H
X
h=1
fh(xhi,t−1) +ui,t−1 (2.7)
for a one period time lag. To cancel out the individual-specific effects γi, we subtract (2.7) from (2.6) and obtain
∆yit=yit−yi,t−1=γi−γi+
H
X
h=1
[fh(xhit)−fh(xhi,t−1)] +uit−ui,t−1 (2.8)
=
H
X
h=1
dh
X
j=1
Bhj(xhit)βhj−
dh
X
j=1
Bhj(xhi,t−1)βhj
+ ∆uit
=
H
X
h=1
[zh(xhit)−zh(xhi,t−1)]Tβh+ ∆uit
=
H
X
h=1
[∆zh(xhit)]Tβh+ ∆uit,
where equation (2.2) is used for the second and third equality and ∆ denotes the first-difference operator over time. Note that onlyT−1 observations per individual are retained after differencing. Accordingly,
4The only reason to refrain from incorporating different observation horizons between persons is notational convenience.
As can be seen in Section 2.4 and Section 2.5, unbalanced panels can be handled without any difficulties in our framework.
as theN T ×dh-dimensional design matrixZhof the evaluated basis functions is now given by
Zh=
Bh1(xh11) . . . Bhdh(xh11) ... . .. ... Bh1(xh1T) . . . Bhdh(xh1T)
... . .. ... Bh1(xhN1) . . . Bhdh(xhN1)
... . .. ... Bh1(xhN T) . . . Bhdh(xhN T)
, (2.9)
we obtain
∆y=
H
X
h=1
∆Zhβh+ ∆u (2.10)
in compact matrix notation, where ∆y= (y12−y11, . . . , y1T−y1,T−1, . . . , yN2−yN1, . . . , yN T−yN,T−1)T is aN(T−1)-dimensional column vector, ∆uis defined analogously and theN(T−1)×dh-dimensional matrix ∆Zh is obtained by building the difference between matrixZh in (2.9) and its one period lagged counterpart:
∆Zh=
Bh1(xh12) . . . Bhdh(xh12) ... . .. ... Bh1(xh1T) . . . Bhdh(xh1T)
... . .. ... Bh1(xhN2) . . . Bhdh(xhN2)
. ..
Bh1(xhN T) . . . Bhdh(xhN T)
−
Bh1(xh11) . . . Bhdh(xh11) ... . .. ... Bh1(xh1,T−1) . . . Bhdh(xh1,T−1)
... . .. ... Bh1(xhN1) . . . Bhdh(xhN1)
... . .. ... Bh1(xhN,T−1) . . . Bhdh(xhN,T−1)
.
Additionally taking into account penalization, a first-difference penalized spline estimator for all βh can be obtained by minimizing the penalized least squares criterion
"
∆y−
H
X
h=1
(∆Zh)βh
#T"
∆y−
H
X
h=1
(∆Zh)βh
# +
H
X
h=1
λhβhT
Khβh. (2.11)
Since the smoothing parameters are unknown, one can again exploit the mixed model representation and using (restricted) maximum likelihood estimation as discussed in the previous subsection. Note that the framework is similar to the cross-sectional data case since only the vector of the dependent variable and the design matrices differ between the equations (2.3) and (2.11). The major difference in comparison to cross-sectional data is the problem of autocorrelated errors which are often encountered in panel data contexts. Krivobokova and Kauermann (2007) show that the restricted maximum likelihood based estimation of a smoothing parameter is robust to modest forms of autocorrelation. Moreover, further adjustments on the design matrices and the dependent variable for addressing serial correlation are also possible, see Section 2.3 for an elaboration of this topic.
We briefly have to refer to the identification problem in case of multiple smooth model components: Our aim is to estimate the functions fh, h= 1, . . . , H. Hence, model (2.6) should be identified such that
N
X
i=1 T
X
t=1
fh(xhit)≈Zhβh= 0
holds for allh= 1, . . . , H.Therefore, we rewrite the design matrices of the evaluated basis function given in (2.9) and the penalty matrices such that ˜Zh=ZhWhand ˜Kh=WThKhWh,proceeding as described in the previous subsection. Furthermore, the identification restriction also implies that a one period lagged design matrix is then constructed directly from ˜Zh by taking its one-period-lagged rows. After building the difference between each ˜Zhand its respective lagged counterpart, the resulting matrices ∆ ˜Zh
and the penalty matrices ˜Kh are plugged into (2.11) to obtain estimators forβhand thus forfh. Another common approach in fixed effects panel data models is time-demeaning, that is, removing the individual-specific effects γi by building the mean over time for each individual in equation (2.6) and subtracting the resulting equation from (2.6). Using the information above, this variant is straightforward to derive.