# Short-term prediction of Time Series based on bounding techniques

Pedro Cadahía<EMAIL_ADDRESS>José M. Bravo  Escuela Técnica
Superior de Ingeniería, Universidad de Huelva,
Carretera Huelva - Palos de la Frontera s/n. 21819. La Rábida - Palos de la
Frontera. Huelva. Spain

###### Abstract

In this paper it is reconsidered the prediction problem in time series
framework by using a new non-parametric approach. Through this
reconsideration, the prediction is obtained by a weighted sum of past observed
data. These weights are obtained by solving a constrained linear optimization
problem that minimizes an outer bound of the prediction error. The innovation
is to consider both deterministic and stochastic assumptions in order to
obtain the upper bound of the prediction error, a tuning parameter is used to
balance these deterministic-stochastic assumptions in order to improve the
predictor performance. A benchmark is included to illustrate that the proposed
predictor can obtain suitable results in a prediction scheme, and can be an
interesting alternative method to the classical non-parametric methods.
Besides, it is shown how this model can outperform the preexisting ones in a
short term forecast.

###### keywords:

Nonparametric methods , Nonlinear models , Optimization , Time series ,
Univariate predicting method

left=1.5in,right=1.5in,top=1.5in,bottom=1.5in

##  1 Introduction

The purpose of this paper is to provide a new model for time series set up on
the observed past values of the time series, by means of a non-parametric
approach. It is well-known fact that in parametric time series analysis the
relationship between observed past values of the time series and the
prediction is defined by specifying a functional form and a fixed finite
number of parameters. Widely studied parametric options are auto-regressive
(AR) models, moving average (MA) models, and different combinations as ARMA or
ARIMA models [1, 2]. In nonlinear time series, some common parametric
structures has been studied, the threshold auto-regressive (TAR) models [3],
the exponential auto-regressive (EXPAR) model and smooth-transition auto-
regressive (STAR) models are some examples [4, 5]. The performance of the
parametric predictor is a consequence of the a priori function form chosen.

By contrast, in non-parametric approaches a more flexible class of functions
is considered. Non-parametric methods avoid the choosing of a specific
functional form. Collected data provides the information to obtain a new
prediction. The price to pay is the ’curse of dimensionality’, that is, a
possible poor performance in high dimensions prediction problems. Local
conditional mean or median method provides a prediction using the mean or the
median of a neighborhood of the interest point [6]. The Nadaraya-Watson
estimator averages past observations by a kernel function to obtain a
prediction [7, 8]. Local linear o polynomial functions of past observations
can be used to approximate a nonlinear relationships [9, 10]. Semi-parametric
models as nonlinear additive auto-regressive (NAAR) models or functional
coefficient auto-regressive (FAR) models have been proposed too [11, 12]. Many
researchers have written an extensive review of non-parametric methods applied
to time series prediction [13, 14, 15, 16].

In this paper a new non-parametric prediction method is proposed. The
prediction is obtained by a weighted sum of past observations. An upper bound
of the prediction error is computed under some deterministic and stochastic
assumptions. A constrained optimization problem is formulated to minimize the
upper bound of the prediction error and to obtain the set of optimal weights
used to compute the prediction. The optimization problem includes a parameter
to balance the deterministic-stochastic assumptions. This is the main novelty
of the proposed method. This parameter can be tuned with training data and a
cross-validation scheme to improve the predictor performance [17]. The
proposed predictor provides a general framework that encompasses some relevant
non-parametrics predictors as the Nadaraya-Watson predictor [7, 8] or
predictors based on local linear regression [10], these models have been widle
used in the literature [18, 19].

The paper is organized as follows. In Section 2, the problem formulation is
addressed. The deterministic and stochastic assumptions are presented in
Section 3. The new predictor is proposed in Section 4. Benchmark results are
illustrated in Section 5. Finally, Section 6 reports some conclusions.

##  2 Formulation

It is considered a discrete111It is assumed a discrete version of data. time
series process $\\{Z_{t}\\}$ with $t\in\\{0,\pm 1,\pm 2,\ldots\\}$. At time
instant $k$ it is assumed that past data $\\{Z_{t}\\}$ with
$t\in\\{k,k-1,k-2,...\\}$ has been observed and there is interest in providing
a forecast for predicting $Z_{k+1}$. Once the detrend is applied to the time
series,222It should be noted that in coherence with the prediction system and
in order to estimate $\mu_{k+1}$, only the past observations can be used,
independently of the detrending method used. the time series is now the series
$\\{y_{t}\\}$ with $t\in\\{k,k-1,...\\}$, where $Z_{t}=y_{t}+\mu_{t}$, being
$\mu_{t}$ the trend component and $y_{k+1}$ the detrended future time series
value.

It is also denoted by $\\{z_{j}\\}$ with ${j=0,1,...,k}$ the set of the
vectors consisting of the observed past values of the time series, that is
$z_{j}=[y_{j},y_{j-1},\ldots,y_{j-p-1}]^{T}$ . Henceforth this $p$-dimensional
vector set will be called embedding vector. This set of data is used to
forecast future values for the time series. It is a must to clarify this point
in order to precise the sense of the parametric and non-parametric models used
in this article. A parametric approach is characterized by the use of the
training set for estimating the parameters of the model and once this
inference is done the data set is not used again. The non-parametric approach
considered in this work, it is a local approach in which each forecast is
obtained by using all the available data set but selecting a neighborhood of
the interest point. In this sense, it is assumed that the time series can be
generated by an unknown local linear model.

###### Assumption 1

Considering the forecast of $y$ modeled as:

$y_{k+1}=r(z_{k})^{T}\Phi_{k}+e_{k}$ (1)

where it is assumed that the existence of an unknown vector of parameters
$\Phi_{k}\in\mathcal{R}^{n}$, a known function $r(\cdot)$ valuated at the
embedding set and an unknown error term $e_{k}$.333This modeling is flexible
enough to admit alternative assumptions about the error term. As discussed
later, the model is presented by using both deterministic and stochastic
bounds for the error term $e_{k}$.

In order to complete the presentation of the model it should be discussed in
more detail the so called regressor generator function $r(\cdot)$. This
function allows transform the original values into vectors of dimension
$n_{r}$ by means of the vectors belonging to the embedding set. A formal
definition of this regressor generator function is as follows.

######  Definition 1 (Regressor generator function)

The function $r(\cdot):\mathcal{R}^{p}\rightarrow\mathcal{R}^{n_{r}}$
specifies the regressor vector components. This function admits any kind of
auto-regressive representation, nonlinear expression of past components and
different functional forms for decomposing the different components of the
time series.444For instance suppose a set $z_{k}=[y_{k},y_{k-1},y_{k-2}]$.
Then $r(z_{k})$ could be the function $r(z_{k})=z_{k}$ that is, an auto-
regressive model. There exist also alternative configurations such as a
nonlinear auto-regressive model
$r(z_{k})=[y_{k}^{2},y_{k-1},y_{k-2},y_{k}\cdot y_{k-2}]$ or any possible
combination.

######  Definition 2 (Linear Prediction)

For an instant $k$, a forecast of $y_{k+1}\in\mathbb{R}$ can be derived
through a linear combination of past data, that is:

$\begin{array}[]{cll}\hat{y}_{k+1}(\Psi)&=&b_{Y}^{T}\Psi\\\
&=&\displaystyle\sum\limits_{j=1}^{v}\Psi_{j}y_{j}\end{array}$ (2)

