Second Moment/Order Approximations by Kernel Smoothers with Application to Volatility Estimation

1. Introduction

Regression analysis can be considered as a set of methodologies for estimating the relationships between a dependent variable y (often called as ’output’ or ’response’) and a vector $\mathbf{x}$ of independent variables (often called ’input’). Essentially, the main goal in a regression problem is to obtain an approximation of the trend defined as the expected value of y given $\mathbf{x}$, i.e., $E\left[y\right|\mathbf{x}]$, that is the first (non-central) moment of the conditional density $p\left(y\right|\mathbf{x})$. We can assert that the most complete regression problem consists of approximating the whole conditional density $p\left(y\right|\mathbf{x})$, whereas the simplest task in a regression problem consists of estimating only the first moment $E\left[y\right|\mathbf{x}]={\int }_{\mathcal{Y}}yp\left(y\right|\mathbf{x})dy$. Intermediate scenarios appear in different applications, where other moments (higher than one) are of interest and hence are also approximated.

Volatility estimation (intended as local variance or local standard deviation) and quantile regression analysis are currently important tasks in statistics, machine learning and econometrics. The problem of the estimation of the volatility has a particular relevance in financial time series analysis, where the volatility represents the degree of variation of a trading price series over time, usually measured by the standard deviation of logarithmic returns. It is still an important active area of research [1,2,3,4]. Local linear regression methods (which can resemble a kernel smoother approach) have been already proposed for volatility estimation [5,6,7]. Quantile regression models study the relationship between an independent input vector $\mathbf{x}$ and some specific quantile or moment of the output variable y[8,9,10]. Therefore, in quantile regression analysis the goal is to estimate higher order moments of the response/output variable, given an input $\mathbf{x}$ [8,9,10,11,12,13]. It is important to remark that many advance regression methods consider a constant variance as initial assumption of the output given the input, $var\left[y\right|\mathbf{x}]={\sigma }_e^2$, i.e., it does vary with the input variable $\mathbf{x}$ [14,15]. In this work, we consider an extended approach (with respect to the assumptions used in classical regression methods) where $var\left[y\right|\mathbf{x}]=v(\mathbf{x})$, i.e., the local variance varies with the input $\mathbf{x}$.

More precisely, in this work, we propose two procedures for estimating the local variance in a regression problem and using a kernel smoother approach [15]-[16]. Note that the kernel smoother schemes contain several well-known techniques as special cases, such as the fixed radius nearest neighbours, as an example [17]. The resulting solution is non-parametric method, hence both complexity and flexibility of the solution grows with the number of data N. This ensures to have the adaquate flexibility to analyze the data.

The first proposed method is based on the Nadaraya-Watson derivation of a linear kernel smoother [14,15]. The second proposed approach draws inspiration from divisive normalization, a function grounded in the activity of brain neurons [18]. This function aims to standardize neuron activity by dividing it by the activity of neighboring neurons. It has demonstrated favorable statistical properties [19,20,21] and has been utilized in various applications [22,23].

The proposed methods can be applied to time series analysis and/or in more general regression problems where the input variables are multidimensional. Therefore, our approach can be implemented in spatial statistical modelling and any other inverse problems [10,24]. Furthermore, another advantage of the proposed scenario, is that the generalization for a multi-output scenarios can be easily designed.

From another point of view, this work gives another relevant contribution. The proposed schemes allow to perform uncertainty estimation with a generic kernel smoother technique. Indeed, a generic kernel smoother method is not generally supported by a generative probabilistic derivation that yields also an uncertainty estimation. The Gaussian processes (GPs) and relevance vector machines (RVMs) regression methods are relevant well-known and virtually unique exceptions [15,25]).

We have tested the proposed methods in five numerical experiments. Four of them involve the application of the different schemes to time series analyis, and the comparison with GARCH models which are considered benchmark techniques for volatility estimation in time series [26]. The last numerical experiment addresses a more general regression problem, with a multidimensional input variable $\mathbf{x}$ of dimension 122, considering a real (emo-soundscapes) database [27,28].

2. Approximating the Trend

