Chapter 6

Kernel estimators in time series analysis

6.1 Non-parametric estimation

Non-parametric methods are widely used in non-linear model building. Such methods are particularly useful if no prior information about the structure is available, since the estimation procedure is free of parameters and model structure (apart from a smoothing constant).

In this chapter we concentrate on kernel estimation. Some other non-parametric methods are based on splines, k nearest neighbour (k-NN) or orthogonal series smoothing. Each method has a specific weighting sequence {Ws(x); s = 1, ..., N}. These weighting sequences are related to each other and it can be argued (Härdle [1990]) that one of the simplest ways of computing a weighting sequence is kernel smoothing. For more information about non-parametric methods see Robinson [1983], Härdle [1990], Silverman [1986], Ruppert et al. [2003] or Hastie et al. [2011].

6.2 Kernel estimators for time series

6.2.1 Introduction

This chapter considers kernel estimation in general and its use in time series analysis. The non-parametric estimators are of particular relevance for non-Gaussian or non-linear time series. From kernel estimates of probability density functions Gaussianity can be verified or the nature of non-Gaussianity can be discovered. Thus non-parametric methods provide information that complements that given, for instance, by higher-order spectral analysis.

Non-parametric estimates of the conditional expectation can be used to detect non-linear structures and to identify a family of relevant time series models. Consider for instance the time series generated by the non-linear model:

$Y_t=g(Y_{t-1},...,Y_{t-p})+h(Y_{t-1},...,Y_{t-p}){\epsilon }_t.(6.1)$

For this model

$\begin{array}{cc}E[Y_t|Y_{t-1},...]=g(Y_{t-1},...,Y_{t-p})& (6.2)\\ \mathrm{Var}[Y_t|Y_{t-1},...]=h^2(Y_{t-1},...,Y_{t-p})& (6.3)\end{array}$

where it is implicitly assumed that the variance of the white noise process εt is unity in order to avoid identifiability issues.

Thus the non-parametric methods can be used to estimate, e.g., g(y) = E[Yt + 1 |Yt = y] in the simple case p = 1. Using this non-parametric estimate a relevant parametrization of g(·) can be suggested.

Furthermore non-parametric estimates can be used to provide predictors, i.e., estimators of future values of the time series. Use of non-parametric estimators for getting insight in the time-varying behaviour and a complex dependency on exogenous variables will also be considered.

In this section we first introduce the most frequently considered non-parametric estimators for time series, the kernel estimators. Next the kernel estimators are described in more detail. Finally, some examples are given on how to use kernel estimation in time series analysis.

6.2.2 Kernel estimator

Let {Yt; t = 0, ±1, ...} be a strictly stationary process. A realization is given as a time series {yt; t = 1, ..., N}, and that single realization of the process is the basis for inference about the process.

Introduce Zt = (Yt + j1, ..., Yt + jn) and $Y_t^{*}=(Y_{t-h_1},...,Y_{t+h_m}).$ . In order to simplify the notation we consider the case$Y_t^{*}=(Y_t)$ in the rest of this section. It is obvious how to generalize the equations.

The basic estimated quantity is

$E[G(Z_t)|Y_t=y]f_Y(y)(6.4)$

where fY is the probability density function (pdf) of Yt and G is a known function.

The estimator of (6.4) is

$[G(Z_t);y]={(N^{\prime }h)}^{-1}{\displaystyle \sum _{t=1}^{N^{\prime }}G(Z_t)k\left(\frac{y-Y_t}{h}\right)}(6.5)$

where k(u) is a real bounded function such that ∫ k(u)du = 1, and h is a real number. Here, N′ is some number which corrects for the fact that not all N observations can be used in the sum, typically N′ = N − jn, and k(u) is the kernel function and h is the bandwidth.

Let us consider two examples of using Equation (6.5). In the first example it is illustrated that Equation (6.5) is a reasonable estimator.

Example 6.1

(Non-parametric estimation of a pdf). The pdf of Yt at y is estimated by

${\widehat{f}}_Y(y)=[1;y]=\frac{1}{Nh}{\displaystyle \sum _{t=1}^{N}k\left(\frac{y+Y_t}{h}\right).}(6.6)$

Assume that Y has pdf f (y), then