where $1\leq v\leq k$ , $\Psi\in\mathbb{R}^{v}$ is a weight vector and
$b_{Y}=[y_{1},\ldots,y_{v}]^{T}$.

When $v=k$, all data is used to forecast $y_{k+1}$. Then, the forecast error
can be explained as the difference between $y_{k+1}$ and the linear prediction
$\hat{y}_{k+1}(\Psi)$.

######  Definition 3 (Prediction error)

It is defined the prediction error $\hat{e}_{k}(\Psi)$, being $k$ the time
instant:

$\hat{e}_{k}(\Psi)=y_{k+1}-\hat{y}_{k+1}(\Psi).$ (3)

Thus, the crux of the matter is how to get not only the weight vector $\Psi$
but also an outer limit of the prediction error. This outer limit is estimated
by using the assumed relationship between $z_{j-1}$ and $y_{j}$, with
$j=1,2,\ldots,k$ in expression (1). Then, a set of past components $z_{j}$
with $j=0,1,...,k$ should be available. Section 3 formulates these key ideas.

##  3 Assumptions

In this section the assumptions are based on some local affine approximations.
In order to construct the proposed predictor, the definition of approximation
error is used. This is, the result of using the vectors $r(z_{j-1})$ and
$\Phi_{k}$ to infer $y_{j}$.

555The reader should note that the point is to relate the k-th prediction
error $e_{k}$ and the prediction errors generated by using the k-th vector of
unknown parameters $\Phi_{k}$ with the i-th regressors $r(z_{i})$, with
$i=0,\ldots,k-1$.

######  Definition 4 (Approximation error)

For a vector $\Phi_{k}$, the approximation error $e_{j-1}$ with the pair
$(z_{j-1},y_{j})$ being $j=1,2,...,k$ can be defined as:

$e_{j-1}=e_{j-1}(\Phi_{k})=y_{j}-r(z_{j-1})^{T}\Phi_{k}.$ (4)

From now on the dependency of $e_{j-1}(\Phi_{k})$ with $\Phi_{k}$ is omitted.
It should be noted that the value of $\Phi_{k}$ is unknown. The prediction
error $\hat{e}_{k}(\Psi)$ may be biased by the selected vector $\Psi$. The
theorem 1 suggests an approach to define the prediction error
$\hat{e}_{k}(\Psi)$ as a function of the vector $\Psi$ and the aforementioned
approximation errors $e_{j}$ .

Theorem 1 proposes an expression to characterize the prediction error
$\hat{e}_{k}(\Psi)$ as a function of vector $\Psi$ and approximation errors
$e_{j}$ previously defined.

###### Theorem 1

For either vector $\Psi\in\mathbb{R}^{v}$ so that

$\displaystyle\sum\limits_{j=1}^{v}\Psi_{j}r(z_{j-1})=r(z_{k}),$ (5)

so the prediction error $\hat{e}_{k}(\Psi)=y_{k+1}-\hat{y}_{k+1}(\Psi)$ is set
as a linear combination of the approximation errors $e_{j}$, this is

$\hat{e}_{k}(\Psi)=-\displaystyle\sum\limits_{j=1}^{v}\Psi_{j}e_{j-1}+e_{k}.$

Remark that $\Psi_{i}$ refers to the j-$th$ item of vector $\Psi$. A proof of
the theorem can be found in the Appendix section 7.1. Matricially, expression
(5) is equivalent to $\Psi\in\\{\Psi\;:\;A^{T}\Psi=r(z_{k})\\}$ where matrix
$A$ is:

$A^{T}=\left[\begin{array}[]{ccccccc}r(z_{0})&r(z_{1})&...&r(z_{v-1})\\\
\end{array}\right].$ (6)

It is necessary to know the vector $\Phi_{k}$ to get an error value $e_{j-1}$.
Alternatively, other properties of $e_{j-1}$ can also be assumed. Both
deterministic and stochastic options are available in the literature. In a
deterministic view, an upper bound of $|e_{j-1}|$ is considered. This idea is
discussed in the section 3.1.

###  3.1 Deterministic error

In methods with bounded-error [20], a parametric model and an unknown but
bounded-error are regarded. An upper limit of this error is expected to
estimate a set of consistent parameters. Similar assumptions are presumed in
this work in order to develop a predictor with deterministic assumptions.

###### Assumption 2

Constants $\sigma,L\geq 0$ are set such that approximation errors $e_{j-1}$
and $e_{k}$ are delimited by expressions

$|e_{j-1}|\leq\sigma+L||z_{j-1}-z_{k}||$ (7)

with $j=1,...,k$ and

$|e_{k}|\leq\sigma$ (8)

being $||\cdot||$ a norm.

The error term is bounded by $|e_{k}|\leq\sigma$. The assumption 2 has been
broadly used in the bounded-error system identification’s context [20]. Remark
that $\sigma$ is the tunning parameter that adds the minimum level of noise
considered and $L$ the tunning parameter of uncertainty due to the local
affine approximation.

###### Remark 1

Historical data can be used to estimate an approximate value of $\sigma$ and
$L$ when no prior knowledge of these constants is available. In [21] a method
based on bounded-error and non-counterfeit data is provided.

###### Lemma 1

Considering Assumptions 1 and 2, for any $\Psi$ such that
$A^{T}\Psi=r(z_{k})$, prediction error
$\hat{e}_{k}(\Psi)=y_{k+1}-\hat{y}_{k+1}(\Psi)$ is bounded by:

$|\hat{e}_{k}(\Psi)|\leq\displaystyle\sum\limits_{j=1}^{v}|\Psi_{j}|(\sigma+L||z_{j-1}-z_{k}||)+\sigma.$
(9)

Proof.Through a straightforward application of Theorem 1 and bound
$|e_{i}|\leq\sigma+L||z_{j}-z_{k}||$ is obtained the expression (9). _QED_

At this point the possibility of considering how to obtain the vector $\Psi$
is established. A wise option is to use the vector that minimizes an upper
bound of $|\hat{e}_{k}(\Psi)|$ using the expression (9).

######  Definition 5 (Deterministic predictor)

The deterministic prediction $\hat{y}_{k+1}(\Psi^{D})$ is defined by

$\hat{y}_{k+1}(\Psi^{D})=\displaystyle\sum\limits_{j=1}^{v}\Psi^{D}_{j}y_{j},$

where vector $\Psi^{D}$ adresses the problem of constrained linear
optimization as follows

$\begin{array}[]{ccc}\Psi^{D}=&arg\min\limits_{\Psi}&||W_{k}\Psi||_{1}\\\
&s.t.&A^{T}\Psi=r(z_{k})\\\ \end{array}$ (10)

where $W_{k}$ is a diagonal matrix with central items
$w^{k}_{j,j}=\sigma+L||z_{j-1}-z_{k}||$ with $j=1,...,v$. Then, an upper bound
of the absolute value of the prediction error is minimized by the vector
$\Psi^{D}$.

It is important to note that the notation $\Psi^{D}$ refers to the
deterministic nature of the estimate. Expression (10) use $L_{1}$-norm to
obtain the vector solution $\Psi^{D}$. In this case, $\Psi^{D}$ is sparse,
that is, most of number of components $\Psi^{D}_{i}$ of vector $\Psi^{D}$ are
zero. As $\Psi^{D}$ is a sparse matrix and considering Definition 5 then it is
deduced that $\hat{y}_{k+1}(\Psi^{D})$ use a relatively short number of
measurements $y_{i}$.