Let us consider a set of N data pairs, such as ${\{{\mathbf{x}}_i,y_i\}}_{i=1}^N$, where ${\mathbf{x}}_i\in {\mathbb{R}}^D$, with $D\ge 1$, and $y_i\in \mathbb{R}$. First of all, we are interested in obtaining a regression function (a.k.a.,“local mean”– trend), i.e., removing the noise in the signal obtaining an estimator $\widehat{f}\left(\mathbf{x}\right)$ for all possible $\mathbf{x}$ in the input space. One possibility is to employ a linear kernel smoother. More specifically, we consider the Nadaraya-Watson estimator [14,15] that has the following form,

$$\begin{array}{cc}\hfill E\left[y\right|\mathbf{x}]\approx \widehat{f}\left(\mathbf{x}\right)& =\hfill \\ \hfill & =\sum _{n=1}^N\frac{h(\mathbf{x},{\mathbf{x}}_n)}{{\sum }_{j=1}^{N}h(\mathbf{x},{\mathbf{x}}_j)}{y}_n=\sum _{n=1}^N\phi (\mathbf{x},{\mathbf{x}}_n)y_n,\hfill \end{array}$$

(1)

where $\phi (\mathbf{x},{\mathbf{x}}_n)=\frac{h(\mathbf{x},{\mathbf{x}}_n)}{{\sum }_{j=1}^{N}h(\mathbf{x},{\mathbf{x}}_j)}$. Note that, by this definition, the nonlinear weights $\phi (\mathbf{x},{\mathbf{x}}_n)$ are normalized, i.e.,

$$\sum _{n=1}^N\phi (\mathbf{x},{\mathbf{x}}_n)=1,\forall \mathbf{x}.$$

(2)

As an example, we could consider

$$h(\mathbf{x},\mathbf{z})=h(\mathbf{x},\mathbf{z}|\lambda )=exp\left(-\frac{\left|\right|\mathbf{x}-{\mathbf{z}\left|\right|}^2}{\lambda }\right),$$

where $\lambda$ is a parameter that should be tuned. Clearly, we have also

$$\phi (\mathbf{x},\mathbf{z})=\phi (\mathbf{x},\mathbf{z}|\lambda ).$$

The form of this estimator above is quite general, and contains different other well-known methods as special cases. For instance, it contains the k-nearest neighbors algorithm (kNN) for regression as a specific case (to be more specific, the fixed-radius near neighbors algorithm). Indeed, with a specific choice of $h(\mathbf{x},{\mathbf{x}}_n|\lambda )$ (as a rectangular function), the expression above can represent the fixed-radius near neighbors estimator [15,17].

Remark 1. 

The resulting regression function is flexible non-parametric method. Both complexity and flexibility of the solution grows with the number of data N.

Remark 2. 

The input variables ${\mathbf{x}}_n\in {\mathbb{R}}^D$ are vectors ($D\ge 1$), in general. Therefore, the described methods have a much wider applicability than the techniques that can be employed only for time series (where the time index is a scalar number, $x=t\in \mathbb{R}$). Clearly, the methodologies described here can be also employed for analyzing time series. Moreover, even in a time series framework, we can obtain a prediction between two consecutive time instants. For instance, if we have a time series with daily data, with a kernel smoother we can have a prediction at each hour (or minute) within a specific day.

Learning 1 ($\lambda$) One possibility for tuning the parameters of the kernel functions is to use a cross-validation (CV) procedure. In this work, we have employed a leave-one-out cross-validation (LOO-CV) [14].

3. Variance Estimation Procedures

Let assume that we have already computed the trend (a.k.a., “local mean”), i.e., the regression function $\widehat{f}\left(\mathbf{x}\right)$. Here we present two methods to obtain an estimation of the local variance (or volatility) at each point $\mathbf{x}$, which is theoretically defined as

$$\begin{array}{cc}\hfill v\left(\mathbf{x}\right)=var\left[y\right|\mathbf{x}]& ={\int }_{\mathcal{Y}}{(y-E\left[y\right|\mathbf{x}])}^{2}p\left(y\right|\mathbf{x})dy,\hfill \\ \hfill & =E\left[y^2\right|\mathbf{x}]-{\left(E\left[y\right|\mathbf{x}]\right)}^2.\hfill \end{array}$$