$f(y)=\underset{h\to 0}{\lim }\frac{1}{2h}\mathbb{P}(y-h<Y\le y+h)(6.7)$

when ℙ is the distribution function at Y. An estimator of f(y) is

$\widehat{f}(y)=\frac{1}{2hn}[Numberofobservationsin(y-h;y+h)].(6.8)$

Define the rectangular kernel

$w(u)=\{\begin{array}{cc}\frac{1}{2}& if|u|<1\\ 0& otherwise.\end{array}(6.9)$

Then the estimator (6.8) can be written as

$\widehat{f}(y)=\frac{1}{N}{\displaystyle \sum _{t=1}^N\frac{1}{h}w\left(\frac{y-Y_t}{h}\right)}.$

Now consider the general kernel (∫ k(u)du = 1). Then the estimator is

$\widehat{f}(y)=\frac{1}{Nh}{\displaystyle \sum _{t=1}^N\left(\frac{y-Y_t}{h}\right)}(6.10)$

which clearly is equal to Equation (6.6).

Example 6.2

The estimator of the conditional expectation of G(Zt), given Yt = y, is

$\begin{array}{ll}\begin{array}{lll}\widehat{E}[G(Z_t)|Y_t=y]& =& \frac{[G(Z_t);y]}{1;y}\\ & =& \frac{\frac{1}{N^{\prime }}\sum G(Z_t)k\left(\frac{y-Y_t}{h}\right)}{\frac{1}{N}{\sum }_k\left(\frac{y-Y_t}{h}\right)}\end{array}& (6.11)\end{array}$

where [1; y] was found in the previous example.

Of special interest is the case G(Zt) = Yt + 1, where (6.11) estimates E[Yt + 1|Yt = y].

6.2.3 Central limit theorems

Central limit theorems can be established under various weak dependence conditions on the process {Yt}. Let ${ℱ}_u^v$ be the σ-field of events generated by Yt; u ≤ t ≤ v. Then introduce the coefficient

${\alpha }_j=\underset{A\in {\sigma }_{-\infty }^t,B\in {\sigma }_{t+j}^{\infty }}{\mathit{\text{sup}} }|P(A\cap B)-P(A)P(B)|:j>0.(6.12)$

Definition 6.1

(Strong mixing condition). The process Xt is said to be strongly mixing if αj, → 0 as j → ∞.

The strong mixing condition has been used frequently in the asymptotic theory of estimators for time series; but the condition is often very difficult to check.

Central limit theorems for non-parametric estimators of pdfs of a process {Yt}, as well as of conditional pdfs and conditional expectations at continuity points, are given by Robinson [1983] under the strong mixing condition.

6.3 Kernel estimation for regression

This section describes in more detail the kernel estimation technique. Since the method is useful for non-parametric regression in general, a notation will be used which is related to non-parametric regression; however, the problem is highly related to the estimator in (6.11), and thus also useful in time series analysis.

6.3.1 Estimator for regression

The kernel based non-parametric regression was first proposed by Nadaraya [1964].

Assume that the theoretical relation is

$\begin{array}{cc}Y=g(X^{(1)}),...,X^{(q)}+\epsilon & (6.13)\end{array}$

where ε is a white noise, and g is a continuous function.

Given n observations

$O=\{(X_1^{(1)},...,X_1^{(q)},Y_1),...,(X_n^{(1)},...,X_n^{(q)},Y_n)\}$

the goal is now to estimate the function g.

If q = 1, then the kernel estimator for g, given the observations, is

$\widehat{g}(x)=\frac{\frac{1}{n}{\displaystyle {\sum }_{s=1}^n{Y}_{s}k\{h^{-1}(x-X_s)\}}}{\frac{1}{n}{\displaystyle {\sum }_{s=1}^{n}k\{h^{-1}(x-X_s)\}}}.(6.14)$

Note that if we compare with Equation (6.5) then G(Zt) = Zt = Yt and the independent variable is now X.

The function k is the kernel, and h is the bandwidth of the kernel, where h → 0 as n → ∞. Various candidate functions for the kernel have been proposed, but it is shown in Härdle [1990] that the parabolic kernel with bounded support, called the Epanechnikov kernel, minimizes the mean square error among all kernels in which the bandwidth is optimally chosen.

The Epanechnikov kernel is given by

$k_h^{Epa}(u)=\frac{3}{4h}\left(1-\frac{u^2}{h^2}\right)I_{\{|u|\le h\}}.(6.15)$

Other possibilities are the rectangular, triangular or the Gaussian kernel given by:

$k_h^{Gau}(u)=\frac{1}{h\sqrt{2\pi }}\exp \left(-\frac{u^2}{2h^2}\right).(6.16)$

These kernels are depicted in Figure 6.1. The Epanechnikov kernel is shown for h = 1, and the Gaussian kernel for h = 0.45.

If we allow q > 1 in (6.14) a kernel of order q has to be used, i.e., a real bounded function kq(x1,..., xq), where ∫ k(u)du = 1.

By a generalization of the bandwidth h to the quadratic matrix h with the dimension q × q the more general kernel estimator for g(X(1),..., X(q)) is

$\widehat{g}(x)=\frac{\frac{1}{n}{\displaystyle {\sum }_{s=1}^n{Y}_s{k}_q[h^{-1}(x-X_s)]}}{\frac{1}{n}{\displaystyle {\sum }_{s=1}^n{k}_q[h^{-1}(x-X_s)]}}(6.17)$

where x = (x1,...,xq) and $X_s=(X_1^{(1)},...,X_1^{(q)})$. This is called the Nadaraya-Watson estimator.

6.3.2 Product kernel

In practice product kernels are used, i.e.,

$k_q(h^{-1}x){\displaystyle \prod _{i=1}^{q}k(x_i/h_i)}(6.18)$

where k is a kernel of order 1.

Assuming that hi = h, i = 1,..., q the following generalization of the one-dimensional case in (6.14) is obtained

$\widehat{g}(x_1,...,x_q)=\frac{\frac{1}{n}{\displaystyle {\sum }_{s=1}^n{Y}_s{\Pi }_{i=1}^{q}k\{h^{-1}(x_i-X_s^{(i)})\}}}{\frac{1}{n}{\displaystyle {\sum }_{s=1}^n{\text{Π}}_{i=1}^{q}k\{h^{-1}(x_i-X_s^{(i)})\}}}.(6.19)$

6.3.3 Non-parametric estimation of the pdf

An estimate of the probability density function (pdf) is Robinson [1983]

$\widehat{f}(x_{i_1},...,x_{i_d})=\frac{1}{n-i_d+i_1}{\displaystyle \sum _{s=i_d+1}^{n+i_1}{\displaystyle \prod _{j=1}^d{h}^{-1}k\{h^{-1}(x_{i_j}-X_{s-i_j})\}.(6.20)}}$

6.3.4 Non-parametric LS

If we define the weight Ws (x) by

$W_s(x)=\frac{{\prod }_{i=1}^{q}k\{h^{-1}(x_i-X_s^{(i)})\}}{\frac{1}{n}{\sum }_{s=1}^n{\prod }_{i=1}^{q}k\{h^{-1}(x_i-X_s^{(i)})\}}(6.21)$

then it is seen that Equation (6.19) corresponds to the local average:

$\widehat{g}(x_1,...,x_q)=\frac{1}{n}{\displaystyle \sum _{s=1}^n{W}_s(x)Y_s}.(6.22)$

The estimate is actually a non-parametric least squares estimate at the point x. This is recognized from the fact that the solution to the least squares problem

$\arg \underset{\theta }{ \min }\frac{1}{n}{\displaystyle \sum _{s=1}^n{W}_s(x){(Y_s-\theta )}^2(6.23)}$

is given by

${\widehat{\theta }}_{LS}(x)=\frac{{\sum }_{s=1}^n{W}_s(x)Y_s}{{\sum }_{s=1}^n{W}_s(x)}.(6.24)$

This shows that at each x, the estimate ĝ is a scaled weighted LS location estimate, i.e.,

$\widehat{g}(X)={\widehat{\theta }}_{LS}\frac{1}{n}{\displaystyle \sum _{s=1}^n{W}_s(x)}(6.25)$

6.3.5 Bandwidth

The bandwidth h determines the smoothness of ĝ. In analogy with smoothing in spectrum analysis:

• If h is small the variance is large but the bias is small.
• If h is large the variance is small but the bias is large.

The limits provide some insight:

• As h → ∞ it is seen that $\widehat{g}(x_1,...,x_q)=\overline{Y}$ as all data are included and given equal weights.
• As h → 0 it is seen that ĝ(x1,...,xq) = Yi for $(x_1,...,x_q)=(X_i^{(1)},...,X_i^{(q)})$ and otherwise undefined (or possibly 0 depending on how ratios of zeros are defined).

6.3.6 Selection of bandwidth — cross validation

The ultimate goal for the selection of h is to minimize the mean square error (MSE)(see Härdle [1990])

$\begin{array}{cc}\text{MSE(}h\text{)=}\frac{1}{\text{n}}{\displaystyle \sum _{i=1}^n{[\widehat{g}(X_i^{(1)},...,X_i^{(q)})-g(X_i^{(1)},...,X_i^{(q)})]}^2\times }& \text{(6}\text{.26)}\\ w(X_i^{(1)},...,X_i^{(q)})& \text{(6}\text{.27)}\end{array}$

where w(...) is a weight function which screens off some of the extreme observations, and g is the unknown function. A solution is the “plug-in” method, where Yi is used as an estimate of $g(X_i^{(1)},...,X_i^{(q)})$ in (6.26).

Hence, the criterion is

$\begin{array}{cc}\widehat{\text{MSE}}(h)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{[\widehat{g}(X_i^{(1)},\mathrm{...},X_i^{(q)})-Y_i]}^{2}w(X_i^{(1)},\mathrm{...},X_i^{(q)})}& (6.28)\end{array}$