###  3.2 Stochastic error

The stochastic view consider the approximation error $e_{j}$ as a random
variable. So there are some assumptions about the mean and the variance of
$e_{j}$. Specifically there are assumptions in the dimension of variance of
$e_{j}$.

###### Assumption 3

The independent variables, approximation error $e_{j-1}$ and error term
$e_{k}$ have zero mean and variances bounded by
$var(e_{j-1})\leq\sigma+L||z_{j-1}-z_{k}||$ and $var(e_{k})\leq\sigma$
accordingly. Positive values of constants $\sigma$ and $L$ is taken as prior
knowledge.

As indicated in Remark 1, if not available previous knowledge of the constants
$\sigma$ and $L$, historical data may be used to obtain an estimation. The
variance of error $e_{j-1}$ consists of a minimum value defined by $\sigma$
and a term depending of the local approximation, i.e. $||z_{j-1}-z_{k}||$. it
is possible to extend that as $e_{j-1}$ and $e_{k}$ are random variables then
$\hat{e}_{k}(\Psi)$ is also random and therefore other properties can be
derived.

###### Assumption 4

Taking into account the previous Assumptions 1 and 3, for any $\Psi$ such that
$A^{T}\Psi=r(z_{k})$, prediction error
$\hat{e}_{k}(\Psi)=y_{k+1}-\hat{y}_{k+1}(\Psi)$ is a random variable with zero
mean and variance, it is defined by:

$\begin{array}[]{cll}var(\hat{e}_{k}(\Psi))&=&\displaystyle\sum\limits_{j=\underline{1}}^{v}\Psi_{j}^{2}var(e_{j-1})+\sigma\\\
&\leq&\displaystyle\sum\limits_{j=1}^{k}\Psi_{j}^{2}(\sigma+L||z_{j-1}-z_{k}||)+\sigma.\\\
\end{array}$ (11)

At this point,it is possible to formulate a predictor that minimize the outer
bound of the variance prediction error.

######  Definition 6 (Stochastic prediction)

The stochastic prediction $\hat{y}_{k+1}(\Psi^{S})$ is defined by:

$\hat{y}_{k+1}(\Psi^{S})=\displaystyle\sum\limits_{j=1}^{v}\Psi^{S}_{j}y_{j},$

being $\Psi^{S}$ a vector that solves a constrained linear optimization
problem as follows:

$\begin{array}[]{ccc}\Psi^{S}=&arg\min\limits_{\Psi}&\Psi^{T}W_{k}\Psi\\\
&s.t.&A^{T}\Psi=r(z_{k}).\\\ \end{array}$ (12)

An explicit notation of this optimization problem is:

$\Psi^{S}=W_{k}^{-1}A(A^{T}W_{k}^{-1}A)^{-1}r(z_{k}).$ (13)

In the same way, $\Psi^{S}$ highlights the stochastic assumptions considered
to get the estimate. The following equality is satisfied

$\hat{y}_{k+1}(\Psi^{S})=b_{Y}^{T}\Psi^{S}=r(z_{k})^{T}\Phi^{*},$

where $\Phi^{*}=(A^{T}W_{k}^{-1}A)^{-1}A^{T}W_{k}^{-1}b_{Y}$ is the argument
of which minimizes the a quadratic prediction-error, with the following cost
function:

$\begin{array}[]{rl}J(\Phi)=&(b_{Y}-A\Phi)^{T}W_{k}^{-1}(b_{Y}-A\Phi)\\\
=&\displaystyle\sum\limits_{j=1}^{k}\frac{(y_{j}-r(z_{j-1})^{T}\Phi)^{2}}{(\sigma+L||z_{j-1}-z_{k}||)}.\end{array}$
(14)

In this way, the stochastic prediction is equivalent to solve a weighted
least-squares problem where the weights are set by the items of the diagonal
of $W_{k}$ squared. Commonly, $\Psi^{S}$ is not a sparse vector, this is that
most items are non zero numbers. So, in order to get the prediction $y_{k+1}$
a great number of $y_{j}$ would be used.

The goal of this paper is to bring a predictor that combines the two
predictions based on the different assumptions obtained from $\Psi^{D}$ and
$\Psi^{S}$ respectively. Section 4 introduces the key points of this paper.

##  4 Proposed predictor

This work proposes to obtain an estimation of the output $y_{k+1}$ by a linear
combination of past data $y_{j}$, with $j=1,2,...,v$ where $v\leq k$ ([22]).
Next, a formal definition of the proposed predictor is provided. This
definition use a constant $\gamma\geq 0$ to balance the deterministic or
stochastic nature of the prediction.

###### Definition 7

Given a constant $\gamma\geq 0$, the predictor $\hat{y}_{k+1}(\Psi^{*})$ is
defined by
$\hat{y}_{k+1}(\Psi^{*})=\displaystyle\sum\limits_{j=1}^{k}\Psi^{*}_{j}y_{j}$
where $\Psi^{*}$ is the optimal solution of:

$\begin{array}[]{ccc}\Psi^{*}(\gamma)=&arg\min\limits_{\Psi}&||W_{k}\Psi||_{1}\\\
&s.t.&A^{T}\Psi=r(z_{k})\\\ &&||\Psi-\Psi^{S}||_{1}\leq\gamma\end{array}$ (15)

and vector $\Psi^{S}$ is defined in (13).

Some qualitative properties of the proposed predictor can be clarified. Note
that, expression (15) is a constrained linear convex optimization problem and
can be solved in an efficient way [23]. Assuming that (15) has a bounded
solution, there is a constant $\bar{\gamma}$ such that if
$\gamma\geq\bar{\gamma}$ then equality $\Psi^{*}=\Psi^{D}$ is obtained. Term
$||\Psi-\Psi^{S}||_{1}$ of expression (15) takes into account the stochastic
Assumption explained in Section 3 to obtain the optimal solutions $\Psi^{*}$.
If $\gamma=0$ then $\Psi^{*}=\Psi^{S}$. So, constant $\gamma$ can be seen as a
tuning parameter to balance the deterministic or stochastic nature of the
considered approximation error.

###### Remark 2

It is important to remark that the proposed predictor encompasses some
relevant nonparametrics predictors. If $\gamma=0$ and $r(z_{k})=1$ the
proposed predictor is equivalent to the Nadaraya-Watson predictor [7, 8]. On
the other hand if $\gamma=0$ and $r(z_{k})=[z_{k}^{T}\;1]$ a predictor based
on Local Linear Regression is obtained. Besides, if $\gamma=0$ and if $L=0$ a
parametric auto-regressive linear regression is performed.

###### Remark 3

It is important to remark that following similar reasoning it is possible to
obtain different forecasting horizons. This is represented by the expression
$y_{k+h}(\Psi^{*})$ being $h\geq 1$ the number of steps ahead.

##  5 Results

In this section results are shown, the predictor was performed in four time
series: the Monthly airline passenger numbers, the Canadian lynx data, the
Monthly critical radio frequencies in Washington, D.C. and the Monthly
pneumonia and influenza deaths time series are used in order to demonstrate
the appropriateness and effectiveness of the proposed predictor. These time
series come from different areas and have different statistical properties, so
is a suitable benchmark to test time series predictors.