(3)

Method 1(M1) If the weights $\phi (\mathbf{x},{\mathbf{x}}_n)$ are adequate for linearly combined the outputs $y_i$ and obtaining an proper approximation of $E\left[y\right|\mathbf{x}]$, one can extend this idea for approximating the second non-central moment $E\left[y^2\right|x]$ as

$$\begin{array}{c}\hfill E\left[y^2\right|\mathbf{x}]\approx \widehat{q}\left(\mathbf{x}\right)=\sum _{n=1}^N\phi (\mathbf{x},{\mathbf{x}}_n)y_n^2.\end{array}$$

(4)

hence

$$\begin{array}{c}\hfill \widehat{v}\left(\mathbf{x}\right)=\widehat{q}\left(\mathbf{x}\right)-{\left(\widehat{f}\left(\mathbf{x}\right)\right)}^2\approx var\left[y\right|\mathbf{x}],\end{array}$$

(5)

which is an estimator of theinstant variance. Note that$\widehat{q\left(\mathbf{x}\right)}$and$\widehat{f}\left(\mathbf{x}\right)$are obtained with the same weights

$$\phi (\mathbf{x},{\mathbf{x}}_n)=\phi (\mathbf{x},{\mathbf{x}}_n|\lambda ),$$

with the same value of$\lambda$(obtained using a LOO-CV procedure). Thus, a variant of this procedure consists into learning another value of${\lambda }_2$, i.e., obtaining other coefficients${\phi }_2(\mathbf{x},\mathbf{z})={\phi }_2(\mathbf{x},\mathbf{z}|{\lambda }_2)$, in the Eq. (4) considering the signal$y_n^2$(instead of$y_n$).

Method 2 (M2). Let us define the signal obtained as the estimated square errors

$$v_n={\left(y_n-\widehat{f}\left(\mathbf{x}\right)\right)}^2,n=1,...,N.$$

If$\widehat{f}\left(\mathbf{x}\right)\approx E\left[y\right|\mathbf{x}]$as assumed,$v_n$is a one-sample estimate of the variance at$\mathbf{x}$. Then, the goal is to approximate the trend of this new signal (i.e., new output)$v_n$,

$$\begin{array}{cc}\hfill \widehat{v}\left(\mathbf{x}\right)& =\sum _{n=1}^N{\phi }_2(\mathbf{x},{\mathbf{x}}_n)v_n,\hfill \end{array}$$

where we consider another parameter${\lambda }_2$, i.e.,

$${\phi }_2(\mathbf{x},\mathbf{z})={\phi }_2(\mathbf{x},\mathbf{z}|{\lambda }_2),$$

that is tuned again with LOO-CV but considering the new signal $v_n$. Note that$\widehat{v}\left(\mathbf{x}\right)$can be interpreted as estimation of the instant variance of the underlying signal. As alternative, also completely different kernel functions (as${\phi }_2$) can be applied (that differs for the functional form with respect to$\phi$, instead just for the choice of$\lambda$).

Remark 3. 

Again, as for estimating the trend, note that the two procedures above can be applied for multivariate inputs ${\mathbf{x}}_n\in {\mathbb{R}}^D$, and not just for scalar inputs (as in the time series).

Remark 4. 

If one divides $\mathbf{x}$ by $\widehat{v}\left(\mathbf{x}\right)$ we get a signal with uniform local variance, i.e. we have removed the (possible) heteroscedasticity. In [21], this procedure was used to define the relation kernels in the divisive normalization, and thus equalize locally the energy of neuron responses.

4. Extensions and Variants

The use of linear kernel smoothers as regression methods is not mandatory. Indeed, the ideas previously described can be employed even applying different regression methods. Below we provide some general steps in the same fashion of a pseudo-code, In order to clarify the application of possible different regression techniques:

• Choose a regression method. Obtain a trend function $\widehat{f}\left(\mathbf{x}\right)\approx E\left[y\right|\mathbf{x}]$ given the dataset ${\{{\mathbf{x}}_i,y_i\}}_{i=1}^N$.
• Choose one method for estimating the instance variance (between the two below):
  • M1. Choose a regression method (the same of the previous step or a different one). Considering the dataset the dataset ${\{{\mathbf{x}}_i,y_i^2\}}_{i=1}^N$ and obtain $\widehat{q}\left(\mathbf{x}\right)\approx E\left[y^2\right|\mathbf{x}]$. Then, compute $\widehat{v}\left(\mathbf{x}\right)=\widehat{q}\left(\mathbf{x}\right)-{\left(\widehat{f}\left(\mathbf{x}\right)\right)}^2$.
  • M2. Choose a regression method (the same of the previous step or a different one). Considering the dataset the dataset ${\{{\mathbf{x}}_i,v_i\}}_{i=1}^N$ where $v_n={(y_n-\widehat{f}\left(\mathbf{x}\right))}^2$, obtain the function $\widehat{v}\left(\mathbf{x}\right)$.

4.1. Multi-Output Scenario and other Extensions

In a multi-output framework, we have N data pairs, such as ${\{{\mathbf{x}}_i,{\mathbf{y}}_i\}}_{i=1}^N$, where ${\mathbf{x}}_i\in {\mathbb{R}}^D$, but in this case also the outputs are vectors ${\mathbf{y}}_i\in {\mathbb{R}}^{D_Y}$. Hence, for the local trend we have also a vectorial hidden function $\widehat{\mathbf{f}}\left(\mathbf{x}\right):{\mathbb{R}}^D\to {\mathbb{R}}^{D_Y}$ for all possible $\mathbf{x}$ in the input space. Then, we could easily write

$$\begin{array}{cc}\hfill \widehat{\mathbf{f}}\left(\mathbf{x}\right)& =\sum _{n=1}^N\phi (\mathbf{x},{\mathbf{x}}_n){\mathbf{y}}_n.\hfill \end{array}$$

(6)

Regarding the estimations of local variances follow the same procedures, defining the following vectorial quantities: $\widehat{\mathbf{q}}\left(\mathbf{x}\right)\approx E\left[{\mathbf{y}}^2\right|\mathbf{x}]$ for M1, $\widehat{\mathbf{v}}\left(\mathbf{x}\right)=\widehat{\mathbf{q}}\left(\mathbf{x}\right)-{\left(\widehat{\mathbf{f}}\left(\mathbf{x}\right)\right)}^2$ and ${\mathbf{v}}_n={({\mathbf{y}}_n-\widehat{\mathbf{f}}\left(\mathbf{x}\right))}^2$ for M2.

Furthermore, we could also consider different local parameters $\lambda$, for instance ${\lambda }_n$ with $n=1,...,N$, or more generally, with $\lambda =\lambda \left(\mathbf{x}\right)$. In this case, we would have different coefficient functions ${\phi }_n(\mathbf{x},{\mathbf{x}}_n)$ (one for each input ${\mathbf{x}}_n$), so that trend could be easily expressed as

$$\begin{array}{cc}\hfill \widehat{\mathbf{f}}\left(\mathbf{x}\right)& =\sum _{n=1}^N{\phi }_n(\mathbf{x},{\mathbf{x}}_n){\mathbf{y}}_n.\hfill \end{array}$$

(7)

4.2. Example of Alternative to the Kernel Smoothers for Time Series Analysis

Let us consider a time series framework, i.e., the input is a scalar time instance, $x=t\in \mathbb{R}$, then the dataset is formed by the following pairs ${\{t_i,y_i\}}_{i=1}^N$. Let also consider that the intervals, $t_i-t_{i-1}$, are constant, in this case, we can skip the t index, and consider the dataset ${\{t_i,y_i\}}_{i=1}^N$. In this scenario, as an alternative, one can use an autoregressive (AR) model,

$$\begin{array}{c}\hfill y_i=a_1{y}_{i-1}+a_2{y}_{i-2}+...+a_W{y}_{i-W}+{ϵ}_y,\end{array}$$

(8)