but it is clear that $\widehat{\text{MSE}}(h)\to 0$ for h → 0, since $y_i-\widehat{g}(x_i^{(1)},...,x_i^{(q)})\to 0$ when h → 0.

It is clear that a modification is needed! In the “leave one out” estimator the idea is to avoid that ($X_i^{(1)},...,X_i^{(q)}$) is used in the estimate for Yi. For every data ($X_i^{(1)},...,X_i^{(q)},Y_i$) we define an estimator ${\widehat{g}}_{(i)}$ for Yi based on all data except ($X_i^{(1)},...,X_i^{(q)}$). The n estimators ${\widehat{g}}_{(1)},...,{\widehat{g}}_{(n)}$ (called the “leave one out” estimators) are written

${\widehat{g}}_{(j)}=(x_1,...,x_q)=\frac{\frac{1}{n-1}{\Sigma }_{s\ne j}{Y}_s{\text{Π}}_{i=1}^{q}k\{h^{-1}(x_i-X_s^{(i)})\}}{\frac{1}{n-1}{\Sigma }_{s\ne j}{\text{Π}}_{i=1}^{q}k\{h^{-1}(x_i-X_s^{(i)})\}}.(6.29)$

Now the cross-validation criterion using the “leave one out” estimates ${\widehat{g}}_{(1)},...,{\widehat{g}}_{(n)}$ is