Subsection 5.1 explain the characteristics of the study performed, hyper-
parameterization, kernels, error measures and more details are exposed below.

###  5.1 Considerations

  * 1.

To simplify the study, the proposed predictor (denoted $CP$) is considered
with values $\sigma=0$ and $L=1$ in all cases. Note that $\sigma$ and $L$
could be considered hyper-parameters in order to improve the results obtained
by the proposed predictor in this study.

  * 2.

The proposed predictor ($CP$) is compared to three Nadaraya-Watson predictors
(denoted as $NW1$, $NW2$ and $NW3$) using Epanechnikov, Gaussian and Tricube
kernel functions respectively and three local linear regression models
(denoted as $LL_{1}$, $LL_{2}$ and $LL_{3}$) using Epanechnikov, Gaussian and
Tricube kernel functions respectively to define the local weights. Table 1
shows the expression of weights $w_{i,i}$ with $i=1,...,N$ for the
aforementioned kernel functions. A bandwidth $\gamma$ is considered in the
non-parametric predictors. Also, an auto-regressive linear regression $(AR)$
is indirectly included in the benchmark, as the proposed predictor also
includes this model according to the hyperparameter combinations as explained
in remark 2 of section 7.

Table 1: Kernel functions Epanechnikov | Gaussian | Tricube
---|---|---
$w_{i,i}=\left\\{\begin{array}[]{cc}1-v_{i}^{2}&if\;|v_{i}|\leq 1\\\ 0&if\;|v_{i}|>1\end{array}\right.$ | $w_{i,i}=e^{-\frac{1}{2}v_{i}^{2}}$ | $w_{i,i}=\left\\{\begin{array}[]{cc}(1-|v_{i}|^{3})^{3}&if\;|v_{i}|\leq 1\\\ 0&if\;|v_{i}|>1\end{array}\right.$
$v_{i}=\frac{||z_{i}-z_{k}||}{\gamma}$
  * 3.

Two different forecast consistency measures are used in order to compare the
predictor performances with the aforementioned models: Mean Absolute Error
(MAPE) and Symmetric Mean Absolute Error (SMAPE) that have been studied by
several authors [24]. Mean Absolute Error is defined by

$\begin{array}[]{ccc}MAPE=\frac{100}{n}\displaystyle\sum_{t=1}^{n}\left|\frac{y_{t}-\hat{y}_{t}}{y_{t}}\right|,\end{array}$
(16)

where $\hat{y}_{t}$ and ${y_{t}}$ are the predicted and observed data,
respectively, and $n$ is the number of data. The second criterion is the
Symmetric Mean Absolute Percentage error (SMAPE), which is

$\begin{array}[]{ccc}SMAPE=\frac{100}{n}\displaystyle\sum_{t=1}^{n}\frac{\left|\hat{y}_{t}-y_{t}\right|}{(|y_{t}|+|\hat{y}_{t}|)/2},\end{array}$
(17)

where $\hat{y}_{t}$ and ${y_{t}}$ are the predicted and observed data,
respectively, and $n$ is the number of data.

  * 4.

In order to train the predictors, the time series has been splitted into a
training and test set with a percentage between 70-90 for a training set and a
percentage between 10-30 of data for a test set. The aforementioned error
measures described in formulas (16) and (17) are selected to benchmark the
predictors, not only because of its interpretability but also because are
scale-independent [25]. The training set, a leave-one-out cross validation
approach, and a grid-search in the hyper-parameter-$\gamma$ space are used to
find a suitable value of hyper-parameter $\gamma$. This value is used in the
test set to evaluate the prediction methods.

  * 5.

As explained in Section 4, the estimation can be used by a linear combination
of the past data, in this way the different models are evaluated by using only
the training set to infer the prediction this is when past data $y_{j}$ is
used with $v<k$. Besides, three different forecast horizons are computed by
predictor (1 step-ahead, 2 step-ahead and 3 step-ahead).

In order to test the proposed predictor, two different type of data are used
to perform the benchmark, the first in subsection 5.2 where some famous and
more academic time series are used to compare the models, the second in
subsection 5.3 a real life time series is used to test the predictor.

###  5.2 Academic Time series

This subsection compare 6 non-parametric models against the proposed model in
four different time series, performing forecasts in three different predicting
horizon lengths. Each time series provides the results in bar-plots by $MAPE$
and $SMAPE$ errors in the test set for the proposed forecasting horizons. A
final sub-subsection averages all time series results in order to extract
general conclusions for the different benchmarked models.

####  5.2.1 Airline passengers dataset

The classic Box and Jenkins airline data contains monthly totals of
international airline passengers from $1949$ to $1960$ [1].

This time series plotted in Figure 1(a) has $144$ observations, the first
$101$ observations were used as training set and the last $43$ as test set. It
has also been very analyzed in the time series literature.

As explained in Section 2, it is considered a detrended time series. In this
sense a log with base 10 and a linear detrend function are applied to
transform the data, it is shown the data set transformed in Figure 1(b).

(a) Original Time series

(b) Transformed Time series

Figure 1: Monthly totals of international airline passengers (1949 $-$ 1960).

The auto-correlation plot of the transformed data set shown in figure 2 at the
appendix represents shows a high correlation between observations of this time
series that are separated by $k=12$ time units, in this sense, the predictor
can be represented as $r(z_{k})=z_{k}=[y_{k-1}\;y_{k-2}\;...\;y_{k-12}]^{T}$.

Figure 2: Auto-correlation function of Airline passengers transformed time
series.

Figure 3 shows the forecasts of the proposed predictor by forecasting horizons
with the hyper-parameters selected in both error measures.

Figure 3: International airline passengers predictions by forecasting horizon
in the test set.

The hyper-parameter $\gamma$ is selected in the training set where the error
is minimum, the value of $\gamma$ is inferred to perform forecasts in the test
set. Besides, depending on the error measure selected in the training set the
results may vary, in this case, same optimal hyper-parameters are found in
both error measures, these results are in table 2.

Table 2: Airline passengers time series optimal gamma. Ahead |  $\gamma$_mape |  $\gamma$_smape
---|---|---
1.00 | 0.12 | 0.12
2.00 | 0.14 | 0.14
3.00 | 0.00 | 0.00

Results of this time series are shown on Table 7 in the appendix. To sum up
the aforementioned table in a graphical way, figure 4 show the error measures
by predictor and prediction horizon.

(a) SMAPE

(b) MAPE

Figure 4: Mean of errors by forecasting horizon in airline passengers time
series in test set.

Results show by forecasting horizon that the proposed predictor outperforms in
both error measures in the proposed forecasting horizons, with the exception
of $SMAPE$ error criteria in the one and two step-ahead prediction, where the
prediction is close to $LL$ results. Besides, results show that there is not a
significant variation of results with the selection of the different kernels.

####  5.2.2 Canadian Lynx

The Annual numbers of lynx trappings in Canada, contains the number of lynx
trapped per year in the Mackenzie River district of Northern Canada from
$1821$ to $1934$ [26].

It has been extensively analyzed in the time series literature with a focus on
the nonlinear modeling. The lynx series plotted in Figure 5(a) shows a
periodicity of approximately 10 years. The lynx series was studied by many
researchers found the best-fitted model is AR(12) model [27]. In this way, the
predictor is based on a auto-regressive model of order $p=12$, this is
$r(z_{k})=z_{k}=[y_{k-1}\;y_{k-2}\;...\;y_{k-12}]^{T}$.

The lynx series plotted in Figure 5(a) has $114$ observations, the first 80
observations of this data set were used as training set and the last 34 as
test set.

A log with base 10 was applied to the series in order to make a symmetrical
data set, the plot of the data set is at Figure 5(b):

(a) Original Time series

(b) Transformed Time series

Figure 5: Annual number of lynx trappings in Canada from 1821 to 1934.

Figure 5 shows the forecasts of the proposed predictor by forecasting horizons
with the hyper-parameters selected in both error measures, in this case,
different optimal hyper-parameters are found in the error measure selection on
the training set.

(a) MAPE selection criterion

(b) SMAPE selection criterion

Figure 6: Canadyan lynx time series predictions by forecasting horizon in the
test set.

Table 3 shows results that comes from different optimal gamma selections by
different error criteria and forecasting horizons.

Table 3: Canadian Lynx time series optimal gamma. Ahead |  $\gamma$_mape |  $\gamma$_smape
---|---|---
1.00 | 0.01 | 0.06
2.00 | 0.00 | 0.02
3.00 | 0.02 | 0.04

Corresponding to this predictions, the error measures are shown on Table
LABEL:table:lynx_table in the appendix. To sum up this table in a graphical
way, the figure 7 plot the error measures by predictor and prediction horizon.

(a) SMAPE

(b) MAPE

Figure 7: Mean of test set errors by forecasting horizon in Canadian lynx time
series.

Results in figure 7 show that the proposed predictor get similar results to
Local Linear regression in both error measures in the two first proposed
forecast horizons, in this case a well-selected kernel on Local Linear
regression could be a good option to take into account in order to get closer
results in a short prediction term to the proposed predictor $CP$, in a three
step-ahead prediction horizon the results mark a tie between $LL$ and $CP$ in
the $SMAPE$, on the contrary by selecting $MAPE$ criterion the $LL$
outperforms.

####  5.2.3 Monthly critical radio frequencies

Monthly critical radio frequencies in Washington, D.C., contains the highest
radio frequency that can be used for broadcasting from May $1934$ to April
$1954$ [28].

This time series plotted in Figure 8 has $240$ observations, the first $216$
observations were used as training set and the last $24$ as a test set.

Figure 8: Monthly critical radio frequencies (1934$-$1954).

According to auto-correlation plot attached in the Figure 9 at the appendix,
the established model is based on an auto-regressive model of order twelve,
which has also been used by many researchers [29, 27]. The auto-regressive
model has the shape like
$r(z_{k})=z_{k}=[y_{k-1}\;y_{k-2}\;...\;y_{k-12}]^{T}$.

Figure 9: Auto-correlation function of Monthly critical radio frequencies time
series.

Prediction are shown in figures 10 for the proposed predictor by forecasting
horizons with the hyper-parameters selected in both error measures, in this
case, different optimal hyper-parameter are found in the different error
measure criteria selected on the training set.

(a) MAPE selection criterion

(b) SMAPE selection criterion

Figure 10: Monthly critical radio frequencies time series prediction by
forecasting horizon in the test set.

Prediction error results of this time series are recorded on Table 9 in the
appendix, table 4 shows the optimal selected gamma.

Table 4: Monthly critical radio frequencies time series optimal gamma. Ahead |  $\gamma$_mape |  $\gamma$_smape
---|---|---
1.00 | 0.00 | 0.00
2.00 | 0.03 | 0.00
3.00 | 0.00 | 0.00

To sum up the table 9 in a graphical way, figure 11 show the error measures by
predictor and forecasting horizon.

(a) SMAPE

(b) MAPE

Figure 11: Mean of test set errors by forecasting horizon in monthly critical
radio frequencies time series.

The proposed predictor get similar results to Local Linear regression in SMAPE
measure in the proposed forecast horizons, on the contrary by selecting $MAPE$
measure the Proposed predictor $CP$ outperforms. Besides, results shows that
there is not a high significant variation of results with the selection of the
different kernels in Local Linear Regression.

####  5.2.4 Monthly pneumonia and influenza deaths

Monthly pneumonia and influenza deaths per $10.000$ people in the United
States for 11 years, $1968$ to $1978$.

This time series plotted in Figure 12 has $132$ observations, the first $84$
observations were used as training set and the last $24$ as a test set.

Figure 12: Monthly pneumonia and influenza deaths (1968$-$1978).

Figure 13 of auto-correlation at the appendix shows a seasonality of
approximately 12 months. In this line the predictor is based on a auto-
regressive model of order $p=12$, this is
$r(z_{k})=z_{k}=[y_{k-1}\;y_{k-2}\;...\;y_{k-12}]^{T}$.

Figure 13: Auto-correlation function of Monthly pneumonia and influenza deaths
time series.

Figure 14 shows the forecasts of the proposed predictor by forecasting
horizons with the hyper-parameters selected in both error measures.

Figure 14: Monthly pneumonia and influenza deaths time series predictions by
forecasting horizon in the test set.

Results of this time series are shown on Table 10 at the appendix. The hyper-
parameter $\gamma$ is selected in the training set where the error is minimum,
the value of $\gamma$ is inferred to perform forecasts in the test set as
shown in table 5.

Table 5: Monthly pneumonia and influenza deaths time series optimal gamma. Ahead |  $\gamma$_mape |  $\gamma$_smape
---|---|---
1.00 | 0.50 | 0.50
2.00 | 0.22 | 0.22
3.00 | 0.08 | 0.08

To sum up table 10 in a graphical way, figure 15 plots the error measures by
predictor and prediction horizon.

(a) SMAPE

(b) MAPE

Figure 15: Mean of test set errors by forecasting horizon in monthly pneumonia
and influenza deaths time series.

Results shows by forecasting horizon that the proposed predictor outperforms
in both error measures in the proposed forecast horizons. Besides, the results
shows that tricube kernel is a suitable option that outperforms between the
other selected non parametric methods for a short term forecast

####  5.2.5 Averaged results

This subsection averages the results shown in tables 7, 8, 9 and 10 attached
at the appendix. This results, averages all time series test error results in
two error measures selection criterion, depending on the selected criteria in
the training set the results could vary. Results shows that in average the
proposed predictor outperforms the selected methods in the different results.

(a) SMAPE

(b) MAPE

Figure 16: Mean of test set errors by forecasting horizon

Figures from 17 to 19 with a different aggregation level show an average by
grouping results in the different forecasting horizons shows that the proposed
predictor is outperforming the proposed methods and offering a competitive
alternative model.

Figure 17: Mean of SMAPE and MAPE results for 1 step-ahead forecasts Figure
18: Mean of SMAPE and MAPE results for 2 step-ahead forecasts Figure 19: Mean
of SMAPE and MAPE results for 3 step-ahead forecasts

###  5.3 Monthly electricity supplied

It is provided in this subsection a real- world dataset in order to perform
the proposed predictor. This is, the IEA provides monthly statistics with
timely and consistent oil, oil price, natural gas and electricity data for all
Organization for Economic Co-operation and Development member countries.

Countries submitted monthly data is adjusted proportionately to maintain
consistency with the most recent annual data for each generation source.

This time series is the electricity supplied for Spain from January of $2000$
to May of $2017$, this data consist in Indigenous production plus Imports
minus Exports. It includes transmission and distribution losses.

Figure 20 plots the aforementioned data set which has $221$ observations, the
first $181$ observations were used as training set and the last $40$ as a test
set.

(a) Original Time series

(b) Transformed Time series

Figure 20: Monthly electricity supplied in spain (2000$-$2017).

Figure 21 of auto-correlation at the appendix shows a seasonality of
approximately 12 months. In this line the predictor is based on a auto-
regressive model of order $p=12$, this is
$r(z_{k})=z_{k}=[y_{k-1}\;y_{k-2}\;...\;y_{k-12}]^{T}$.

Figure 21: Auto-correlation function of Monthly electricity supplied in spain.

The forecasts for the proposed predictor are in figure 22, these are plotted
by forecasting horizons with the hyper-parameters selected in both error
measures.

(a) MAPE selection criterion

(b) SMAPE selection criterion

Figure 22: Monthly electricity supplied time series predictions by forecasting
horizon in the test set.

Results of this time series are shown on Table 11 at the appendix. The hyper-
parameter $\gamma$ is selected in the training set where the error is minimum,
the value of $\gamma$ is inferred to perform forecasts in the test set as
shown in table 6.

Table 6: Monthly electricity supplied time series optimal gamma. Ahead |  $\gamma$_mape |  $\gamma$_smape
---|---|---
1.00 | 0.06 | 0.07
2.00 | 0.26 | 0.25
3.00 | 0.09 | 0.09

To sum up table 11 in a graphical way, figure 23 plots the error measures by
predictor and prediction horizon.

(a) SMAPE

(b) MAPE

Figure 23: Mean of test set errors by forecasting horizon in Monthly
electricity supplied time series.

Results from Figure 23 shows that the proposed predictor outperforms in both
error measures in the proposed forecast horizons. The proposed predictor is a
suitable option to consider in a real life problem.

##  6 Conclusions

A novel non-parametric Time Series forecasting method has been proposed. The
prediction is obtained by a weighted sum of past observations. A combination
of deterministic and stochastic assumptions are used to obtain an expression
of the outer bound of the prediction error. The weights are obtained solving a
convex optimization problem that minimizes the upper bound of the prediction
error. The method includes a tuning hyper-parameter. This hyper-parameter may
balance the deterministic and stochastic considered assumptions. By a cross-
validation scheme, a suitable hyper-parameter can be obtained. The performance
of the proposed predictor is exposed by some datasets.

##  7 Appendix

The following section contains the mathematical proofs, as well as research
results contained in tables and some summary plots of this tables to make it
easier for the reader to navigate through the document.

###  7.1 Mathematical derivations

Taking into account Assumption 1 and Definitions 3 and 2 the following
equalities can be inferred

$\displaystyle\hat{e}_{k}(\Psi)$ $\displaystyle=$ $\displaystyle
y_{k+1}-\hat{y}_{k+1}(\Psi)$ (18) $\displaystyle=$ $\displaystyle
y_{k+1}-\Psi^{T}b_{Y}$ (19) $\displaystyle=$ $\displaystyle
r(z_{k})^{T}\Phi_{k}-\Psi^{T}b_{Y}+e_{k}$ (20) $\displaystyle=$
$\displaystyle(A^{T}\Psi)^{T}\Phi_{k}-\Psi^{T}b_{Y}+e_{k}$ (21)
$\displaystyle=$ $\displaystyle\Psi^{T}(A\Phi_{k}-b_{Y})+e_{k}$ (22)
$\displaystyle=$
$\displaystyle\displaystyle\sum\limits_{j=1}^{k}\Psi_{j}(r(z_{j-1})^{T}\Phi_{k}-y_{j})+e_{k}$
(23) $\displaystyle=$
$\displaystyle-\displaystyle\sum\limits_{j=1}^{k}\Psi_{j}e_{j-1}+e_{k}.$ (24)

_QED_

###  7.2 Tables

Model | Ahead | $\gamma$ | tr_MAPE | te_MAPE | $\gamma$ | tr_SMAPE | te_SMAPE
---|---|---|---|---|---|---|---
CP | 1.00 | 0.12 | 15.75 | 13.04 | 0.12 | 15.46 | 14.97
LL1 | 1.00 | 1.57 | 15.19 | 14.65 | 1.84 | 14.34 | 15.69
LL2 | 1.00 | 1.55 | 14.82 | 17.11 | 1.62 | 13.77 | 14.13
LL3 | 1.00 | 1.89 | 16.42 | 15.47 | 2.33 | 14.04 | 16.98
NW1 | 1.00 | 1.08 | 17.69 | 27.26 | 1.08 | 13.82 | 21.05
NW2 | 1.00 | 1.02 | 17.36 | 28.38 | 1.02 | 13.48 | 21.90
NW3 | 1.00 | 1.24 | 17.91 | 27.16 | 1.26 | 13.98 | 21.09
CP | 2.00 | 0.14 | 17.36 | 14.05 | 0.14 | 16.88 | 16.58
LL1 | 2.00 | 1.58 | 13.69 | 18.89 | 1.58 | 13.56 | 16.45
LL2 | 2.00 | 1.54 | 14.22 | 22.02 | 1.69 | 13.24 | 17.71
LL3 | 2.00 | 1.67 | 14.79 | 18.37 | 2.09 | 14.38 | 18.44
NW1 | 2.00 | 1.12 | 17.06 | 29.65 | 1.14 | 13.29 | 22.91
NW2 | 2.00 | 0.98 | 16.71 | 30.40 | 0.98 | 13.05 | 23.22
NW3 | 2.00 | 1.26 | 17.33 | 29.41 | 1.28 | 13.40 | 22.61
CP | 3.00 | 0.00 | 17.43 | 15.45 | 0.00 | 16.50 | 18.66
LL1 | 3.00 | 1.64 | 16.06 | 24.99 | 2.15 | 15.59 | 21.13
LL2 | 3.00 | 1.61 | 17.89 | 25.89 | 2.01 | 15.37 | 22.14
LL3 | 3.00 | 1.78 | 16.49 | 25.33 | 2.34 | 15.72 | 20.72
NW1 | 3.00 | 1.12 | 17.16 | 32.13 | 1.14 | 13.23 | 24.05
NW2 | 3.00 | 0.98 | 16.88 | 32.80 | 0.98 | 13.13 | 24.09
NW3 | 3.00 | 1.24 | 17.55 | 31.12 | 1.30 | 13.37 | 23.75
Table 7: Airline passengers time series results

Model | Ahead | $\gamma$ | tr_MAPE | te_MAPE | $\gamma$ | tr_SMAPE | te_SMAPE
---|---|---|---|---|---|---|---
CP | 1.00 | 0.01 | 6.75 | 5.09 | 0.06 | 6.73 | 5.13
LL1 | 1.00 | 13.01 | 6.66 | 4.94 | 13.01 | 6.62 | 5.03
LL2 | 1.00 | 12.01 | 6.54 | 5.00 | 12.01 | 6.50 | 5.09
LL3 | 1.00 | 16.01 | 6.72 | 4.86 | 16.01 | 6.68 | 4.94
NW1 | 1.00 | 4.33 | 8.42 | 10.53 | 4.47 | 7.99 | 10.98
NW2 | 1.00 | 4.03 | 8.28 | 10.40 | 4.03 | 7.87 | 10.88
NW3 | 1.00 | 5.11 | 8.42 | 10.69 | 5.27 | 7.99 | 11.18
CP | 2.00 | 0.00 | 10.08 | 9.31 | 0.02 | 10.07 | 8.88
LL1 | 2.00 | 11.01 | 9.94 | 9.41 | 11.01 | 9.85 | 9.79
LL2 | 2.00 | 11.01 | 9.82 | 9.29 | 12.01 | 9.66 | 8.83
LL3 | 2.00 | 12.01 | 9.99 | 9.67 | 12.01 | 9.89 | 10.11
NW1 | 2.00 | 4.76 | 9.45 | 12.29 | 4.77 | 8.91 | 12.87
NW2 | 2.00 | 4.31 | 9.33 | 11.83 | 4.31 | 8.79 | 12.39
NW3 | 2.00 | 5.49 | 9.46 | 12.41 | 5.66 | 8.92 | 13.04
CP | 3.00 | 0.02 | 11.18 | 12.49 | 0.04 | 11.27 | 11.73
LL1 | 3.00 | 12.01 | 11.24 | 11.29 | 12.01 | 10.92 | 11.84
LL2 | 3.00 | 10.01 | 10.98 | 11.32 | 12.01 | 10.72 | 11.40
LL3 | 3.00 | 14.01 | 11.36 | 11.20 | 16.01 | 11.00 | 11.58
NW1 | 3.00 | 4.96 | 9.88 | 13.04 | 4.96 | 9.30 | 13.71
NW2 | 3.00 | 4.35 | 9.58 | 12.76 | 4.35 | 9.03 | 13.43
NW3 | 3.00 | 5.70 | 9.92 | 13.07 | 5.70 | 9.34 | 13.75
Table 8: Canadian Lynx time series results

Model | Ahead | $\gamma$ | tr_MAPE | te_MAPE | $\gamma$ | tr_SMAPE | te_SMAPE
---|---|---|---|---|---|---|---
CP | 1.00 | 0.00 | 7.33 | 6.70 | 0.00 | 7.26 | 6.89
LL1 | 1.00 | 50.00 | 7.25 | 7.29 | 50.00 | 7.23 | 7.09
LL2 | 1.00 | 30.00 | 7.17 | 6.98 | 30.00 | 7.18 | 6.74
LL3 | 1.00 | 60.00 | 7.26 | 7.31 | 60.00 | 7.24 | 7.10
NW1 | 1.00 | 13.00 | 8.92 | 10.06 | 13.00 | 8.92 | 9.36
NW2 | 1.00 | 13.00 | 9.09 | 11.68 | 13.00 | 9.04 | 10.79
NW3 | 1.00 | 15.00 | 8.86 | 9.91 | 15.00 | 8.85 | 9.24
CP | 2.00 | 0.03 | 11.16 | 8.71 | 0.00 | 10.92 | 9.58
LL1 | 2.00 | 40.00 | 10.75 | 10.13 | 40.00 | 10.74 | 9.50
LL2 | 2.00 | 40.00 | 10.71 | 10.27 | 40.00 | 10.69 | 9.66
LL3 | 2.00 | 50.00 | 10.74 | 10.11 | 50.00 | 10.74 | 9.49
NW1 | 2.00 | 14.00 | 10.51 | 14.07 | 14.00 | 10.45 | 12.84
NW2 | 2.00 | 12.00 | 10.60 | 14.68 | 13.00 | 10.52 | 13.78
NW3 | 2.00 | 16.00 | 10.49 | 13.69 | 16.00 | 10.44 | 12.53
CP | 3.00 | 0.00 | 12.42 | 9.63 | 0.00 | 12.25 | 10.34
LL1 | 3.00 | 50.00 | 12.26 | 11.05 | 50.00 | 12.14 | 10.22
LL2 | 3.00 | 40.00 | 12.18 | 11.08 | 40.00 | 12.07 | 10.25
LL3 | 3.00 | 50.00 | 12.24 | 10.98 | 50.00 | 12.14 | 10.15
NW1 | 3.00 | 14.00 | 11.65 | 16.46 | 14.00 | 11.53 | 14.87
NW2 | 3.00 | 13.00 | 11.74 | 17.64 | 13.00 | 11.57 | 15.81
NW3 | 3.00 | 15.00 | 11.61 | 15.00 | 16.00 | 11.52 | 14.56
Table 9: Monthly critical radio frequencies time series results Model | Ahead | $\gamma$ | tr_MAPE | te_MAPE | $\gamma$ | tr_SMAPE | te_SMAPE
---|---|---|---|---|---|---|---
CP | 1.00 | 0.50 | 10.30 | 11.89 | 0.50 | 9.70 | 12.06
LL1 | 1.00 | 1.47 | 10.47 | 18.44 | 1.48 | 10.31 | 17.19
LL2 | 1.00 | 1.41 | 10.69 | 19.08 | 1.41 | 10.56 | 17.60
LL3 | 1.00 | 1.61 | 10.80 | 17.63 | 1.61 | 10.68 | 16.38
NW1 | 1.00 | 0.92 | 10.87 | 20.04 | 0.92 | 10.96 | 18.46
NW2 | 1.00 | 0.92 | 11.38 | 21.41 | 0.92 | 11.45 | 19.72
NW3 | 1.00 | 0.92 | 10.58 | 18.84 | 0.92 | 10.69 | 17.31
CP | 2.00 | 0.22 | 12.81 | 15.14 | 0.22 | 12.04 | 15.85
LL1 | 2.00 | 1.95 | 13.99 | 23.72 | 1.48 | 13.95 | 21.40
LL2 | 2.00 | 1.90 | 13.54 | 23.27 | 1.90 | 13.77 | 21.02
LL3 | 2.00 | 1.61 | 14.70 | 22.66 | 1.55 | 14.08 | 19.38
NW1 | 2.00 | 1.01 | 11.22 | 22.56 | 0.96 | 11.24 | 20.07
NW2 | 2.00 | 0.95 | 11.65 | 23.46 | 0.95 | 11.67 | 21.50
NW3 | 2.00 | 0.92 | 10.84 | 20.67 | 0.92 | 10.78 | 18.73
CP | 3.00 | 0.08 | 13.62 | 15.11 | 0.08 | 12.97 | 15.71
LL1 | 3.00 | 1.91 | 13.52 | 23.40 | 1.49 | 13.38 | 19.94
LL2 | 3.00 | 1.88 | 13.09 | 23.12 | 1.88 | 13.11 | 21.08
LL3 | 3.00 | 2.00 | 13.88 | 23.09 | 1.63 | 13.43 | 19.07
NW1 | 3.00 | 0.92 | 11.08 | 21.20 | 0.92 | 11.17 | 19.53
NW2 | 3.00 | 0.92 | 11.65 | 21.58 | 0.92 | 11.84 | 19.99
NW3 | 3.00 | 0.94 | 10.92 | 20.81 | 0.94 | 10.96 | 19.09
Table 10: Monthly pneumonia and influenza deaths time series results Model | Ahead | $\gamma$ | tr_MAPE | te_MAPE | $\gamma$ | tr_SMAPE | te_SMAPE
---|---|---|---|---|---|---|---
CP | 1.00 | 0.06 | 13.88 | 21.81 | 0.07 | 14.16 | 23.88
LL1 | 1.00 | 1.90 | 17.08 | 30.84 | 2.00 | 14.52 | 24.42
LL2 | 1.00 | 1.87 | 16.93 | 30.80 | 1.50 | 14.40 | 24.14
LL3 | 1.00 | 2.00 | 17.14 | 30.86 | 2.00 | 14.62 | 24.38
NW1 | 1.00 | 0.68 | 17.76 | 31.16 | 0.73 | 15.00 | 25.00
NW2 | 1.00 | 0.65 | 17.58 | 31.60 | 0.65 | 15.04 | 24.65
NW3 | 1.00 | 0.78 | 17.81 | 31.21 | 0.81 | 14.98 | 24.76
CP | 2.00 | 0.26 | 15.22 | 24.32 | 0.25 | 15.23 | 26.22
LL1 | 2.00 | 2.00 | 18.01 | 35.40 | 2.00 | 15.94 | 27.67
LL2 | 2.00 | 1.98 | 17.77 | 35.36 | 1.85 | 15.85 | 27.53
LL3 | 2.00 | 2.00 | 18.18 | 35.46 | 2.00 | 16.09 | 27.71
NW1 | 2.00 | 0.70 | 19.01 | 35.19 | 0.72 | 16.14 | 27.07
NW2 | 2.00 | 0.68 | 18.96 | 37.16 | 0.70 | 16.24 | 27.75
NW3 | 2.00 | 0.78 | 19.04 | 34.48 | 0.80 | 16.19 | 26.74
CP | 3.00 | 0.09 | 15.95 | 21.80 | 0.09 | 16.13 | 24.51
LL1 | 3.00 | 1.96 | 18.62 | 35.73 | 1.92 | 16.30 | 26.50
LL2 | 3.00 | 1.79 | 18.33 | 35.24 | 1.59 | 16.02 | 26.80
LL3 | 3.00 | 2.00 | 18.76 | 35.87 | 2.00 | 16.48 | 26.58
NW1 | 3.00 | 0.73 | 19.92 | 37.01 | 0.73 | 16.87 | 27.81
NW2 | 3.00 | 0.65 | 19.24 | 38.11 | 0.65 | 16.70 | 28.31
NW3 | 3.00 | 0.83 | 20.01 | 36.74 | 0.85 | 16.92 | 27.79
Table 11: Monthly electricity supplied in spain time series results

## References

  * [1] G. Box and G. Jenkins, Time Series Analysis: Forecasting and Control.  San Francisco, CA: Holden-Day, 1976.
  * [2] J. Hamilton, Time series analysis.  Princeton, NJ: Princeton Univ. Press, 1994.
  * [3] H. Tong, Threshold Models in Nonlinear Time Series Analysis, vol. 21 of Lecture Notes in Statistics.  Heidelberg: Springer, 1983.
  * [4] V. Haggan and T. Ozaki, “Modeling nonlinear vibrations using an amplitude-dependent autoregressive time series model,” Biometrika, vol. 68, pp. 186–196, 1981.
  * [5] K. Chang and H. Tong, “On estimating thresholds in autoregressive models,” Journal of Time Series Analysis, vol. 7, pp. 179–190, 1986.
  * [6] Y. Truong, A nonparametric framework for time series analysis.  New Directions in Time Series Analysis, New York: Springer, 1993.
  * [7] E. A. Nadaraya, “On estimating regression,” Theory of Probability & Its Applications, vol. 9, no. 1, pp. 141–142, 1964.
  * [8] G. S. Watson, “Smooth regression analysis,” Sankhyā Ser., vol. 26, pp. 359–372, 1964.
  * [9] W. Härdle, Applied nonparametric regression.  No. 19 in Econometric Society monographs, Cambridge u. a.: Cambridge University Pr., 1990.
  * [10] J. Fan and I. Gijbels, Local polynomial modelling and its applications.  No. 66 in Monographs on statistics and applied probability series, London [u.a.]: Chapman and Hall, 1996.
  * [11] T. Hastie and R. Tibshirani, Generalized Additive Models.  Monographs on Statistics and Applied Probability, Chapman and Hall, 1990.
  * [12] W. Härdle, H. Lütkepohl, and R. Chen, “A review of nonparametric time series analysis,” International Statistical Review, vol. 65, no. 1, pp. 49–73, 1997.
  * [13] J. Fan and Q. Yao, Nonlinear Time Series: Nonparametric Methods and Parametric Methods.  Springer Series in Statistics, New York: Springer, 2003.
  * [14] J. Gao, Nonlinear Time Series: Semiparametric and Nonparametric Methods.  Chapman and Hall/CRC, 2007.
  * [15] J. G. Gooijer and A. Gannoun, “Nonparametric conditional predictive regions for time series,” Computational Statistics & Data Analysis, vol. 33, no. 3, pp. 259 – 275, 2000.
  * [16] Y. Yin and P. Shang, “Forecasting traffic time series with multivariate predicting method,” Applied Mathematics and Computation, vol. 291, pp. 266 – 278, 2016.
  * [17] C. Bergmeir, R. J. Hyndman, and B. Koo, “A note on the validity of cross-validation for evaluating autoregressive time series prediction,” Computational Statistics & Data Analysis, vol. 120, pp. 70 – 83, 2018.
  * [18] E. Mangalova and O. Shesterneva, “Sequence of nonparametric models for gefcom 2014 - probabilistic electric load forecasting,” International Journal of Forecasting, vol. 32, no. 3, pp. 1023 – 1028, 2016.
  * [19] R. J. Hyndman, M. L. King, I. Pitrun, and B. Billah, “Local linear forecasts using cubic smoothing splines,” Australian & New Zealand Journal of Statistics, vol. 47, no. 1, pp. 87–99, 2005.
  * [20] M. Milanese, J. Norton, H. Piet-Lahanier, and E. Walter, Bounding Approaches to System Identification.  Plenum Press, New York, 1996.
  * [21] J. M. Bravo, T. Alamo, M. Vasallo, and M. E. Gegúndez, “A general framework for predictors based on bounding techniques and local approximation,” IEEE Transactions on Automatic Control, vol. 62, pp. 3430–3435, July 2017.
  * [22] J. Roll, A. Nazin, and L. Ljung, “Nonlinear system identification via direct weight optimization,” Automatica, vol. 41, no. 3, pp. 475–490, 2005.
  * [23] S. Boyd and L. Vandenberghe, Convex Optimization.  Cambridge University Press, 2004.
  * [24] J. S. Armstrong, Long-range Forecasting: From Crystal Ball to Computer, vol. 2.  Wiley, 1985.
  * [25] R. J. Hyndman and A. B. Koehler, “Another look at measures of forecast accuracy,” International Journal of Forecasting, vol. 22, no. 4, pp. 679 – 688, 2006.
  * [26] M. J. A. M. W. Campbell, “A survey of statistical work on the mackenzie river series of annual canadian lynx trappings for the years 1821?1934 and a new analysis,” Journal of the Royal Statistical Society series, vol. A, 140, pp. 411–431, 1977.
  * [27] G. Zhang, “Time series forecasting using a hybrid arima and neural network model,” Neurocomputing, vol. 50, pp. 159–175, jan 2003.
  * [28] Newton, “Monthly critical radio frequencies in washington,” 1988.  data retrieved from datamarket.
  * [29] T. S. Rao and M. Gabr, An Introduction to Bispectral Analysis and Bilinear Time Series Models, vol. 24.  New York: Springer-Verlag, 1984.