where ${ϵ}_y$ is a noise perturbation and the coefficients $a_i$, for $i=1,...,W$, are obtained by a least squares (LS) minimization and the length of the temporal window W (i.e., the order of the AR filter) is obtained using a spectral information criterion (SIC) [29,30]. Then, considering the estimated coefficients ${\widehat{a}}_i$ by LS, and $\widehat{W}$ by SIC, we have

$$\begin{array}{c}\hfill \widehat{f}\left(t_i\right)={\widehat{f}}_i={\widehat{a}}_1{y}_{i-1}+{\widehat{a}}_2{y}_{i-2}+...+{\widehat{a}}_{\widehat{W}}{y}_{i-\widehat{W}}.\end{array}$$

(9)

Just as an example, in order to apply M2 for the estimating the instant variance, we can set $v_i={(y_i-{\widehat{f}}_i)}^2$, and assume another AR model over $v_i$ with window length now denoted as H,

$$\begin{array}{c}\hfill v_i=b_1{v}_{i-1}+b_2{v}_{i-2}+...+b_H{v}_{i-H}+{ϵ}_v,\end{array}$$

(10)

where again the estimations ${\widehat{b}}_i$ can be obtained by LS, and $\widehat{H}$ by SIC. The resulting instant variance function is

$$\begin{array}{c}\hfill v\left(t_i\right)=v_i={\widehat{b}}_1{v}_{i-1}+{\widehat{b}}_2{v}_{i-2}+...+{\widehat{b}}_{\widehat{H}}{v}_{i-\widehat{H}}.\end{array}$$

(11)

Clearly, the application of M1 could be performed in a similar way.

Remark 5. 

The resulting estimators in this section are still a linear combination of (a portion) of the outputs. Therefore, we still have linear smoothers (or better linear filters, since only considering combinations of past samples). However, note that in Eqs. (9)-(11) the coefficients of the linear combinations are obtained by least squares (LS) method. Hence, we have a window of length W (and then H) but this differs from the use of a rectangular kernel function for two reasons: (a) the window only considers past samples; (b) the coefficients are not all equals (to the ratio 1 over the number of samples included in the window) but they are tuned by LS.

5. Numerical Experiments

5.1. Applications to Time Series

In this section, we consider 5 different numerical experiments with different (true) generating models. Here, we focus on time series ($x=t\in \mathbb{R}$) for two main reasons: (a) first of all, for the sake of simplicity (we can easily show a realization of the data in one plot), (b) and last but not the least, we have relevant benchmark competitors such as the GARCH models [31,32,33], that we can also test. Note that GARCH models have been specifically designed for handling volatility in time series.

More specifically, in this section, we test our methodologies and different GARCH models in four easy reproducible examples. These four examples differ for the data generating model. Indeed, the data are generated as the sum of a mean function, $f\left(t\right)$, and an error term, $s\left(t\right)\epsilon \left(t\right)$, i.e.,

$$\begin{array}{cc}\hfill y\left(t\right)& =f\left(t\right)+s\left(t\right)\epsilon \left(t\right),\epsilon \left(t\right)\sim \mathcal{N}(0,1).\hfill \end{array}$$

(12)

where we have used the fact that we are analyzing time series, i.e., $x=t\in \mathbb{R}$. The functions $f\left(t\right)$ and $s\left(t\right)$ are deterministic, and represent for the trend and the standard deviation of the observations $y\left(t\right)$. Whereas $\epsilon \left(t\right)\sim \mathcal{N}(0,1)$ represent a standard Gaussian random noise. Note also that $v\left(t\right)=s{\left(t\right)}^2$. In order to keep the notation of the rest of the paper, note that the dataset ${\{t_n,y_n\}}_{n=1}^N$ will be formed by the N pairs,

$$x_n=t_n,y_n=y\left(t_n\right),$$

where $y_n$ will be generated according to the model above. The four experiments consider:

• $f\left(t\right)=sin\left(0.5t\right)$, $s\left(t\right)=0.2$.
• $f\left(t\right)=t$, $s\left(t\right)=2\left|t\right|$.
• $f\left(t\right)=0.5t^2$, $s\left(t\right)=5sin\left(0.3t\right)$.
• $f\left(t\right)=0$ and $s\left(t\right)$ generated by a GARCH model.