$CV(h)=\frac{1}{n}{\displaystyle \sum _{i=1}^n}[I{[{\widehat{g}}_{(i)}(X_i^{(1)},...,X_i^{(q)})-Y_i]}^{2}w(X_i^{(1)},...,X_i^{(q)}).(6.30)$

It can be shown that under weak assumptions the estimate of the bandwidth $\widehat{h}$ that is obtained by minimizing the cross-validation criterion is asymptotic optimal, i.e., it minimizes (6.26) (Härdle [1990]). An example of using the CV procedure is shown in Figure 6.7.

6.3.7 Variance of the non-parametric estimates

To assess the variance of the curve estimate at point x Härdle [1990] proposes the pointwise estimator given by

${\widehat{\sigma }}^2(x)=\frac{1}{n}{\displaystyle \sum _{s=1}^n}{W}_s(x){(\widehat{g}(X_s)-Y_s)}^2(6.31)$

where the weights are given as shown previously by

$W_s(x)=\frac{{\text{Π}}_{i=1}^{q}k\{h^{-1}(x_i-X_s^{(i)})\}}{\frac{1}{n}{\Sigma }_{s=1}^n{\text{Π}}_{i=1}^{q}k\{h^{-1}(x_i-X_s^{(i)})\}}.(6.32)$

Remembering the WLS interpretation this estimate seems reasonable.

6.4 Applications of kernel estimators

6.4.1 Non-parametric estimation of the conditional mean and variance

Assume a realization of a stochastic process {X1,..., Xn}.

Goal: Use the realization to estimate the functions g(·) and h(·) in the model

$X_t=g(X_{t-1},...,X_{t-p})+h(X_{t-1},...,X_{t-q}){\epsilon }_t.(6.33)$

For the conditional mean and the conditional variance we shall use the notation

$\begin{array}{cc}M(x_{i_1},...,x_{i_d})=E[X_t|i{X}_{t-i_1}=x_{i_1},...,X_{t-i_d}=x_{i_d}],& (6.34)\\ V(x_{i_1},...,x_{i_d})=\mathrm{Var}[X_t|X_{t-i_1}=x_{i_1},...,X_{t-i_d}=x_{i_d}].& (6.35)\end{array}$

One solution is to use the kernel estimator (other possibilities are splines, nearest neighbour or neural network estimates).

Using a product kernel we obtain

$\widehat{M}(x_{x_1},...,x_{i_d})=\frac{\frac{1}{n-i_d}{\Sigma }_{s=i_d+1}^n{X}_s{\text{Π}}_{j=1}^{d}k\{h^{-1}(x_{i_j}-X_{s-i_j})\}}{\frac{1}{n-i_d+i_1}{\Sigma }_{s=i_d+1}^{n+i_1}{\text{Π}}_{j=1}^{d}k\{h^{-1}(x_{i_j}-X_{s-i_j})\}}.(6.36)$

Assuming that E[Xt] = 0 the estimator for V(·) is

$\widehat{V}(x_{i_1},...,x_{i_d})=\frac{\frac{1}{n-i_d}{\Sigma }_{s=i_d+1}^n{X}_s^2{\text{Π}}_{j=1}^{d}k\{h^{-1}(x_{i_j}-X_{s-i_j})\}}{\frac{1}{n-i_d+i_1}{\Sigma }_{s=i_d+1}^{n+i_1}{\text{Π}}_{j=1}^{d}k\{h^{-1}(x_{i_j}-X_{s-i_j})\}}\text{.}(6.37)$

If E[Xt] ≠ 0 it is clear that the above estimator is changed to

$\begin{array}{cc}\begin{array}{c}\widehat{V}(x_{i_1},...,x_{i_d})=\frac{\frac{1}{n-i_d}{\Sigma }_{s=i_d+1}^n{X}_s^2{\text{Π}}_{j=1}^{d}k\{h^{-1}(x_{i_j}-X_{s-i_j})\}}{\frac{1}{n-i_d+i_1}{\Sigma }_{s=i_d+1}^{n+i_1}{\text{Π}}_{j=1}^{d}k\{h^{-1}(x_{i_j}-X_{s-i_j})\}}\\ -{\widehat{M}}^2(x_{i_1},...,x_{i_d}).\end{array}& (6.38)\end{array}$

Example 6.3.

Simulations of 10 independent time series have been generated from the following non-linear models:

$\begin{array}{llll}A:\hfill & X_t\hfill & =\hfill & \{\begin{array}{cc}-0.8X_{t-1}+\epsilon & x_{t-1}\ge 0\\ 0.8X_{t-1}+1.5+{\epsilon }_t,& x_{t-1}<0\end{array}\hfill \\ B:\hfill & X_t\hfill & =\hfill & \sqrt{1+X_{t-1}^2}{\epsilon }_t.\hfill \end{array}$

The first model is a SETAR model with non-linear conditional mean and constant conditional variance whereas the second model is an ARCH model with constant conditional mean and state dependent conditional variance.

A Gaussian kernel with h = 0.6 is used to estimate the conditional mean (Figure 6.2) and variance (Figure 6.3) in order to detect the nonlinearities. The theoretical conditional mean and variance is shown with "+”.

6.4.2 Non-parametric estimation of non-stationarity — an example

Non-parametric methods can be applied to identify the structure of the existing relationships leading to proposals for parametric model classes. An example of identifying the diurnal dependence in a non-stationary time series is given in the following.

Measurements of heat supply from 16 terrace houses in Kulladal, a suburb of Malmö in Sweden, are used to estimate the heat load as a function of the time of the day.

The heat supply was measured every 15 minutes for 27 days. The transport delay between the consumers and the measuring instruments is not more than a couple of minutes.

It is assumed that the heat supply can be related to the time of day, and that there is no difference between the days for the considered period. Then the regression curve is a function only of the time of day, and it is this functional relationship, which will be considered.

Figure 6.6 shows the Epanechnikov curve estimate calculated for a spectrum of bandwidths. It is clear that the characteristics of the non-smoothed averages gradually disappear, when the bandwidth is increased.

When the cross-validation function is brought into action for bandwidth selection on the original data we get the result shown in Figure 6.7. Obviously the minimum of the CV-function is obtained for a bandwidth close to 1.3.

6.4.3 Non-parametric estimation of dependence on external variables — an example

This example illustrates the use of kernel estimators for exploring the (static) dependence on external variables. Data were collected in the district heating system in Esbjerg, Denmark, during the period August 14 to December 10, 1989.

The interest in this example is the heat load in a district heating system, when the two most influential explanatory variables in district heating systems, i.e., ambient air temperature and supply temperature, are included in the study.

For the estimation of the regression surface the kernel method is applied using the Epanechnikov kernel (6.15). A product kernel is used, i.e., the power load estimate, p, at time t, at ambient air temperature a and, for supply, temperature s is estimated as

$\widehat{p}{h}_t,h_a,h_s(t,a,s)=\frac{{\Sigma }_{i=1009}^{3843}p(i)K_{h_t}(t-t(i))K_{h_a}(a-a(i))K_{h_s}(s-s(i))}{{\Sigma }_{i=1009}^{3843}{K}_{h_t}(t-t(i))K_{h_a}(a-a(i))K_{h_s}(s-s(i))}(6.39)$

where Kh(u) = k(u/h). The bandwidths were chosen using a combination of Cross-Validation and visual inspection.

6.4.4 Non-parametric GARCH models

The very large number of parametric GARCH models makes model selection non-trivial (Bollerslev [2008] for a recent overview of different models). Non-parametric GARCH models were introduced by Bühlmann and McNeil [2002] as an alternative to the non-parametric models. Their model is defined as

$\begin{array}{llll}{X}_t\hfill & =\hfill & {\sigma }_t{Z}_t\hfill & (6.40)\hfill \\ {\sigma }_t^2\hfill & =\hfill & f(X_{t-1},{\sigma }_{t-1}^2)\hfill & (6.41)\hfill \end{array}$

where {Zt} is an iid sequence of zero mean, unit variance random variables. Their basic idea is to start from a parametric GARCH model, and iteratively improve the estimate of the volatility function. Their algorithm is given by:

1. Estimate a simple GARCH model, recovering an approximation of the volatility {σt,0: 1 ≤ t ≤n}.
2. Regress $X_t^2$ against Xt−T and ${\sigma }_{t-1,m-1}^2$ to estimate the unknown function f(·,·), denoted ${\widehat{f}}_m(\cdot ,\cdot )$, for m = 1,...M.
3. Calculate ${\widehat{\sigma }}_{t,m}^2$ from the estimated function ${\widehat{f}}_m(\cdot ,\cdot )$.

The algorithm iterates between step 2 and 3 until the estimated function $\widehat{f}(\cdot ,\cdot )$ has converged.

The method was evaluated in a simulation, where the volatility is given by

$f(x,{\sigma }^2)=5+0.2x^2+(0.751_{\{x>0\}}+0.1\cdot 1_{\{x<0\}}){\sigma }^2.(6.42)$

It was shown in Bühlmann and McNeil [2002] that the algorithm typically converges in just a few iterations, but they also argue that it can be worthwhile to iterate a few extra times, and use the average over the estimated volatility surfaces as the final volatility forecast.

The non-parametric GARCH model was successfully used to forecast crude oil price return volatility in Hou and Suardi [2012]. This is an excellent application of the model, as the price dynamics of many commodities often are more complex than those of, say, equities.

6.5 Notes

The reader is encouraged to dig into Härdle [1990] for a more complete overview of non-parametric methods and Robinson [1983] for details on how this can be applied to dependent data.

Alternatively, semiparametric methods are nowadays an options (Ruppert et al. [2003]). The explanatory properties of semiparametric methods are similar to non-parametric methods (Härdle [2004], Hastie et al. [2009]), but there are times when these are easier to apply to data.