In all cases, the generated data $y\left(t\right)$ and the underlying standard deviation $s\left(t\right)$ (or variance $v\left(t\right)=s{\left(t\right)}^2$) are known, hence we can compute errors and evaluate the performance of the used algorithms. For the application of the proposed methods we consider a nonlinear kernel function as

$$\begin{array}{c}\hfill h(t,t_n)=exp\left(-\frac{{(t-t_n)}^{\alpha }}{\lambda }\right),\end{array}$$

(13)

with $\alpha \in \{1,2,10\}$.1 The parameter $\lambda$ is tuned by LOO-CV. We test 7 different combinations, each one mixing a specific trend approximation method with a specific variance estimation technique:

• Based on kernel smoothers: we consider different regression methods based in Eqs. (1) and (13). Furthermore, for the variance estimation we consider both procedures M1 and M2 described in Section 3. More specifically:

  -

  KS-1-1: $\alpha =1$ and M1,

  -

  KS-1-2: $\alpha =1$ and M2,

  -

  KS-2-1: $\alpha =2$ and M1,

  -

  KS-2-2: $\alpha =2$ and M2,

  -

  KS-10-1: $\alpha =10$ and M1,

  -

  KS-10-2: $\alpha =10$ and M2.

• Based on AR models: we apply the method described in Section 4.2, that employs two auto-regressive models, one for the trend and one for the variance, jointly with procedure M2.

Table 1 summarizes all the seven considered specific implementations. Moreover, we test several GARCH(p,q) models: GARCH(1,1), GARCH(2,2), GARCH(5,5) and GARCH(10,10) (using a Matlab implementation of GARCH models).

Remark 6. 

The GARCH models work directly with the signal without trend, i.e., $|y_n-{\widehat{f}}_n|$. In order to have a fair comparison, we use a kernel smoother with $\alpha =2$, i.e., applying Eqs. (1)-(13) with $\alpha =2$ (as in the KS-2 methods).

We recall that GARCH models have been specific designed for making inference in heteroscedastic time series.

5.1.1. Experiment 1

As previously said, in this first example we consider

$$\begin{array}{c}\hfill y\left(t\right)=sin\left(0.5t\right)+0.2\epsilon \left(t\right),\epsilon \left(t\right)\sim \mathcal{N}(0,1),\end{array}$$

i.e., where $f\left(t\right)=sin\left(0.5t\right)$ and $s\left(t\right)=0.2$ ( in this case we have homoscedastic series). We generate $N={10}^4$ data with $t\in [-10,10]$. We apply all the algorithms, compute the mean absolute error (MAE), and repeat and average the results over ${10}^4$ independent runs. Figure 1 depicts one realization of the generated data with the model above.

Results 1.Figure 1b,c provide the results in this example. Figure 1 shows the MAE obtained in estimating the trend, and Figure 1 depicts the MAE in estimating the standard deviation. In this case, regarding the trad, we see that the AR-based has an MAE higher than the rest when estimating the trend, but the other MAEs are quite similar. Similar outcomes are also obtained for approximating with the standard deviation, although the smaller MAE are provided by the GARCH models. Regarding the estimating of the standard deviation, the M1 seems to provide better results than M2, in this example.

5.1.2. Experiment 2

in this second example, we consider

$$\begin{array}{c}\hfill y\left(t\right)=t+2\left|t\right|\epsilon \left(t\right),\epsilon \left(t\right)\sim \mathcal{N}(0,1),\end{array}$$

where $f\left(t\right)=t$ and $s\left(t\right)=2\left|t\right|$ (i.e., in this case, we are heteroscedastic scenario). As previously, we generate $N={10}^4$ data with $t\in [-10,10]$. We apply all the algorithms, compute the mean absolute error (MAE), and again we repeat and average the results over ${10}^4$ independent runs. Figure 2 depicts one realization of the generated data according to the model above.

Results 2.The results are presented in Figure 2b,c. Again the MAE for the AR-based method seems to be higher than those for the rest of the algorithms when estimating the trend and standard deviation of the time series. Regarding the estimation of the standard deviation the best results are provided by the kernel methods which use the M2 (the MAE is particularly small with$\alpha =10$). Thus, it is remarkable that, even in this heteroscedastic scenario, the proposed methods provide similar and, in same cases, better results than GARCH models.

5.1.3. Experiment 3

In this third example, we consider a standard deviation that varies periodically, and a second order polynomial as a trend function, i.e.,

$$\begin{array}{c}\hfill y\left(t\right)=0.5t^2+5sin\left(0.3t\right)\epsilon \left(t\right),\epsilon \left(t\right)\sim \mathcal{N}(0,1),\end{array}$$

where $f\left(t\right)=0.5t^2$ and $s\left(x\right)=5sin\left(0.3t\right)$. Again, we have heterostodastic scenario. As previously, we generate $N={10}^4$ data with $t\in [-10,10]$. We compute the mean absolute error (MAE), and again we repeat and average the results over ${10}^4$ independent runs. Figure 3 depicts one realization of the generated data according to the described model.

Results 3.Figure 3b,c depict the box-plots of the MAE in estimating the trend and the standard deviation. As with the rest of the examples, the AR-based gives the worst results when calculating the trend of the time series, meanwhile the rest of the algorithms have very similar error values (recall that the GARCH models employ the KS-2 method for the trend estimation). However, the AR-based method provides good performance in the approximation of the standard deviation. The AR-based technique uses M2 for the variance. Moreover, for the standard deviation estimation, the methods KS-1-1-1, KS-2-1, KS-10-1 gives the smaller errors. Hence, it seems that in this case M1 works better even than GARCH models.

5.1.4. Experiment 4

Let us consider

$$\begin{array}{cc}\hfill y\left(t\right)& =s\left(t\right)\epsilon \left(t\right),\epsilon \left(t\right)\sim \mathcal{N}(0,1),\hfill \\ \hfill y_t& =s_t{ϵ}_t=z_t,\hfill \end{array}$$

i.e., $f\left(t\right)=f_t=0$, and we have denoted $z_t=s_t{ϵ}_t$. The process for the standard deviation $s_t=\sqrt{v_t}$ will be a GARCH($P,Q$), i.e.,

$$\begin{array}{cc}\hfill v_t=& \kappa +{\gamma }_1{v}_{t-1}+\dots +{\gamma }_P{v}_{t-P}+\hfill \\ \hfill & +{\alpha }_1{z}_{t-1}^2+\dots +{\alpha }_Q{z}_{t-Q}^2.\hfill \end{array}$$

In this example, we generate the variance from a GARCH($2,3$)

$$\begin{array}{cc}\hfill v_t=& 2+0.2·v_{t-1}+0.1·v_{t-2}+0.2·z_{t-1}^2+\hfill \\ \hfill & +0.2·z_{t-2}^2+0.2·z_{t-3}^2\hfill \end{array}$$

(14)

Clearly, again, it is a heterostodastic scenario. we generate $N={10}^4$ data and apply the different techniques. We compute the mean absolute error (MAE), and we repeat and average the results over ${10}^4$ independent runs. One realization of the generated data is given in Figure 4. Note that the trend is linear but the volatility presents very fast changes since it generated by a GARCH(2,3) model (see blue line in Figure 4).

Results 4.The results for this example are presented in Figure 4b,c. The best performance in the estimation of$s\left(t\right)$are obtained by KS-1-2 and all the GARCH models (as expected in this case). The AR-based and KS-2-2 also provided close values of MAE. In this example, clearly M2 outperforms M1. We remark that, even in this scenario that is particularly favorable for the GARCH models, all the proposed methods provide competitive results. The highest (hence worst) MAEs values are given by the KS with$\alpha =10$which is a quite extreme tuning/choice of this parameter, specially when the volatility follows an autoregressive rule. Indeed, a Laplacian kernel with$\alpha =1$is more adequate in this framework. For further considerations see [15].

5.1.5. Final Comments about the Applications to Time Series

All the proposed methods based on kernel smoothers (KS) and M1, M2 always provide competitive results (with closer or smaller MAE values) with respect to the benchemark algorithms such as GARCH methods, which have been specifically designed for the estimation of the standard deviation in time series.

These numerical experiments show that the proposed M1 and M2 can even outperform GARCH models the estimation of the standard deviation in time series. The GARCH models seems to provide more robust results when the signal has been generated truly by a GARCH model in Section 5.1.4 and, surprisingly, in the homoscedasticscenario in Section 5.1.1. The methods based on kernel smoothers virtually outperform the AR-based due to the fact that they incorporate the “future samples” in their estimators. A rectangular window, including both past and future samples, should improve the performance of the AR-based method.

5.2. Application in Higher input Dimension: Real Data on Soundscape Emotions

GARCH models have been specifically designed for estimating the volatility in time series, hence $x=t\in \mathbb{R}$. However, the methods proposed in this work can be applied to problems with multidimensional inputs $\mathbf{x}\in {\mathbb{R}}^D$ with $D\ge 1$.

More specifically, here we focus on analyzing a soundscape emotion database. In the last decade, soundscapes have become one of the most active topics in acoustics. In fact, the number of related research projects and scientific articles grows exponentially [34,35,36]. In urban planning and environmental acoustics, the general scheme consists of (a) soundscapes recording, (b) computing of acoustic and psychoacoustic indicators of the signals, (c) including other context indicators (e.g. visual information [37]), and (d) ranking of soundscapes audio signals employing emotional descriptors. Finally, the model can be developed [24].

Here, we consider the emo-soundscapes database (EMO) [27,28], where $N=1200$ and $D=122$, i.e., $\mathbf{x}\in {\mathbb{R}}^{122}$, which is the largest available soundscape database with annotations of emotion labels, and also the most recent, up to now [27,38]. Among the 122 audio features that form the inputs $\mathbf{x}$, we can distingush three main groups of features of the audio signals: psychoacoustic features, time-domain features and frequency-domain features. The output y that we consider is called arousal. For further information, the complete database can be found at https://metacreation.net/emo-soundscapes/. We employ again a kernel function of type

$$h(\mathbf{x},\mathbf{z})=h(\mathbf{x},\mathbf{z}|\lambda )=exp\left(-\frac{\left|\right|\mathbf{x}-{\mathbf{z}\left|\right|}^{\alpha }}{\lambda }\right),$$

with $\alpha =2$. Therefore, considering the EMO database ${\{{\mathbf{x}}_i,y_i\}}_{i=1}^{1200}$ We apply M1 and M2 and CV-LOO for tuning $\lambda$. The MAE in estimating the trend is $0.0388$ obtained with an optimal parameter ${\lambda }^{*}=0.12$. Moreover, in both cases, we find a sensible increase of the volatility of the output in the last half of the samples, for i from 621 to 1213. This could be an interesting result, that should be discussed with experts on the field and in the EMO database.

6. Conclusions

In this work, we have proposed two methods for estimating the instant/local variance in regression problems based on kernel smoothers. From another point of view, this work can be also seen as a way to perform quantile regression (at least, a second order quantile regression) using kernel smoothers.

With respect to other procedures given in the literature, the proposed methods have a much wider applicability than other techniques, since they can be employed only for time series, as the well-known GARCH models (where the time index is a scalar number, $x=t\in \mathbb{R}$). Indeed, the proposed schemes can be employed with multidimensional input variables ${\mathbf{x}}_n\in {\mathbb{R}}^D$. Moreover, the proposed techniques can be also applied to time series ($D=1$) and also can obtain a prediction between two consecutive time instants (i.e., with an higher resolution that the received data). More generally, estimation of the local variance can given at any the training data ${\mathbf{x}}_i$ but also in any generic test input $\mathbf{x}$.

Furthermore, analyzing time series, the numerical simulations have shown that the proposed schemes provides competitive results even with respect to GARCH models (which are specifically design for analyzing time series). An application to real data and multidimensional inputs (with dimension 122) have been also provided. As a future research line, we plan to design kernel smoother functions with local parameters $\lambda$, or more generally with $\lambda =\lambda \left(\mathbf{x}\right)$, to better handle the heteroscedasticity. Moreover, the application to multioutput scenarios can be designed and considered.