chapter [0em] section[1.72em]. *[1.2pc].

# OFTER: An online pipeline for time series forecasting

Nikolas Michael† Corresponding author. Email<EMAIL_ADDRESS>${\dagger}$Department of Statistics, University of Oxford, 24-29 St Giles’,
Oxford OX1 3LB, UK
${\ddagger}$Mathematical Institute, University of Oxford, Andrew Wiles
Building, Woodstock Rd, Oxford OX2 6GG
$\S$ Oxford-Man Institute of Quantitative Finance, University of Oxford,
Oxford, UK
$\odot$ The Alan Turing Institute, John Dodson House, 96 Euston Rd, London NW1
2DB  Mihai Cucuringu†‡⊙§ and Sam Howison‡ ${\dagger}$Department of Statistics,
University of Oxford, 24-29 St Giles’, Oxford OX1 3LB, UK
${\ddagger}$Mathematical Institute, University of Oxford, Andrew Wiles
Building, Woodstock Rd, Oxford OX2 6GG
$\S$ Oxford-Man Institute of Quantitative Finance, University of Oxford,
Oxford, UK
$\odot$ The Alan Turing Institute, John Dodson House, 96 Euston Rd, London NW1
2DB

(April 2023)

###### Abstract

We introduce OFTER, a time series forecasting pipeline tailored for mid-sized
multivariate time series. OFTER utilizes the non-parametric models of
k-nearest neighbors and Generalized Regression Neural Networks, integrated
with a dimensionality reduction component. To circumvent the curse of
dimensionality, we employ a weighted norm based on a modified version of the
maximal correlation coefficient. The pipeline we introduce is specifically
designed for online tasks, has an interpretable output, and is able to
outperform several state-of-the art baselines. The computational efficacy of
the algorithm, its online nature, and its ability to operate in low signal-to-
noise regimes, render OFTER an ideal approach for financial multivariate time
series problems, such as daily equity forecasting. Our work demonstrates that
while deep learning models hold significant promise for time series
forecasting, traditional methods carefully integrating mainstream tools remain
very competitive alternatives with the added benefits of scalability and
interpretability.

Keywords Multivariate time series forecasting; non-parametric regression;
dimensionality reduction; maximal correlation; rank-one update; equity price
and volume forecasting.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries and related works
    1. 2.1 K-nearest neighbors and Generalized Regression Neural Network
    2. 2.2 Dimensionality reduction
    3. 2.3 Maximal correlation
  3. 3 Methodology
    1. 3.1 Online update of the embedding
    2. 3.2 Choice of scaling factor
    3. 3.3 Overcoming the shortcomings of kNN and GRNN
    4. 3.4 Algorithm outline
  4. 4 Empirical results
    1. 4.1 Baseline models
    2. 4.2 Synthetic data sets
    3. 4.3 Real data sets
    4. 4.4 Numerical results
  5. 5 Applications to feature importance and outlier detection.
    1. 5.1 Feature importance
    2. 5.2 Outlier detection
  6. 6 Conclusion and future work
  7. A A convexifying approach for solving the Secular equation
  8. B Brief summary of baseline methods
  9. C Robustness analysis

##  1 Introduction

Time series are temporally structured data sets that appear in most
disciplines, ranging from physics to the social sciences. For forecasting
tasks, the temporal structure can be leveraged to improve performance but also
to avoid look-ahead biases in the training phase. For univariate time series,
parametric linear time series models such as the Autoregressive Integrated
Moving Average (ARIMA) and the Autoregressive Conditional Heteroscedasticity
(ARCH) have attained sufficiently good results in most applications of
interest. While these algorithms can be extended to the setting of
multivariate time series, the rapid rise in computing power in recent decades
has rendered machine learning techniques an attractive alternative to
traditional time series methods. An early comparison between different machine
learning models was [1], where the authors concluded that the choice of method
had a significant effect in performance.

Deep learning models have been a popular choice for multivariate time series
forecasting, including Long-Short Term Memory (LSTM) [2], Bidirectional LSTM
[3] and Gated Recursion Unit [4] (GRU). Various expansions of such deep
learning models have been proposed in recent years, such as DeepAR [5] and
Temporal Fusion Transformers [6], which generally achieve state-of-the-art
performances across various domains. However, such approaches often employ
thousands of trainable parameters, which can be difficult to train with a
limited number of observations, thus leading to degraded out-of-sample
forecasting performance. Despite such advancements in deep learning, simple
regression-based models are still widely used and sought after. This is
especially true for mid-low frequency financial tasks, such predicting daily
returns, volumes, and realized volatilities, partly due to the complexity,
unpredictability and constantly shifting dynamics of the financial markets.
For example, Ordinary Least Squares (OLS) is still widely used for both
forecasting and interpretive statistical analysis due to its simplicity,
closed-form solution, and the fact that it requires no choice of hyper-
parameters. Its regularized versions, such as the LASSO and the Elastic Net
are some of the most widely utilized regression models in finance, as they
require only a few hyperparameters.

On the other hand, various recent lines of work have demonstrated the merit of
combining traditional time series models and machine learning models, such as
[7, 8, 9, 10, 11]. In this paper, we consider one of the simplest types of
hybrid models, which builds on the results of a univariate ARIMA model and
enhances it by performing a novel regression method on the residuals of an
ARIMA model.

##### Non-parametric models and dimensionality reduction

Two simple models that have been relatively overlooked by the machine learning
community, are k-Nearest neighbors (kNN) and Generalized Regression Neural
Networks (GRNN), both of which are simple, non-parametric, fast methods
requiring only a single hyperparameter. This mainly stems from the curse of
dimensionality, which affects both techniques. GRNNs have been utilized
primarily in engineering applications, such as [12] for logo recognition, [13]
for forecasting feeding patterns, and [14] for modeling the pile-bearing
capacity of soil. An overview of its different applications for dynamical
systems has been pursued by [15]. In this work, we aim to circumvent the curse
of dimensionality that these two methods suffer from by using the following
two techniques. Firstly, instead of the Euclidean norm employed by kNN and
GRNN, we consider a thresholded weighted norm that is driven by the feature’s
predictive potential. Secondly, before employing these well-known models, we
consider a dimensionality reduction component, which is by far the most common
remedy for the curse of dimensionality.

Dimensionality reduction is one of the main unsupervised machine learning
paradigms, with the past couple of decades having witnessed seminal
contributions in terms of non-linear dimensionality reduction methods, such as
locally linear embedding (LLE) [16], ISOMAP [17], Hessian LLE [18], Laplacian
eigenmaps [19], and diffusion maps [20]. In more recent years, deep learning
models have also been increasingly used as a means of dimensionality
reduction, such as the class of Autoencoders [21], which can be seen as a non-
linear extension of PCA. Performing dimensionality reduction before employing
a kNN technique has shown a significant improvement in performance [22].
Partial Least Squares is a popular method of incorporating information from
the target into the initial features, during the embedding stage. More
recently, [23] achieved improvements in the performance of PCA by scaling each
factor based on its predictive slope with the target.

##### Uses in financial applications

Examples of dimension reduction use cases in time series forecasting tasks
include [24], which used a truncated diffusion map to construct a global
manifold for the slow dynamics of a reduced kinetic model, and [25] which
utilizes manifold learning techniques in conjunction with various supervised
learning models to build a profitable strategy for the Forex market. Along the
same lines, [26] proposed a hybrid approach integrating the Midas-RV model and
PCA to achieve better results in forecasting the volatility of crude oil
futures in the Chinese market. One financial application we consider in this
paper is the forecasting of daily traded volumes. Being able to accurately
predict next-day volumes is of interest to practitioners and especially to
market makers. Previous works include [27], who showed how search engine data
can be leveraged to forecast future trading volumes. More similarly to our
work is [28], who utilizes Support Vector Machines to predict the daily volume
percentages with promising performance, in the spirit of [29] who combine SVM
and LSTM after residualizing with an AR model.

##### Main contribution: OFTER \- a novel integrated pipeline

We consider the task of one-step time series forecasting and propose an
approach that integrates ARIMA models, dimensionality reduction, and simple
non-parametric methods into a single pipeline, which we denote as Online
Forecasting via Temporal Embedding Regression (OFTER). To alleviate the curse
of high dimensionality and enhance our dimensionality reduction component, we
conduct feature selection on the embedding features based on the concept of
maximal correlation. We apply OFTER to the financial tasks of forecasting
next-day volume and returns by building an embedding based on various features
of 10 Exchange Traded Funds (ETFs). Through an extensive set of numerical
simulations on both synthetic and real data, we demonstrate that OFTER
provides a powerful alternative to other mainstream machine learning models
mid-sized data sets, even performing better than established time series-
tailored deep learning models. In particular, it can work effectively in low
signal-to-noise regimes for mid-sized data sets, while still being
computationally efficient.

Furthermore, our proposed methodology provides the following additional
benefits.

  * •

Once an embedding has been computed, it is straightforward for the algorithm
to perform out-of-sample extension and incorporate newly arrived time
observations; therefore, OFTER can naturally be construed as an online
algorithm that operates efficiently in an online setting.

  * •

Certain configurations of our method fall under the unsupervised learning
setting, as the embedding is independent of the target. In other words, the
same embedding can potentially be employed to forecast different target
variables of interest, and thus OFTER is amenable to extensions to the multi-
target setting.

  * •

Our method can further provide an avenue for change-point detection and
anomaly detection tasks, which can be seen as natural byproducts. A previous
or newly arrived observation without any local neighboring points in the
embedding, could indicate either an extreme measurement or a change-point in
the time series.

##### Paper Structure

The rest of this paper is structured as follows. In Section 2 we give an
overview of some of the essential concepts used in the proposed pipeline,
including the kNN and GRNN algorithms, principal component analysis, and
maximal correlation. In Section 3, we introduce our proposed methodology,
along with its main components and their integration into an algorithm. In
Section 4 we analyze the performance of OFTER in both synthetic experiments
and real-world data. We also discuss the additional properties of OFTER in
Section 5. Finally, we conclude in Section 6.

##  2 Preliminaries and related works

Consider a $d$-dimensional time series represented as a matrix
$X=\left(x_{1},\dots,x_{T}\right)\in\mathbb{R}^{T\times d}$ with rows denoting
observations measures along $T$ time steps, which we refer to the feature
matrix. Assume for the rest of this paper that we are interested in
forecasting a single univariate target variable, which we denote by
$y=\left(y_{1},\dots,y_{T}\right)\in\mathbb{R}^{T}$ and refer to as the target
variable. Furthermore, denote by $x^{\prime}\in\mathbb{R}^{d}$ and
$y^{\prime}\in\mathbb{R}$, the value of the features and the target variable,
respectively, associated with a new observation. We denote by $e_{d}$ the all-
ones vector of length $d$. In terms of matrix indexing, we generally denote
the (i,j) entry of $X$ by $X_{i,j}$, the submatrix of $X$ including rows
starting from $t_{1}$ up to $t_{2}$ by $X_{t_{1}:t_{2}}$ , i.e
$\left(x_{t_{1}},\dots,x_{t_{2}}\right)$ and by $X_{t_{1}:t_{2},j}$ the
$j^{th}$ column of $X_{t_{1}:t_{2}}$.

###  2.1 K-nearest neighbors and Generalized Regression Neural Network

K-Nearest Neighbors (kNN) is a non-parametric supervised learning algorithm
that dates back to the work of Fix and Hodges in 1951 [30], and has grown to
become one of the most popular machine learning techniques. The main idea is
that a similarity of values in the feature space implies a similarity in the
values for the target variable. Denote by $x_{(1)},\dots,x_{(T)}$ a reordering
of the time series entries, such that

$ ||x_{(1)}-x^{\prime}|| \leq\dots\leq||x_{(T)}-x^{\prime}||,$ (1)

where $||.||$ denotes the Euclidean norm. The k-nearest neighbors set for the
new entry $x$ is then defined as

$\textrm{NN}_{k}(x^{\prime})=\\{ i\in\\{1,\dots,T
\\}:||x_{(i)}-x^{\prime}||\leq||x_{(k)}-x^{\prime}||\\}. $ (2)

Based on $\textrm{NN}_{k}(x^{\prime})$, the forecast for $y^{\prime}$ is
computed as

$ \hat{y}=\frac{1}{k}\sum_{i\in\textrm{NN}_{k}(x^{\prime})}y_{i}.$ (3)

A closely related method is the Generalized Regression Neural Network (GRNN),
introduced in [31]. When forecasting $y^{\prime}$, GRNN attaches a weight to
each of the previous time series entries, determined by their distance to
$x^{\prime}$. In particular, define the weight $w_{i}$ of the $i^{th}$
observation as

$ w_{i}=\frac{\exp\left(\frac{-||x_{i}-x^{\prime}||^{2}}{h}
\right)}{\sum\limits_{j=1}^{T}\exp\left(\frac{-||x_{j}-x^{\prime}||^{2}}{h}
\right)},$ (4)

where $h$ is a bandwidth hyper-parameter to be optimized. The prediction for
the response of a new observation $\hat{y}$ is then simply the weighted
average of the responses of the training observations, defined as

$ \hat{y}=\sum_{i=1}^{T}w_{i}y_{i}.$ (5)

Consider the observations $(x_{1},y_{1}),\dots,(x_{T},y_{T})$, to be
realizations of an $\mathbb{R}^{d+1}$ random vector, and let $f$ denote the
corresponding true underlying density of the random variable. The theoretical
motivation of GRNN stems from the Parzen estimators [31], used to estimate
unknown probability density functions $f$ as

$ \hat{f}_{h}(x^{\prime},y^{\prime})=\frac{1}{T}\frac{1}{\left(\pi
h\right)^{\frac{d+1}{2}}}\sum_{i=1}^{T}\textrm{exp}\left[-\frac{(x^{\prime}-x_{i})^{T}(x^{\prime}-x_{i})}{h}\right]\textrm{exp}\left[-\frac{(y^{\prime}-y_{i})^{2}}{h}\right].$
(6)

If $h$ varies with $T$ such that $\lim_{T\to\infty}h(T)=0$ and
$\lim_{T\to\infty}Th^{d}(T)=0$, these estimators are shown to be consistent,
i.e

$plim_{T\to\infty}\hat{f}_{h(T)}(x^{\prime},y^{\prime})=f(x^{\prime},y^{\prime}),\
\ \forall(x^{\prime},y^{\prime})\in\mathbb{R}^{d+1}. $ (7)

Informally, [31] argues that if $\hat{f}(x,y)\approx f(x,y)$, then

$\hat{y}=\sum_{i=1}^{T}w_{i}y_{i}=\int
y^{\prime}\frac{\hat{f}(x,y^{\prime})}{\int\hat{f}(x,y)dy}dy^{\prime}\approx\int
y^{\prime}\frac{f(x,y^{\prime})}{\int f(x,y)dy}dy^{\prime}=
\mathbb{E}\left[y^{\prime}|x^{\prime}\right].$ (8)

We caution the reader that, despite GRNN bearing “Neural Network” in its name,
it is quite different from other popular variants. Specifically, aside from
the width hyper-parameter, there are no ad-hoc hyperparameters to be chosen or
optimized. Furthermore, it is a single-pass algorithm, as the weights are not
optimized. Thus, unlike other neural network structures, it is much faster and
easier to train but lacks the expressive power of neural networks.

###  2.2 Dimensionality reduction

Both kNN and GRNN suffer from the so-called _curse of dimensionality_ which
refers to the phenomenon that for certain methods, the amount of data points
required to estimate an unknown function with a certain accuracy, grows
exponentially with the dimensionality of the data set.111For a common
illustration of this problem, consider points that lie uniformly on a
$d$-dimensional hypercube of side length $(1+\epsilon)r$. In this setup, only
a proportion of $(1+\epsilon)^{-d}$ data points will lie in the central inner
hypercube of side length $r$, with the rest of the points lying in the outer
shell. This demonstrates that the distance between data points becomes
unsuitable for a meaningful description of the structure of high-dimensional
data sets To alleviate this issue, kNN or GRNN is often applied after
projecting the data points from $\mathbb{R}^{d}$ into a lower-dimensional
space $\mathbb{R}^{p}$ ($p\leq d$), while somewhat retaining the structure of
the original data set. By far the most popular such method is Principal
Component Analysis (PCA), which applies a linear transformation to the
original features, chosen to maximize the variance of the projected features
in the newly reduced space while ensuring the resulting features remain
uncorrelated with each other.

Denote by $C(X)$ the covariance matrix of $X$, and its eigen-decomposition by
$C(X)=U\Lambda U^{T}$, where $\Lambda$ is a diagonal matrix capturing the
$\lambda_{1}\geq\dots\geq\lambda_{d}\geq 0$ eigenvalues of $C(X)$, and the
columns of $U$ are given by the $u_{1},\dots,u_{d}$ orthonormal eigenvectors
(in the same order as their corresponding eigenvalues). The principal
component scores are given by the columns of $XU$, with the variance of each
of the new components given by its corresponding eigenvalue (component $i$ has
variance $\lambda_{i}$). Selecting the first $p$ components retains a
proportion of the total variance in the feature matrix equal to
$\sum_{i=1}^{p}\lambda_{i}/\sum_{i=1}^{d}\lambda_{i}$. Given the proportion of
variance to retain $\delta\in(0,1)$, we set

$p=\textrm{min}\left\\{i:\frac{\sum\limits_{k=1}^{i}\lambda_{k}}{\sum\limits_{k=1}^{d}\lambda_{k}}\geq\delta\right\\}.$
(9)

We remark that random-matrix theory based techniques could also be employed
for delimiting the informative part of the spectrum from the bulk of the
eigenvalues that correspond to the noise component [32]. Note that before
calculating $C(X)$, we ensure that $X$ is of full rank by removing any columns
of $X$ that result in rank deficiency.222Using a QR-decomposition algorithm,
we can express $X=QR$, where $Q$ is an orthogonal matrix and $R$ is an upper-
triangular matrix. Then we remove any feature $i$, where $|R_{i,i}|<\epsilon$
for some small $\epsilon>0$ Before applying PCA, we standardize $X$
(subtracting by each feature its mean and dividing by its standard deviation).
We shall denote the transformed features as $\tilde{X}=XU$.

###  2.3 Maximal correlation

Consider two univariate random variables denoted as $V_{1}$ and $V_{2}$.
Loosely stated, the correlation between $V_{1}$ and $V_{2}$ is a measure of
their association. Popular types of correlation include the Pearson, Spearman
and Kendall correlations. Another type is the maximal correlation (see e.g
[33]) defined as

$\underset{\phi,\psi\in\mathcal{F}}{\textrm{argmax}}\
\textrm{Cor}(\phi(V_{1}),\psi(V_{2})),$ (10)

where $\mathcal{F}$ is a suitable function class and Cor is the Pearson
correlation defined as

$\textrm{Cor}(V_{1},V_{2})=\frac{\mathbb{E}\left[V_{1}V_{2}\right]-\mathbb{E}\left[V_{1}\right]\mathbb{E}\left[V_{2}\right]}{\sqrt{\mathbb{E}\left[V_{1}^{2}\right]-\mathbb{E}\left[V_{1}\right]^{2}}\sqrt{\mathbb{E}\left[V_{2}^{2}\right]-\mathbb{E}\left[V_{2}\right]^{2}}}.$
(11)

A solution to the problem of determining maximal correlation can be computed
by optimizing the functions $\phi$ and $\psi$ iteratively. We consider a
related but simpler problem, where we aim to calculate

$\underset{\phi\in\mathcal{F}}{\textrm{argmax}}\
\textrm{Cor}(\phi(V_{1}),V_{2}).$ (12)

The reason for considering (12) instead of (10), is that in a regression
setting, we are only interested in transforming the feature, rather than the
target. By solving (10), one can transform the features in the embedding space
so as to be maximally correlated to the target.

We consider a solution in a restricted function class $\mathcal{F}$.333If
$\mathcal{F}$ is chosen to be the class of all functions, then the solution to
the above problem can easily be shown to have the trivial solution
$\Phi(V_{1})=\frac{V_{2}-\mathbb{E}\left[V_{2}\right]}{\textrm{Var}(V_{2})}$.
In particular, we consider the class of Bernstein polynomials (see e.g [34]),
which are defined as

$b_{\nu,n}(x)=\binom{\nu}{n}x^{\nu}(1-x)^{n-\nu},\quad\nu=0,1,\dots,n.$ (13)

The choice of this basis is justified via the Stone-Weistrass theorem.444In
its most usual form it says that for any real continuous function $f$ on the
interval $\left[a,b\right]$ and $\epsilon>0,\exists$ a polynomial function $g$
such that $\forall x\in\left[a,b\right]$, we have $|f(x)-g(x)|<\epsilon$.
Bernstein polynomials have been widely used, especially in the context of
Bezier curves [35], and B-splines [36], often arising in the development of
modeling tools in computer-aided geometric design.

For the solution of Equation 12, we assume that we observe vectors $v_{1}$ and
$v_{2}$, which capture $N$ observations of the random variables $V_{1}$ and
$V_{2}$, respectively. Denote by $\bar{v}_{1}$ and $\bar{v}_{2}$ the mean of
$v_{1}$ and $v_{2}$ respectively. Similarly denote by $\sigma_{v_{1}}$ and
$\sigma_{v_{2}}$ standard deviation of $v_{1}$ and $v_{2}$ respectively.
Consider the function basis $\left(\phi_{1},\dots,\phi_{K}\right)$, and define
the $N\times K$ matrix $\Phi$ s.t $\Phi_{i,j}=\phi_{j}(v_{1,i})$. We assume
$\Phi$ is centered, i.e $\Phi^{T}e_{N}=\textbf{0}$.

With the function basis defined, we consider the function class

$\mathcal{F}=\left\\{\sum_{k=1}^{K}c_{k}\phi_{k}:c_{1},\dots,c_{K}\in\mathbb{R}\right\\},$
(14)

so that for $f\in\mathcal{F}$,
$f(v_{1,i})=\sum_{k=1}^{K}c_{k}\phi_{k}(v_{1,i}),$ and problem (12) reduces to
finding the suitable vector $c$.

To simplify calculations, we also require the transformed vector to be
standardized, i.e $c^{T}\Phi^{T}e_{N}=0,\frac{1}{N-1}c^{T}\Phi^{T}\Phi c=1$.
Thus, we aim to solve

$\underset{c\in\mathbb{R}^{K},\ c^{T}\Phi^{T}e_{N}=0,\ c^{T}\Phi^{T}\Phi
c=N-1}{\textrm{argmax}}\textrm{Cor}(\Phi c,v_{2}).$ (15)

Note that

$\textrm{Cor}(\Phi c,v_{2})=\frac{\textrm{Cov}(\Phi
c,v_{2})}{\sigma_{v_{2}}\sigma_{v_{1}}}=\frac{\frac{(\Phi
c)^{T}v_{2}}{N-1}-\left[\frac{(\Phi
c)^{T}e_{N}}{N}\right]\left[\frac{v_{2}^{T}e_{N}}{N}\right]}{\sigma_{v_{2}}\sigma_{v_{1}}},$
(16)

and thus, under the constraints of (15), we have

$\textrm{Cor}(\Phi c,v_{2})=\frac{(\Phi c)^{T}v_{2}}{(N-1)\sigma_{v_{2}}},$
(17)

for which a solution can be found by considering the Lagrangian

$L=\frac{c^{T}\Phi^{T}v_{2}}{(N-1)\sigma_{v_{2}}}-\lambda(c^{T}\Phi^{T}e_{N})-\mu(c^{T}\Phi^{T}\Phi
c-(N-1)).$ (18)

Extrimizers of the above Lagrangian satisfy

$\mu=\frac{c^{T}\Phi^{T}v_{2}}{2(N-1)^{2}\sigma_{v_{2}}},\textrm{ and}$ (19)
$\lambda=\frac{\bar{v}_{2}}{(N-1)\sigma_{v_{2}}}.$ (20)

Substituting back to Equation 18, and writing $z=\Phi^{T}v_{2}$,
$\bar{z}=z-\bar{v}_{2}\Phi^{T}e_{N}$ and $A=\frac{\Phi^{T}\Phi}{N-1}$, we
reduce the problem to solving

$\bar{z}=\left(c^{T}z\right)Ac.$ (21)

By considering a solution of the form $c=\tau A^{-1}\bar{z}$ for
$\tau\in\mathbb{R}$, we arrive at the solution

$\tau=\frac{1}{\sqrt{\bar{z}^{T}(A^{-1})^{T}z}}.$ (22)

Overall, the solution to (15) is given by

$c=\frac{A^{-1}\bar{z}}{\sqrt{\bar{z}^{T}(A^{-1})^{T}z}}=\frac{(\Phi^{T}\Phi)^{-1}\Phi^{T}\left(v_{2}-\bar{v}_{2}e_{N}\right)}{\sqrt{\left(v_{2}-\bar{v}_{2}e_{N}\right)^{T}\Phi\bar{z}^{T}\left[(\Phi^{T}\Phi)^{-1}\right]^{T}\Phi^{T}v_{2}}}.$
(23)

We will refer to the correlation between the target $y$ and the transformed
$x$ as the _one-sided maximal correlation_ (OSMC) of x and $y$, as

$\textrm{OSMC}_{K}(x,y)=\textrm{Cor}(\Phi c,y),$ (24)

where $\Phi$ and $c$ are depended on $x$ and the choice of the number of basis
functions used $K$.

##  3 Methodology

Figure 1 is a pictorial diagram of the main components of the OFTER pipeline,
each of which is described in detail in the remainder of this section.

Figure 1: Pictorial illustration of OFTER, where the main components of its
iterative part are displayed.

###  3.1 Online update of the embedding

In order for the algorithm to run in an online setting, the embedding needs to
be updated as new data arrives. The initial embedding will be calculated based
on the first $L_{0}<T$ entries, and then updated at every time step. However,
re-computing the embedding at every single iteration can be quite expensive
(in fact, of complexity $O(nd^{2}+d^{3})$ at each step). For PCA, there exists
a computationally efficient approach that updates $U$ and $\Lambda$. In
particular, if $C(X_{t})$ denotes the covariance matrix built from the first
$t$ data points, one can compute the eigenvectors of the updated covariance

$C(X_{1:t+1})=\frac{t-1}{t}C(X_{1:t})+\frac{1}{t}||x_{t+1}||^{2}\frac{x_{t+1}}{||x_{t+1}||}\frac{x_{t+1}^{T}}{||x_{t+1}||},$
(25)

by using the previously computed eigenvectors of $C(X_{1:t})$ (we assume
centered data for now). The online algorithm proposed in [37] achieves just
that, with a complexity of $O(p^{2}d)$ at each iteration (recall that $p$
denotes the chosen dimension of the embedding).

The above is an instance of a more general problem, which is the rank-one
update of an eigensystem [38]. Consider again the symmetric matrix $A$, whose
eigen-decomposition is known, with $\lambda_{1}\geq\dots\geq\lambda_{d}$
denoting its eigenvalues, and $u_{i}$ be the eigenvector corresponding to
$\lambda_{i},\forall i\in\\{1,\dots,d\\}$ such that $A=U\Lambda U^{T}$. The
task at hand is to evaluate the eigensystem of $A+\rho vv^{T}$, by using the
known eigensystem of $A$. As shown by [38], the eigenvalues of $A+\rho vv^{T}$
are the roots $\kappa_{1},\dots,\kappa_{d}$ of the _secular equation_

$1+\rho\sum_{i=1}^{d}\frac{z_{i}^{2}}{\lambda_{i}-\kappa}=0,\ \textrm{where}\
z=U^{T}v.$ (26)

Once the updated eigenvalues are computed, each updated eigenvector $u_{k}$
can be calculated as

$u_{k}=\frac{U\left(\Lambda-\kappa_{k}I\right)^{-1}z}{||U\left(\Lambda-\kappa_{k}I\right)^{-1}z||}.$
(27)

Going a step further, one could reduce computations by only updating the first
few eigenvalues. A good approximation of the solution is given by [39], who
shows how this can be achieved by solving

$1+\rho\left(\sum_{i=1}^{p}\frac{z_{i}^{2}}{\lambda_{i}-\kappa}+\frac{1-\sum\limits_{i=1}^{p}z_{i}^{2}}{\mu-\kappa}\right)=0.$
(28)

The solutions of this equation are shown to approximate the first $m$
eigenvalues, with an error of the order $O(\underset{m+1\leq i\leq
p}{\textrm{max}}|\lambda_{i}-\mu|)$. Note that $\mu$ is another parameter to
be chosen, and should be made so as to minimize $O(\underset{m+1\leq i\leq
p}{\textrm{max}}|\lambda_{i}-\mu|)$, for example
$\mu=\textrm{median}(\lambda_{m+1},\dots,\lambda_{p})$. The parameter $p$ will
have already been chosen at an earlier part of the algorithm using Equation 9,
and does not change during the update step. While there are various existing
techniques to solve equation Equation 28, we employ the convexifying method
described in [40], adopting the earlier work of [38].

Note that initially the feature matrix is standardized based on the initial
$L_{0}$ entries (i.e we standardize $X_{(1:L_{0})}$). Denote the mean vector
of $X_{1:t}$ by $\mu_{t}\in\mathbb{R}^{d}$, and the correctly centered
$X_{1:t}$ by $\bar{X}_{1:t}\in\mathbb{R}^{d}$ Assuming that $X_{(1:t-1)}$ has
been centered by $\mu_{t-1}$, we need to re-center based on $\mu_{t}$.
Following [37], the update takes the form
$\bar{X}_{1:t}=X_{1:t}+e_{t}(\mu_{t-1}-\mu_{t})^{T}$. It turns out that we can
adjust for the re-centering, via two more rank-one updates of the covariance
matrix, by expressing $\bar{X}^{T}_{1:t}\bar{X}_{1:t}$ in the following way.

$\displaystyle\bar{X}^{T}_{1:t}\bar{X}_{1:t}$
$\displaystyle=X_{1:t}^{T}X_{1:t}+X_{1:t}^{T}e_{N}(\mu_{t-1}-\mu_{t})^{T}+(\mu_{t-1}-\mu_{t})e_{N}^{T}+t(\mu_{t-1}-\mu_{t})(\mu_{t-1}-\mu_{t})^{T}$
$\displaystyle=X_{1:t}^{T}X_{1:t}+\begin{bmatrix}\mu_{t-1}-\mu_{t}&X_{1:t}^{T}e_{N}\end{bmatrix}\begin{bmatrix}t&1\\\
1&0\end{bmatrix}\begin{bmatrix}(\mu_{t-1}-\mu_{t})^{T}\\\
e_{N}^{T}X_{1:t}\end{bmatrix}$
$\displaystyle=X_{1:t}^{T}X_{1:t}+\begin{bmatrix}\mu_{t-1}-\mu_{t}&X_{1:t}^{T}e_{N}\end{bmatrix}\begin{bmatrix}\frac{\rho_{1}}{\sqrt{1+\rho_{1}^{2}}}&\frac{\rho_{2}}{\sqrt{1+\rho_{2}^{2}}}\\\
\frac{1}{\sqrt{1+\rho_{1}^{2}}}&\frac{1}{\sqrt{1+\rho_{2}^{2}}}\end{bmatrix}\begin{bmatrix}\rho_{1}&0\\\
0&\rho_{2}\end{bmatrix}\begin{bmatrix}\frac{\rho_{1}}{\sqrt{1+\rho_{1}^{2}}}&\frac{\rho_{2}}{\sqrt{1+\rho_{2}^{2}}}\\\
\frac{1}{\sqrt{1+\rho_{1}^{2}}}&\frac{1}{\sqrt{1+\rho_{2}^{2}}}\end{bmatrix}^{T}\begin{bmatrix}(\mu_{t-1}-\mu_{t})^{T}\\\
e_{N}^{T}X_{1:t}\end{bmatrix}$
$\displaystyle=X_{1:t}^{T}X_{1:t}+\begin{bmatrix}\frac{\rho_{1}(\mu_{t-1}-\mu_{t})+X_{1:t}^{T}e_{N}}{\sqrt{1+\rho_{1}^{2}}}&\frac{\rho_{2}(\mu_{t-1}-\mu_{t})+X_{1:t}^{T}e_{N}}{\sqrt{1+\rho_{2}^{2}}}\end{bmatrix}\begin{bmatrix}\rho_{1}&0\\\
0&\rho_{2}\end{bmatrix}\begin{bmatrix}\frac{\rho_{1}(\mu_{t-1}-\mu_{t})+X_{1:t}^{T}e_{N}}{\sqrt{1+\rho_{1}^{2}}}\\\
\frac{\rho_{2}(\mu_{t-1}-\mu_{t})+X_{1:t}^{T}e_{N}}{\sqrt{1+\rho_{2}^{2}}}\end{bmatrix}$
$\displaystyle=X_{1:t}^{T}X_{1:t}+\rho_{1}b_{1}b_{1}^{T}+\rho_{2}b_{2}b_{2}^{T},$

where in the above

$\rho_{1}=\frac{t-\sqrt{t^{2}+4}}{2},$ (29)
$\rho_{2}=\frac{t+\sqrt{t^{2}+4}}{2},$ (30)
$b_{1}=\frac{\rho_{1}(\mu_{t-1}-\mu_{t})+X_{1:t}^{T}e_{N}}{\sqrt{1+\rho_{1}^{2}}},$
(31)
$b_{2}=\frac{\rho_{2}(\mu_{t-1}-\mu_{t})+X_{1:t}^{T}e_{N}}{\sqrt{1+\rho_{2}^{2}}}.$
(32)

###  3.2 Choice of scaling factor

The hyperparameter $k$ of kNN in (5) and the hyperparameter $h$ in (4) are
both of crucial to the performance of their respective algorithms. Notice
that, as $h\to\infty$, $w_{i}\to 1/T$ (Equation 4) for all $i$. On the other
hand, as $h\to 0$, GRNN becomes equivalent to $1\textrm{NN}$ (kNN with $k=1$).
To demonstrate the effect of the scaling factor on the predictive task, we
consider the following toy model

$y_{t}=0.5+\cos\left(\frac{\pi}{64}t\right)+\epsilon_{t},\textrm{where
}\epsilon_{t}\sim N(0,\sigma^{2}).$ (33)

Figure 2: Visualization of the forecasted values when using GRNN with
different scaling parameters for the model (33). In a high signal-to-noise
problem like this, as $h$ becomes larger GRNN becomes almost perfectly capable
of forecasting future values.

In Figure 2 we demonstrate how the choice of the scaling factor affects the
nature of the obtained fit. For low values of $h$, the obtained fit is close
to giving constant fitted values equal to the mean, while for high scaling
factors the obtained fit matches the observed data points almost
perfectly.555Note that this is a pattern for this specific toy data set and
not a general pattern

To choose a suitable value of $h$, [31] proposes the use of the holdout
method, where the optimal parameter is found by calculating a loss function by
removing a sample point at a time and taking the average. More modern
applications use K-fold validation instead, and these methods can be modified
for time series data (see [41] for a discussion on the different approaches).

To remove any scale dependencies, we re-parameterize $h$ as

$h=\frac{\textrm{median}(\\{||x_{j}-x^{\prime}||\\}_{j=1}^{T})}{s}.$ (34)

Since this is an online algorithm, we could update the scaling parameter
accordingly. Instead of choosing a single bandwidth parameter from a candidate
set, we use each of the candidate bandwidths to fit a different model,
updating the performance of all models at each time point. The model with the
highest performance up to the previous time point is chosen for the forecast
of the current time point.666Alternatively, we can take a weighted average of
the forecasts of each model, but the performance was inferior.

The choice of the scaling factor can be performed by the use of a loss
function. If the fitted value corresponding to scaling parameter $s$ for day
$t$ is $\hat{y}^{s}_{t}$, the forecast $\hat{y_{t}}$ is calculated as

$\hat{y}_{t}=\sum_{s=1}^{|S|}\eta_{t}^{s}\hat{y}_{t}^{s},\ \textrm{where}\ \
\eta_{t}^{s}=\begin{cases}1,\ \
&R_{t}^{s}=\textrm{min}(R_{t}^{1},\dots,R_{t}^{S})\\\ 0,\ \
&\textrm{otherwise}\\\ \end{cases}.$ (35)

By $R_{t}^{s}$, we define the loss corresponding to parameter $s$, i.e
$R_{t}^{s}=\sum_{i=1}^{t-1}\textrm{Loss}(\hat{y}_{t}^{s},y_{t})$. The loss
function can be the Mean Squared Error (MSE) or Mean Absfolute Error (MAE) for
general applications. For the application of forecasting future returns, we
use the loss function of minus P&L (we define P&L in Section 4.3).

Similarly, a large of $k$ for kNN algorithm would produce similar forecasts to
using low values of $h$. Instead of using either GRNN or kNN, we use both
models, each with their own hyperparameter candidate set, and then select
based on performance, in a similar manner to (35). Aside from these two types
of models, we also consider a standard OLS model.

Denote by $S^{\textrm{GRNN}}$ the hyperparameter canditate for GRNN and
$S^{\textrm{kNN}}$ the hyperparameter canditate set for kNN. Denote by
$\hat{y}^{\textrm{GRNN},s}_{t}$, the forecast when using GRNN with parameter
$s$, $\hat{y}^{\textrm{kNN},k}_{t}$ the forecast when using kNN (with
parameter $k$), and $\hat{y}_{t}^{\textrm{OLS}}$ the forecast when using an
OLS model. Then define

$\begin{gathered}R^{\textrm{GRNN},s}_{t}=\sum_{i=1}^{t-1}\textrm{Loss}(\hat{y}^{\textrm{GRNN},s}_{i},y_{i}).\\\
R^{\textrm{kNN},k}_{t}=\sum_{i=1}^{t-1}\textrm{Loss}(\hat{y}^{\textrm{kNN},k}_{i},y_{i}).\\\
R^{\textrm{OLS}}_{t}=\sum_{i=1}^{t-1}\textrm{Loss}(\hat{y}^{\textrm{OLS}}_{i},y_{i}).\end{gathered}$
(36)

Based on the losses up to time point $t$, the minimum loss is calculated as

$R^{\textrm{optimal}}_{t}=\textrm{min}(\\{R^{\textrm{GRNN},s}_{t}:s\in
S^{\textrm{GRNN}}\\}\cup\\{R^{\textrm{kNN},k}_{t}:k\in
S^{\textrm{kNN}}\\}\cup\\{R_{t}^{\textrm{OLS}}\\}),$ (37)

and the number of models achieving the optimal performance as

$N^{\textrm{optimal}}_{t}=|\\{s:R^{\textrm{GRNN},sx}=R^{\textrm{optimal}}_{t}\\}|+|\\{k:R^{\textrm{kNN},k}=R^{\textrm{optimal}}_{t}\\}|+\mathbbm{1}({R_{t}^{\textrm{OLS}}=R^{\textrm{optimal}}_{t}}),$
(38)

where $\mathbbm{1}(A)$ is an indicator function taking the value $1$, if
condition $A$ is satisfied and $0$ otherwise. Then we reformulate (35), to
include the kNN and OLS models as

$\hat{y}_{t}=\sum_{s\in
S^{\textrm{GRNN}}}\eta_{t}^{\textrm{GRNN},s}\hat{y}^{\textrm{GRNN},s}_{t}+\sum_{k\in
S^{\textrm{kNN}}}\eta_{t}^{\textrm{kNN},k}\hat{y}_{t}^{\textrm{kNN},k}+\eta_{t}^{\textrm{OLS}}\hat{y}_{t}^{\textrm{OLS}},$
(39)

where

$\begin{gathered}\eta^{\textrm{GRNN},s}_{t}=\begin{cases}\frac{1}{N^{\textrm{optimal}}_{t}},\
\ &R^{\textrm{GRNN},s}_{t}=R^{\textrm{optimal}}_{t}\\\ 0,\ \
&\textrm{otherwise}\\\ \end{cases},\\\
\eta^{\textrm{kNN},k}_{t}=\begin{cases}\frac{1}{N^{\textrm{optimal}}_{t}},\ \
&R^{\textrm{kNN},k}_{t}=R^{\textrm{optimal}}_{t}\\\ 0,\ \
&\textrm{otherwise}\\\ \end{cases},\\\
\eta^{\textrm{OLS}}_{t}=\begin{cases}\frac{1}{N^{\textrm{optimal}}_{t}},\ \
&R^{\textrm{OLS}}_{t}=R^{\textrm{optimal}}_{t}\\\ 0,\ \ &\textrm{otherwise}\\\
\end{cases}.\end{gathered}$ (40)

###  3.3 Overcoming the shortcomings of kNN and GRNN

One downside of applying dimensionality reduction methods is that they are
agnostic of the target to be forecasted. In transforming the initial space in
order to retain a large proportion of the variability with a reduced number of
features, an initial feature that is already highly correlated with the target
but also weakly correlated with other features could end up having a
relatively small weight in the final forecasting step. To circumvent this
shortcoming, we employ two additional steps.

Firstly, we make sure that all features from the original feature space $X$
that have a high enough correlation with the target are retained in the
forecasting step. This way, we ensure that selecting the components that
capture most of the variability is not done at the expense of dropping
original features that provide predictive power on their own. Therefore, we
define the following set, after calculating the OSMC between each of the
original features and the target

$I^{\textrm{Original}}=\\{j:\textrm{OSMC}(X_{1:L_{0},j},y_{1:L_{0}})\geq
c_{\textrm{Original}}\\}.$ (41)

During the pre-processing part of the algorithm, we replace $\tilde{X}$ with
$(\tilde{X},X_{1:T,I^{\textrm{Original}}})$, with a number of features equal
to $p+|I^{\textrm{Original}}|$. Note that for mid-sized, low signal-to-error
data sets like the financial data set we use in Section 4.3,
$I_{\textrm{Original}}$ is usually null.

Secondly, certain features in the embedding space ought to have more influence
on the neighbor selection. This can be achieved by using a weighted metric
rather than the standard Euclidean norm. In particular, if $v$ represents the
vector of weights, the weighted Euclidean distance between $X_{t_{1}}$ and
$X_{t_{2}}$ is given by

$||\tilde{X}_{t_{1}}-\tilde{X}_{t_{2}}||_{v}=\sqrt{\sum_{j=1}^{p+|I^{\textrm{Original}}|}v_{j}(\tilde{X}_{t_{1},j}-\tilde{X}_{t_{2},j})^{2}}.$
(42)

The weights are based on the OSMC between each feature and the target

$v_{j}=\frac{c_{j}^{2}\mathbbm{1}(|c_{j}|\geq
c_{\textrm{min}})}{\sum\limits_{l=1}^{p+|I^{\textrm{Original}}|}c_{l}^{2}\mathbbm{1}(|c_{l}|\geq
c_{\textrm{min}})},$ (43)

where $c_{j}=\textrm{OSMC}(\tilde{X}_{1:L_{0},j},y_{1:L_{0}})\ \ \forall
j\in\\{1,\dots,p,p+1,\dots,|I^{\textrm{Original}}|\\}$.

###  3.4 Algorithm outline

To ensure stationarity of the target $y$, we first perform an ADF test with
significance level $p_{\textrm{ADF}}$, and if the null hypothesis is not
rejected, we difference the time series, and repeat until we reject the null
hypothesis. We then perform ACF and PACF tests to check if there are any
significant ARMA effects.777We use a significance level of 0.001. We consider
a specific lag effect significant if the confidence interval of either the ACF
or the PACF is above or below zero. If such effects are identified, we perform
a grid search over the orders of the AR and MA orders (p and q respectively)
to find the optimal ARIMA(p,r,q) model, where $r$ denotes the number of times
differencing was applied. Otherwise, we simply fit ARIMA(0,r,0).

Using this procedure, we decompose the target $y$ into a time series and a
residual component

$y=\underbrace{y^{\textrm{TS}}}_{\textrm{Time series
component}}+\underbrace{\epsilon}_{\textrm{Residual}}.$ (44)

Overall, the algorithm consists of three main steps: (1) computing an
embedding, (2) choosing features based on the OSMC, and (3) forecasting based
on kNN and GRNN. The complete algorithm is given in algorithm 2.

Input: $\tilde{X}$ (Feature Matrix), $y$ (Target), $L$ (Lookback Window
Length), $v$ (distance vector), $S^{\textrm{kNN}},S^{\textrm{GRNN}}$
(parameter candidate sets)

Output: $\hat{y}_{t}$

1

1: Compute $d_{j,t}=||\tilde{X}_{t,I}-\tilde{X}_{j,I}||_{v}$ for
$j=t-L,\dots,t-1$

2: for $s\in S^{\textrm{GRNN}}$ do

3: Set $h_{s}=\frac{\textrm{median}(\\{d_{j,t}\\}_{j=t-L}^{t-1})}{s}$

4: Compute
$w_{j,t}^{s}=\frac{\exp\left(\frac{-d_{j,t}^{2}}{h_{s}}\right)}{\sum\limits_{m=t-L}^{t-1}\exp\left(\frac{-d_{m,t}^{2}}{h_{s}}\right)}\forall
j\in\\{t-L,\dots,t-1\\}$

5: Compute
$\hat{y}^{\textrm{GRNN},s}=\sum\limits_{j=t-L}^{t-1}w_{j,t}^{s}y_{i}$

6: end for

7: for $k\in S^{\textrm{kNN}}$ do

8: Find the set $\textrm{NN}_{k}(\tilde{X}_{t,})$ of the $k$-nearest neighbors
of $\tilde{X}_{t}$, in $\tilde{X}_{(t-L):(t-1)}$

9: Compute
$\hat{y}^{\textrm{kNN},k}=\frac{1}{k}\sum\limits_{i\in\textrm{NN}_{k}(\tilde{X}_{t})}y_{i}$

10: end for

11: Compute $\hat{y}^{\textrm{OLS}}=\beta_{0}+\beta^{T}\tilde{X}_{t}$

12: Estimate the fitted value $\hat{y}_{t}$ using (39)

Algorithm 1 Forecasting for time point $t$

Input: $X$ (Feature Matrix), y (target), $L_{0}$ (Training Window Length), $L$
(Lookback Window Length), $c_{\textrm{min}}$ , $c_{\textrm{Original}}$
(correlation thresholds), $S^{\textrm{kNN}},S^{\textrm{GRNN}}$ (parameter
candidate sets), $\delta$ (proportion of variance to
retain),$p_{\textrm{ADF}}$ (p-value for ADF), LOSS (loss function)

Output: $\hat{y}$

1

1: Use the procedure described in Section 3.4, to arrive at the optimal
ARIMA($p$,$r$,$q$) model for $y$

2: Use the optimal model to obtain fitted values $y^{\textrm{TS}}$. If
ARIMA(0,0,0) is chosen set $y^{\textrm{TS}}=(0,\dots,0)$

3: Transform $X$ so as to make $X_{1:L_{0}}$ a full column-rank (by QR
decomposition) and standardized matrix

4: Obtain the first $p$ principal components of $X_{1:L_{0}}$, with $p$ chosen
via Equation 9, and set $U_{p}=(u_{1},\dots,u_{p})$ and
$\Lambda_{p}=(\lambda_{1},\dots,\lambda_{p})$

5: Obtain the $L_{0}\times p$ matrix $\tilde{X}=XU_{1:L_{0},1:p}$

6: Set
$I^{\textrm{Original}}=\\{i:\textrm{OSMC}(X_{1:L_{0},i},y_{1:L_{0}})\geq
c_{\textrm{Original}}\\}$

7: Replace $\tilde{X}$ by $(\tilde{X},X_{1:T,I^{\textrm{Original}}})$

28: Compute $R^{\textrm{GRNN},s}_{L_{0}}\ \forall s\in
S^{\textrm{GRNN}},R^{\textrm{kNN},k}_{L_{0}}\ \forall k\in S^{\textrm{kNN}}$,
and $R^{\textrm{OLS}}_{L_{0}}$ using (36)

39: Compute $\eta^{\textrm{GRNN},s}_{L_{0}}\ \forall s\in
S^{\textrm{GRNN}},\eta^{\textrm{kNN},k}_{L_{0}}\ \forall k\in
S^{\textrm{kNN}}$, and $\eta^{\textrm{OLS}}_{L_{0}}$ using (40)

10: for $j\leftarrow 1$ to $p+|I^{\textrm{Original}}|$ do

11: Compute $c_{j}=\textrm{OSMC}(\tilde{X}_{1:L_{0},j},y_{1:L_{0}})$

12: Set $v_{j}=\frac{c_{j}^{2}\mathbbm{1}(|c_{j}|\geq
c_{\textrm{min}})}{\sum_{l=1}^{p+|I^{\textrm{Original}}|}c_{l}^{2}\mathbbm{1}(|c_{l}|\geq
c_{\textrm{min}})}$

13: end for

14: Regress $y_{1:L_{0}}$ on $\tilde{X}_{1:L_{0},}$ using OLS regression, to
obtain the coefficient vector $\beta$ and the intercept $\beta_{0}$

15: for $t\leftarrow L_{0}+1$ to $D$ do

16: Compute $\hat{y}_{t}$ using algorithm 1

17: Update $U_{p}$ and $\Lambda_{p}$, using the secular equation method
(detailed in Section 3.1)

18: If t is divisible by 100, re-fit the OLS model using $\tilde{X}_{(t-L):t}$
and $y_{(t-L):t}$, to obtain $\beta$ and $\beta_{0}$

19: end for

20: Update $\hat{y}$ to $\hat{y}+y^{\textrm{TS}}_{(L_{0}+1):D}$

Algorithm 2 OFTER

##  4 Empirical results

We consider and compare four different configurations of OFTER, which we
denote as follows

  * •

OFTER (no dimensionality reduction and no feature transformation). This
entails setting $\tilde{X}=X$, using the Pearson correlation instead of the
OSMC, and skipping steps 4,5, 6, 7 and 17 of algorithm 2.

  * •

OFTER-DR (dimensionality reduction and no feature transformation). This is
almost the same as algorithm 2, but we use the Pearson correlation instead of
the OSMC.

  * •

OFTER-FT (feature transformation and no dimensionality reduction). This
entails setting $\tilde{X}=X$ as well as skipping steps 4,5, 6, 7 and 17 of
algorithm 2.

  * •

OFTER-DR-FT (feature transformation and dimensionality reduction). This is the
exact algorithm as given by algorithm 2.

Unless stated otherwise we use $p_{\textrm{ADF}}=0.05$. When performing
dimensionality reduction, we consider $\delta=0.9$. For the candidate sets we
use $S^{GRNN}=\\{0.001,0.005,0.01,0.05,0.1,0.5,1,5,10,50,100\\}$ and
$S^{kNN}=\\{1,2,3,5,10,15,20,30,50\\}$. In terms of features, we include those
corresponding to the current period up to a lag of $3$. Lastly, we set $L_{0}$
equal to $\left\lfloor 0.7T\right\rfloor$. Note that OFTER is generally robust
for different choices of hyper-parameters, as evidenced by our robustness
analysis reported in Appendix C.

###  4.1 Baseline models

We compare the performance of OFTER against that of several mainstream machine
learning models. We include the linear models OLS, PCR, LASSO, but also the
widely used SVM and Random Forests. We also include some of the most widely
used time series tailored deep learning models, such as LSTM, Bi-LSTM, and
GRU. Unlike OFTER, most of these methods are not naturally amenable to
efficient implementation in an online setting and thus cannot be updated. For
comparison purposes, we retrain each of these models every 100 steps. In the
case of LSTM, Bi-LSTM, and GRU, we re-fit the model with only five epochs,
starting from the previous weight matrix. For the sake of a fair comparison,
instead of following the embedding updating step, we also re-compute the
embedding of OFTER every 100 steps. We provide an overview of these methods,
as well as details of their implementation for this task, in Appendix B.

###  4.2 Synthetic data sets

Before testing OFTER on real life data sets, we first check whether it can
perform consistently well on some simple low dimensional toy models. For this
purpose, the three five-dimensional models used by [25], which are simplistic,
that are crafted to resemble electroencephalography signals, which is suitable
as a data set for this forecasting task.

$\textrm{M1:}\begin{cases}y_{t}^{(1)}=0.2y_{t-1}^{(1)}-0.4y_{t-1}^{(2)}+\epsilon_{t}^{(1)}\\\
y_{t}^{(2)}=-0.5y_{t-1}^{(1)}+0.15y_{t-1}^{(2)}+\epsilon_{t}^{(2)}\\\
y_{t}^{(3)}=-0.14y_{t-1}^{(2)}+\epsilon_{t}^{(3)}\\\
y_{t}^{(4)}=0.5y_{t-1}^{(1)}-0.25y_{t-1}^{(2)}+\epsilon_{t}^{(4)}\\\
y_{t}^{(5)}=0.15y_{t-1}^{(1)}+\epsilon_{t}^{(5)}\end{cases}$ (45)

The first data set (M1) shown in Equation 45 is generated by a set of five
linear equations.

$\textrm{M2:}\begin{cases}y_{t}^{(1)}=3.4y_{t-1}^{(1)}\left(1-\textrm{exp}\left(-\left(y_{t-1}^{(1)}\right)^{2}\right)\right)\textrm{exp}\left({-y_{t-1}^{(1)}}^{2}\right)+\epsilon_{t}^{(1)}\\\
y_{t}^{(2)}=3.4y_{t-1}^{(2)}\left(1-\textrm{exp}\left(-\left(y_{t-1}^{(2)}\right)^{2}\right)\right)\textrm{exp}\left({-y_{t-1}^{(2)}}^{2}\right)+0.5y_{t-1}^{(1)}y_{t-1}^{(2)}+\epsilon_{t}^{(2)}\\\
y_{t}^{(3)}=3.4y_{t-1}^{(3)}\left(1-\textrm{exp}\left(-\left(y_{t-1}^{(3)}\right)^{2}\right)\right)\textrm{exp}\left({-y_{t-1}^{(3)}}^{2}\right)+0.3y_{t-1}^{(2)}+0.5{y_{t-1}^{(1)}}^{2}+\epsilon_{t}^{(3)}\\\
y_{t}^{(4)}=0.5y_{t-1}^{(1)}-0.25y_{t-1}^{(2)}+\epsilon_{t}^{(4)}\\\
y_{t}^{(5)}=0.15y_{t-1}^{(1)}+\epsilon_{t}^{(5)}\end{cases}$ (46)

The second data set (M2) shown in Equation 46 is generated by a set of five
equations, some of which are non-linear. [25] based these equations on some
models appearing in [42].

$\textrm{M3:}\begin{cases}y_{t}^{(1)}=0.1y_{t-1}^{(1)}-0.6y_{t-3}^{(2)}+\epsilon_{t}^{(1)}\\\
y_{t}^{(2)}=-0.15y_{t-3}^{(1)}+0.8y_{t-3}^{(2)}+\epsilon_{t}^{(2)}\\\
y_{t}^{(3)}=-0.45y_{t-3}^{(2)}+\epsilon_{t}^{(3)}\\\
y_{t}^{(4)}=0.45y_{t-3}^{(1)}-0.85y_{t-3}^{(2)}+\epsilon_{t}^{(4)}\\\
y_{t}^{(5)}=0.95y_{t-2}^{(1)}+\epsilon_{t}^{(5)}\end{cases}$ (47)

The third data set (M3) resembles M1, but is a time series with a higher
order, making M3 a VAR model.

In each of the above cases $\epsilon_{t}^{(1)},\dots,\epsilon_{t}^{(5)}$ are
uncorrelated normally distributed random variables with unit variance.
$y_{1}^{(1)},\dots,y_{1}^{(5)}$ are also equal to zero.

###  4.3 Real data sets

#### Financial equity data set

We consider a financial data set for two forecasting applications, namely
daily returns and volumes. We consider a universe of 210 instruments, which
consists of 10 ETFs (SPY, XLF, XLB, XLK, XLV, XLI, XLU, XLY, XLP, XLE), and
200 equities from the financial sector.

Let $P_{i,d}$ denote the closing price for asset $i$ at day $d$. We consider
the close to close returns888The returns considered are adjusted for splits
and dividends for asset $i$, defined as

$R_{i,d}^{k}=\frac{P_{i,d+k}-P_{i,d}}{P_{i,d}}.$ (48)

As a more representative quantity, grounded within financial applications, we
consider the market excess returns, that is the difference between the returns
of an instrument and that of the S&P-500 index, a widely used proxy for market
returns.

$\tilde{R}_{i,t}^{k}=\tilde{R}_{i,t}^{k}-\tilde{R}_{\textrm{SPY},t}^{k}.$ (49)

Denote the volume of instrument $i$ as $V_{i,t}$, which is defined as the
total number of shares traded during day $t$. Note that this quantity is non-
stationary and generally increases over time, across most instruments. To this
end, we instead consider the ’log-volume-returns’, defined as

$\delta^{k}_{i,t}=\frac{\textrm{log}(V_{i,t+k})-\textrm{log}(V_{i,t})}{\textrm{log}(V_{i,t})}.$
(50)

Predicting daily returns is one of the most popular yet difficult problems in
finance. Due to the low signal-to-noise ratio of this task, standard metrics
like Mean Squared Error and correlation, are considered unsuitable. A more
suitable economic-driven evaluation method is following a hypothetical
transaction in the underlying, based on the forecasted signs (with $+$ and $-$
corresponding to buy and sell respectively). As this is simply for evaluation
purposes, we will ignore transaction costs and liquidity constraints.

Following this hypothetical strategy, we consider the (excess-returns-
adjusted) Profit and Loss (P&L) metric defined by

$\textrm{P}\&\textrm{L}^{k}_{t}=\frac{1}{N}\sum_{i=1}^{N_{t}}\widetilde{R}^{k}_{i,t}\times\textrm{sign}\left(s_{i,t}\right),$
(51)

The P&L is defined above for each day $d$, which yields a vector
$\textbf{P\&L}=\left(\textrm{P\&L}_{1},\textrm{P\&L}_{2},\dots,\textrm{P\&L}_{T-1}\right)$
of daily P&L’s, which measures the outcome of the hypothetical strategy over
the entire time horizon under consideration. A measure of the long-term
performance of a strategy is the annualized Sharpe Ratio (hereafter SR) [43],
defined as

$\textrm{SR}=\frac{\textrm{mean}(\textbf{P\&L})\times\sqrt{252}}{\textrm{sd}(\textbf{P\&L})}.$
(52)

Lastly, we define the P&L per dollar traded (PPD) as

$\textrm{PPD}=\frac{1}{T-1}\sum_{t=1}^{T-1}\textrm{P\&L},$ (53)

which assumes that every day, the total allocated capital to the portfolio
remains constant at 1 USD.

We employ the SR significance test by [44]. As returns are known to exhibit
heavier tails than a normal distribution, this test is an extension of a
traditional $t$-test, adjusted to account for the skewness and kurtosis of the
distribution of the returns.

A higher magnitude forecast should naturally signify a higher confidence. For
each trading day, we consider the quantiles of the forecast magnitudes
$(\\{|s_{i,t}|\\}_{i=1}^{N}$, in order to define five quintile portfolios, of
decreasing size. In particular, we define Q1 to be the portfolio consisting of
all assets, Q2 to be the portfolio consisting of the top 80% assets with the
highest magnitude, and so on, up to Q5 which consists of the 20% highest
magnitude assets.

The feature matrix consists of the following indicators

  * •

1, 2 and 3-day returns for each of the instruments, i.e $R^{1}_{i,t}$,
$R^{2}_{i,t}$, $R^{3}_{i,t}$

  * •

1, 2 and 3-day absolute returns for each of the instruments, i.e
$|R^{1}_{i,t}|$, $|R^{2}_{i,t}|$, $|R^{3}_{i,t}|$

  * •

1, 2 and 3-day log-volume returns for each of the instruments, i.e
$\delta^{1}_{i,t}$, $\delta^{2}_{i,t}$, $\delta^{3}_{i,t}$

  * •

the three Fama & French daily market factors

  * •

the daily returns of 49 different industry portfolios999The industry portfolio
returns and the daily market factors are taken directly from the Kenneth R.
French library [45]. Based on a market portfolio, the three factors include
small minus big (SMB) based on firm size, high minus low (HML) based on the
book-to-market ratio, and the market excess return of the portfolio.

  * •

daily returns of the VIX index

We build an embedding based on all the above features. We set the signal
$s_{i,t}$ equal to the output of the OFTER algorithm for each of the
instruments $i$ and day $t$.

#### Air quality data set

As a non-financial application we also consider the Beijing air quality data
set from the UCI Machine Learning Repository. Data is collected by the Beijing
Municipal Environmental Monitoring Center, and considers various air quality
measurements from 12 monitoring points. We aim at forecasting each of the six
types of pollution (PM2.5,PM10,03,C0,S02,N02) concentrations (in $ug/m^{3}$),
for each of the stations. Other features in the data set, used as predictors,
include Temperature, Dewpoint Temperature, Precipitation, Pressure as well as
Wind speed and Direction.

###  4.4 Numerical results

Data Set | T | d | LM | PCR | LASSO |  | Random
---
Forests
SVM | LSTM | Bi-LSTM | GRU | OFTER | OFTER-DR | OFTER-FT | OFTER-DR-FT
Synthetic | 3000 | 15 | 0.01 | 0,02 | 0.64 | 33.90 | 1.81 | 173.45 | 202.99 | 184.81 | 7.87 | 9.68 | 37.52 | 24.53
Financial | 5022 | 154 | 0.37 | 0.48 | 5.38 | 3,144.02 | 0.17 | 442.86 | 503.78 | 549.30 | 17.58 | 22.25 | 222.87 | 153.06
Air Quality | 1461 | 228 | 0.10 | 0.31 | 2.23 | 132.05 | 1.00 | 165.48 | 152.22 | 184.85 | 3.25 | 20.18 | 157.95 | 62.64

Table 1: The average running time of all 12 methods on all of the forecasted
quantities for each of the three multivariate time series. $T$ denotes the
length of the time series, while $d$ denotes the dimensionality of each data
set.

We start with a discussion of the running time of OFTER. In Table 1, we
provide the average running time of the algorithm for each of the three
considered data sets. We notice that all four configurations of OFTER are
significantly faster than LSTM, Bi-LSTM and GRU, while at the same time
significantly slower than the OLS, LASSO, PCR and SVM, as expected. It is also
clear that among the OFTER configurations, OFTER-FT is the slowest as each
feature is transformed separately without first reducing the dimensionality of
the feature space. The time complexity of OFTER and OFTER-DR are more
comparable, with the latter being slightly more computationally intensive.

For the numerical experiments concerning forecasting returns, we employ
economic-driven metrics typical of this literature, as discussed in Section
4.3. For all the remaining experiments, we base the performance comparisons on
the (Pearson) correlation metric, as it is easily interpretable, in comparison
to mean squared error and other similar metrics.101010Note that the comparison
between the performances of different models is almost identical if we instead
use the mean squared error or the mean absolute error.

####  $\bullet$ Synthetic data

For this task we use $L=800,C_{\textrm{min}}=0.05$ and
$C_{\textrm{Original}}=0.05$. In this experiment, OFTER performs consistently
well. In almost all of the features, one of the OFTER configurations is among
the two best-performing algorithms (with the GRU generally having a very high
performance).

The results on the synthetic data set confirm the expectations from our model.
On M1, GRU appears to outperform the other methods for four out of five
features as shown in Figure 3. Following GRU, OFTER and OFTER-DR appear to
have the highest performance for three out of five features.

In the case of M2, which is a non-linear multivariate non-linear model, the
configurations of OFTER that employ the feature transformation component, are
outperforming all other models, with the exception of Random Forests, in three
out of five features (in particular, the three features which are defined in a
non-linear manner) as can be seen in Figure 4.

Figure 3: Grouped barplot showing a performance comparison of the 12 methods
considered, on each of the five features of the M1 data set described in (45).
The metric of the comparison is the correlation between the target values and
the forecasts.  Figure 4: Grouped barplot showing a performance comparison of
the 12 methods considered, on each of the five features of the M2 data set
described in (46). The metric of the comparison is the correlation between the
target values and the forecasts.  Figure 5: Grouped barplot showing a
performance comparison of the 12 methods considered, on each of the five
features of the M3 model described in (47). The metric of the comparison is
the correlation between the target values and the forecasts.

Lastly, in the case of M3, which regards a vector autoregressive (VAR) model,
OFTER is performing particularly well, outperforming all other models, with
the exception of GRU, as can be seen in Figure 5. Similarly to M1, the two
OFTER configurations that do not employ the feature transformation component,
fare slightly better.

####  $\bullet$ ETF volumes

For both financial tasks, we use $L=800$ and set
$C_{\textrm{min}}=C_{\textrm{Original}}=0.01$, as typically, these are
problems with a lower signal-to-noise ratio.

Figure 6: Groupped barplot showing the correlation between the delta-volume
and its forecasts, by a selection of the baselines and OFTER configurations.
The selection was such that the highest performing methods were selected from
each category of methods.

For the task of predicting the log-delta-volumes of ETF instruments (50),
OFTER outperforms all the other benchmark methods. As demonstrated in Figure
6, both OFTER-FT and OFTER-DR configurations outperform all of the baseline
methods, across all 10 ETFs.

####  $\bullet$ ETF market excess returns

| Q1 | Q2 | Q3 | Q4 | Q5
---|---|---|---|---|---
| SR | PPD | p-value | SR | PPD | p-value | SR | PPD | p-value | SR | PPD | p-value | SR | PPD | p-value
LM | -0.745 | -1.211 | 0.967 | -0.729 | -1.270 | 0.964 | -0.873 | -1.783 | 0.984 | -0.740 | -1.689 | 0.965 | -0.814 | -2.664 | 0.977
PCR | -0.863 | -1.324 | 0.981 | -0.901 | -1.523 | 0.985 | -0.841 | -1.737 | 0.979 | -0.952 | -2.234 | 0.989 | -0.809 | -2.728 | 0.976
LASSO | -0.890 | -1.293 | 0.985 | -0.920 | -1.411 | 0.987 | -1.052 | -1.905 | 0.995 | -1.192 | -2.466 | 0.998 | -0.995 | -3.163 | 0.993
Random Forests | 0.371 | 0.534 | 0.182 | 0.282 | 0.436 | 0.246 | 0.174 | 0.339 | 0.335 | 0.233 | 0.521 | 0.285 | 0.285 | 1.000 | 0.243
SVM | -0.370 | -0.528 | 0.807 | -0.314 | -0.471 | 0.763 | 0.309 | 0.429 | 0.204 | 0.613 | 0.984 | 0.058 | 0.706 | 1.267 | 0.039
LSTM | 0.308 | 0.433 | 0.227 | 0.261 | 0.396 | 0.262 | 0.421 | 0.771 | 0.152 | 0.403 | 0.818 | 0.163 | 0.059 | 0.176 | 0.442
Bi-LSTM | -0.810 | -0.868 | 0.976 | -0.824 | -1.017 | 0.977 | -0.388 | -0.641 | 0.828 | -0.238 | -0.462 | 0.719 | -0.456 | -1.377 | 0.866
GRU | 0.424 | 0.610 | 0.152 | 0.286 | 0.448 | 0.243 | 0.148 | 0.268 | 0.359 | 0.144 | 0.294 | 0.363 | 0.129 | 0.403 | 0.376
OFTER | 0.534 | 0.811 | 0.089 | 0.653 | 1.082 | 0.052 | 0.616 | 1.203 | 0.064 | 0.290 | 0.646 | 0.233 | 0.444 | 1.463 | 0.140
OFTER-DR | 1.616 | 2.346 | 0.000 | 1.704 | 2.607 | 0.000 | 1.229 | 2.194 | 0.001 | 1.258 | 2.441 | 0.001 | 1.261 | 3.765 | 0.001
OFTER-FT | 1.143 | 1.737 | 0.002 | 1.344 | 2.197 | 0.000 | 1.261 | 2.476 | 0.001 | 0.988 | 2.167 | 0.007 | 0.491 | 1.545 | 0.110
OFTER-DR-FT | 1.468 | 2.125 | 0.000 | 1.532 | 2.409 | 0.000 | 1.562 | 2.795 | 0.000 | 0.884 | 1.760 | 0.015 | 0.826 | 2.371 | 0.022

Table 2: In this table we summarize the performance of 12 models for the ETF
price movement forecasting task. SR denotes the annualized Sharpe Ratio, PPD
denotes the Profit Per Dollar, and p-value denotes the p-value obtained from
the one-tailed SR significance test of [44]. Q1-Q5 correspond to the quantile
portfolios as defined in Section 4.3. Figure 7: Cumulative P&L plot for the
daily OFTER-DR generated strategy. SR denotes the annualized Sharpe Ratio, PPD
denotes the Profit Per Dollar, and p-value denotes the p-value obtained from
the one-tailed SR significance test of [44]. Q1-Q5 correspond to the quantile
portfolios as defined in Section 4.3. Figure 8: Cumulative P&L plot of the
daily strategies corresponding to the 12 methods for the ETF price movement
forecasting task. SR denotes the annualized Sharpe Ratio, PPD denotes the
Profit Per Dollar, and p-value denotes the p-value as obtained from the one-
tailed Sharpe Ratio significance test of [44]. The performance metrics for
each method are given in Table 2.

We simulate a daily strategy based on trading ten ETF instruments for the
Q1-Q5 portfolios, and summarize the results in Table 2. For three out of four
OFTER strategies, the SR is considered statistically significant.111111Using
the Sharpe Ratio significance test of [44]. This is based on a generalization
of the t-test, that accounts for the non-normality of returns, by taking into
account the kurtosis and skewness of the P&L distribution. As a visual
demonstration of the results, we also plot the cumulative P&L for the Q1
strategy of each method in Figure 8, clearly demonstrating that OFTER-DR and
OFTER-DR-FT are outperforming the other methods. With the exception of the
Q5-SVM strategy, none of the other benchmark methods result in strategies with
a statistically significant SR. Across most quantiles, OFTER-DR scores the
highest performance, both in terms of SR and PPD. In Figure 7, we plot the
cumulative P&L for the OFTER-DR strategy, clearly demonstrating a consistently
well-performing set of strategies across different quantiles.

####  $\bullet$ Air pollution data set

Figure 9: This heatmap is a visualization of the comparison between the
different forecasting models for the air pollution data set. The proportion of
the number of cases where a method (label of the vertical axis) outperforms
another method (label of the horizontal axis) in terms of correlation is given
in the corresponding cell.

Finally, for the air pollution task, we forecast a total of 6 features, for 12
stations. In Figure 9, we show a grid of the proportions of cases where one
method outperforms another based on the correlation. As it can be seen, all
four OFTER configurations are outperforming all of the other eight benchmarks.
In this particular case, OFTER-DR-FT outperforms the other three versions of
OFTER.

##  5 Applications to feature importance and outlier detection.

In this section, we discuss additional applications of OFTER to two tasks of
interest in the machine learning community. In particular, we consider a
feature importance component that boosts the interpretative power of OFTER,
and an application to outlier detection.

###  5.1 Feature importance

Figure 10:  Figure 11: In this 3x3 plot grid, we show the feature importance
weights, as obtained in Section 5.1. The first feature of each of M1, M2, and
M3 is forecasted in the first, second and third row, respectively. In each
case, there are 15 features corresponding to three lags of each of the five
features. This is repeated for three of OFTER’s configurations, corresponding
to the three columns of the grid.

In the OFTER pipeline, one can put forth a natural approach for estimating the
contribution of each feature to the forecasting task. In particular, we aim to
find a proxy for the importance of each of the input features. To this end,
consider two points in the initial feature space $x_{1}$ and $x_{2}$, and
denote by $x_{1}^{\prime}=U^{T}x_{1}$ and $x_{2}^{\prime}=U^{T}x_{2}$ the
transformed variables. Denote by
$D^{\prime}=||x_{1}^{\prime}-x_{2}^{\prime}||_{v}$ the distance in the reduced
space (used for calculating the nearest neighbors), and let
$D_{j}=x_{1,j}-x_{2,j}$ such that
$||x_{1}-x_{2}||_{2}=\sqrt{\sum_{j}D_{j}^{2}}$.

A metric for the importance of a variable in determining the relative distance
of data points in the embedding space can be modeled as the sensitivity of the
distance $D^{\prime}$ on the individual distance $d_{j}$, i.e
$\left|\frac{\partial D^{\prime}}{\partial D_{j}}\right|.$ For our goal, it is
natural to take the derivative at a point where $D_{1}=\dots=D_{d}=1$. We
begin by writing

$D^{\prime
2}=||x_{1}^{\prime}-x_{2}^{\prime}||_{v}=||U^{T}(x_{1}-x_{2})||_{v}=\sum_{j=1}^{p}v_{j}\left(\sum_{k=1}^{d}U_{k,j}(x_{1}-x_{2})_{k}\right)^{2}=\sum_{j=1}^{p}v_{j}\left(\sum_{k=1}^{d}U_{k,j}D_{k}\right)^{2}.$
(54)

We then compute

$ \frac{\partial D^{\prime 2}}{\partial D_{k}}\Bigr{|}_{d=e_{d}} =
2\sum_{j=1}^{p}v_{j}\left(\sum_{k^{\prime}=1}^{d}U_{k^{\prime},j}D_{k^{\prime}}
\right)U_{k,j},$ (55)

which, using $\odot$ to denote the Hadamard product, can be written as

$ \nabla D^{\prime 2}(e_{d})=2\left|U \left[(U^{T}e_{d})\odot v\right]\right|.
$ (56)

Keeping in mind (41), we also augment the embedding space with the original
features that are highly correlated with the target. Consider the vector
$v_{2}$ containing the weight given to the original feature if used, and zero
otherwise, i.e

$v_{2,j}=\begin{cases}v_{I^{\textrm{Original}}_{k}}&\text{$j=I^{\textrm{Original}}_{k}$}\\\
0&\text{$j\notin I^{\textrm{Original}}$}\end{cases}.$

This now leads to

$D^{\prime
2}=\sum_{j=1}^{p}v_{j}\left(\sum_{k=1}^{d}U_{k,j}D_{k}\right)^{2}+\sum_{l=1}^{d}v_{2,l}D_{l}^{2},$
(57)

and thus the derivative is instead given by

$ \nabla D^{\prime 2}(e_{d})=2\left|U \left[(U^{T}e_{d})\odot
v\right]+v_{2}\right|. $ (58)

As an illustration, we consider the task of forecasting the first feature of
each of the three synthetic data sets (M1-M3), using OFTER-DR, OFTER-FT, and
OFTER-DR-FT. In Figure 11, we show a grid of the feature importance of each
individual feature. As we ignore the unimportant features via the calculation
of $v$ in Equation 43, OFTER-FT attaches a non-zero importance weight only to
a subset of the features. On the other hand, OFTER-DR assigns a non-zero
weight to all of the original features, as they all contribute to the
estimation of each principal component. This also holds true for OFTER-DR-FT,
where dimensionality reduction is also performed. However, in the case of M2
where most of the predictability is present in a single feature (M2.1 - Lag1),
OFTER-DR-FT assigns a zero importance weight to most of the features,
similarly to OFTER-FT. Therefore, for the task of feature importance, we
recommend the user to rely on either OFTER–FT or OFTER-DR-FT.

###  5.2 Outlier detection

Figure 12: Plot showing the minimum distance $d_{j}^{\textrm{min}}$ obtained
using the suggested outlier detection algorithm for OFTER-DR, $L=600$ and
$\kappa=5$. All numbered points are identified as outliers.
1\. 18 September 2001 - Aftermath of the 9/11 Attack
2\. 24 July 2002 - Stock market downturn of 2002 ($12^{\textrm{th}}$
consecutive day of falling prices)
3\. 29 July 2002 - Stock market downturn of 2002 ($3^{\textrm{rd}}$ day of
rising prices after sharp price rises following (2))
4\. 28 February 2007 - Big drop in US and China stocks
5\. 10 October 2008 - Market Crash across Europe and Asia
6\. 7 May 2010 - Following the Flash Crash of May 6
7\. 4 August 2011 - Preceding the Black Monday (8 August)
8\. 15 April 2013 - One of the biggest drops in 2013 for DJIA, SPX
9\. 8 November 2016 - US election

A byproduct of the OFTER pipeline is that it can facilitate the task of
outlier detection, by considering the pairwise distance between time points in
the embedding space. More specifically, points whose nearest neighbors are
still quite “far” in the embedding space could be labeled as outliers, in the
same spirit as the one-class SVM algorithm [46]. We demonstrate this using the
following example.

Define the minimum distance between the current time point $t$ and the time
points in its recent history in the reduced space as

$d_{t}^{\textrm{min}}=\textrm{min}(\\{d_{t-L,t},\dots,d_{t-1,t}\\}).$ (59)

Comparing $d_{t}^{\textrm{min}}$ with the distribution of
$(d_{t-L}^{\textrm{min}},\dots,d_{t-1}^{\textrm{min}})$ may constitute input
for any outlier detection algorithm.121212Instead of $d_{t}^{\textrm{min}}$,
one could employ the average distance
$\bar{d}_{t}=\frac{\sum_{j=t-L}^{t-1}d_{j,t}}{t}$. Denoting by
$q^{\textrm{min}}_{t,1}$ and $q^{\textrm{min}}_{t,3}$ the first and third
quantiles of $\\{d^{\textrm{min}}_{t-L},\dots,d^{\textrm{min}}_{t-1}\\}$, we
classify a time point $t$ as an outlier if
$d_{j}^{\textrm{min}}>q^{\textrm{min}}_{t,3}+\kappa\left(q^{\textrm{min}}_{t,3}-q^{\textrm{min}}_{t,1}\right)$
for some $\kappa>0$.

In Figure 12, we plot the minimum distance between each trading day and its
$L=600$ preceding days, by using $\kappa=5$. From this figure, one could infer
several different regimes, as there are periods where the levels of
$d_{\textrm{min}}$ are low, and periods where these are high, such as the
periods around the 2008 financial crisis. The fact that almost all labeled
outliers fall in the immediate vicinity of days significant in the financial
history, does indeed indicate that the algorithm can successfully be used to
uncover abnormality signals.

Finally, we remark that selecting a larger value of $L$ will result in a large
number of anomalies being highlighted. In practice, if this outlier detection
approach is used in a trading application, a smaller $\kappa$ would be
desirable, as the outlier warning could possibly be qualitatively analyzed.
Note that the above quantile-based outlier detection approach is just a sample
algorithm to illustrate the potential of OFTER, as the pipeline is amenable to
be incorporated into any other outlier detection algorithm.

##  6 Conclusion and future work

In this paper, we have introduced OFTER, a complete pipeline for online
multivariate time series forecasting. By leveraging kNN and GRNN, both simple
well-understood methods with trivial training complexities, OFTER manages to
outperform a number of baselines, for mid-sized data sets. In particular, an
extensive set of numerical experiments demonstrated that OFTER performed
competitively with Random Forests and GRU in the synthetic data set, and
outperformed all other methods in the real-world applications considered, with
a running time much lower than those of deep learning methods. The
shortcomings due to the curse of dimensionality are overcome via
dimensionality reduction and the use of a weighted norm when comparing
distances in the embedding space.

Furthermore, OFTER brings additional advantages, including its flexible
nature, with the ability to use or omit many of its components. As a
byproduct, OFTER can be integrated seamlessly into an outlier detection
pipeline, in an online manner. Lastly, OFTER is naturally equipped to return
the user a feature importance score for all the input variables.

Future work. Future extensions of this work could include the utilization of
other machine learning methods in the model selection step of the algorithm,
for example, LASSO models with different regularization parameters. One
disadvantage of OFTER is that it struggles to discover non-linear interaction
effects between covariates, even though this is a general issue with models
and data sets of this size. One potential extension would be to include a pre-
screening component, where a number of potential interaction terms are tested
in order to augment the feature space, either prior to or after the embedding
step. Additionally, it would be interesting to see the performance of OFTER in
other mid-sized data sets and assess whether other dimensionality reduction
algorithms are able to attain improved performance, relative to the use of PCA
(this was not the case for the data sets used in this project). Yet another
interesting direction concerns an adaptation of OFTER for probabilistic
forecasting. Probabilistic models have recently gained in popularity, being
able to provide not only forecasts but also associated uncertainty estimates,
which is of crucial importance in certain applications. Lastly, we note that
the OSMC metric we introduced is of independent interest, both as an
alternative to the maximal correlation but also as a computationally fast
additional step for other simple regression methods.

## Acknowledgments

This work was supported by the UK Engineering and Physical Sciences Research
Council (EPSRC) grants EP/R513295/1 and EP/V520202/1.

## References

  * [1] Ahmed NK, Atiya AF, Gayar NE, El-Shishiny H.  An Empirical Comparison of Machine Learning Models for Time Series Forecasting.  Econometric Reviews. 2010;29(5-6):594-621.  Available from: https://doi.org/10.1080/07474938.2010.481556.
  * [2] Hochreiter S, Schmidhuber J.  Long Short-Term Memory.  Neural Comput. 1997 nov;9(8):1735–1780.  Available from: https://doi.org/10.1162/neco.1997.9.8.1735.
  * [3] Cornegruta S, Bakewell R, Withey S, Montana G.  Modelling Radiological Language with Bidirectional Long Short-Term Memory Networks.  In: Proceedings of the Seventh International Workshop on Health Text Mining and Information Analysis. Auxtin, TX: Association for Computational Linguistics; 2016. p. 17-27.  Available from: https://aclanthology.org/W16-6103.
  * [4] Chung J, Gulcehre C, Cho K, Bengio Y.  Empirical evaluation of gated recurrent neural networks on sequence modeling.  In: NIPS 2014 Workshop on Deep Learning, December 2014; 2014. .
  * [5] Flunkert V, Salinas D, Gasthaus J.  DeepAR: Probabilistic Forecasting with Autoregressive Recurrent Networks.  International Journal of Forecasting. 2017 04;36.
  * [6] Lim B, Arık SO, Loeff N, Pfister T.  Temporal Fusion Transformers for interpretable multi-horizon time series forecasting.  International Journal of Forecasting. 2021;37(4):1748-64.  Available from: https://www.sciencedirect.com/science/article/pii/S0169207021000637.
  * [7] de O Santos Júnior DS, de Oliveira JFL, de Mattos Neto PSG.  An intelligent hybridization of ARIMA with machine learning models for time series forecasting.  Knowledge-Based Systems. 2019;175:72-86.  Available from: https://www.sciencedirect.com/science/article/pii/S0950705119301327.
  * [8] Firmino PRA, de Mattos Neto PSG, Ferreira TAE.  Correcting and combining time series forecasters.  Neural Networks. 2014;50:1-11.  Available from: https://www.sciencedirect.com/science/article/pii/S0893608013002517.
  * [9] Panigrahi S, Behera HS.  A hybrid ETS–ANN model for time series forecasting.  Engineering Applications of Artificial Intelligence. 2017;66:49-59.  Available from: https://www.sciencedirect.com/science/article/pii/S0952197617301550.
  * [10] Pai PF, Lin CS.  A hybrid ARIMA and support vector machines model in stock price forecasting.  Omega. 2005;33(6):497-505.  Available from: https://www.sciencedirect.com/science/article/pii/S0305048304001082.
  * [11] Zhang GP.  Time series forecasting using a hybrid ARIMA and neural network model.  Neurocomputing. 2003;50:159-75.  Available from: https://www.sciencedirect.com/science/article/pii/S0925231201007020.
  * [12] Zyga K, Price R, Williams B.  A generalized regression neural network for logo recognition.  In: KES’2000. Fourth International Conference on Knowledge-Based Intelligent Engineering Systems and Allied Technologies. Proceedings (Cat. No.00TH8516). vol. 2; 2000. p. 475-8 vol.2.
  * [13] Amanda JC, Gary AR, Tami MBB, Joseph PC, Brittney NK.  Feed-forward and generalised regression neural networks in modelling feeding behaviour of pigs in the grow-finish phase.  Biosystems Engineering. 2018;173:124-33.  Advances in the Engineering of Sensor-based Monitoring and Management Systems for Precision Livestock Farming.  Available from: https://www.sciencedirect.com/science/article/pii/S1537511017304750.
  * [14] Pal M.  Modelling pile capacity using generalised regression neural network.  In: Proceedings of Indian geotechnical conference December; 2011. p. 15-7.
  * [15] Al-Mahasneh AJ, Anavatti S, Garratt M, Pratama M.  Applications of General Regression Neural Networks in Dynamic Systems.  In: Asadpour V, editor. Digital Systems. Rijeka: IntechOpen; 2018. p. 133-54.  Available from: https://doi.org/10.5772/intechopen.80258.
  * [16] J Wang JZ, M Xu HW.  Classification of Imbalanced Data by Using the SMOTE Algorithm and Locally Linear Embedding.  Signal Processing, 2006 8th International Conference on Vol 3 IEEE. 2007\.
  * [17] Tenenbaum JB, Silva V, Langford JC.  A Global Geometric Framework for Nonlinear Dimensionality Reduction.  Science. 2000;290(5500):2319-23.
  * [18] Donoho DL, Grimes C.  Hessian Eigenmaps: Locally Linear Embedding Techniques for High-Dimensional Data.  Proceedings of the National Academy of Sciences. 2003;100(10):5591-6.
  * [19] Belkin M, Niyogi P.  Laplacian Eigenmaps for Dimensionality Reduction and Data Representation.  Neural Comput. 2003 Jun;15(6):1373-96.  Available from: http://dx.doi.org/10.1162/089976603321780317.
  * [20] Coifman RR, S L.  Diffusion maps.  Applied and Computational Harmonic Analysis. 2006;21(1):5 30.  Special Issue: Diffusion Maps and Wavelets.  Available from: http://www.sciencedirect.com/science/article/pii/S1063520306000546.
  * [21] Kramer MA.  Nonlinear principal component analysis using autoassociative neural networks.  AIChE Journal. 1991;37(2):233-43.  Available from: https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.690370209.
  * [22] Pulgar FJ, Charte F, Rivera AJ, del Jesus MJ.  AEkNN: An AutoEncoder kNN-Based Classifier With Built-in Dimensionality Reduction.  International Journal of Computational Intelligence Systems. 2019;12:436-52.  Available from: https://doi.org/10.2991/ijcis.2018.125905686.
  * [23] Huang D, Jiang F, Li K, Tong G, Zhou G.  Scaled PCA: A New Approach to Dimension Reduction.  Management Science. 2022;68(3):1678-95.  Available from: https://doi.org/10.1287/mnsc.2021.4020.
  * [24] Chiavazzo E, Gear C, Dsilva C, Rabin N, Kevrekidis I.  Reduced Models in Chemical Kinetics via Nonlinear Data-Mining.  Processes. 2013 07;2.
  * [25] Papaioannou PG, Talmon R, Kevrekidis IG, Siettos C.  Time-series forecasting using manifold learning, radial basis function interpolation, and geometric harmonics.  Chaos: An Interdisciplinary Journal of Nonlinear Science. 2022;32(8):083113.  Available from: https://doi.org/10.1063/5.0094887.
  * [26] Yan X, Bai J, Li X, Chen Z.  Can dimensional reduction technology make better use of the information of uncertainty indices when predicting volatility of Chinese crude oil futures?  Resources Policy. 2022;75:102521.  Available from: https://www.sciencedirect.com/science/article/pii/S0301420721005286.
  * [27] Joseph K, Babajide Wintoki M, Zhang Z.  Forecasting abnormal stock returns and trading volume using investor sentiment: Evidence from online search.  International Journal of Forecasting. 2011;27(4):1116-27.  Available from: https://www.sciencedirect.com/science/article/pii/S0169207011000021.
  * [28] Liu X, keung Lai K.  Intraday volume percentages forecasting using a dynamic SVM-based approach.  Journal of Systems Science and Complexity. 2017;30:421-33.
  * [29] Libman D, Haber S, Schaps M.  Volume Prediction With Neural Networks.  Frontiers in Artificial Intelligence. 2019;2.  Available from: https://www.frontiersin.org/article/10.3389/frai.2019.00021.
  * [30] Fix E, Hodges JL.  Discriminatory Analysis. Nonparametric Discrimination: Consistency Properties.  International Statistical Review / Revue Internationale de Statistique. 1989;57(3):238-47.  Available from: http://www.jstor.org/stable/1403797.
  * [31] Specht DF.  A general regression neural network.  IEEE Transactions on Neural Networks. 1991;2(6):568-76.
  * [32] Bouchaud JP, Potters M.  Financial applications of random matrix theory: a short review.  In: The Oxford Handbook of Random Matrix Theory. Oxford University Press; 2015. Available from: https://doi.org/10.1093/oxfordhb/9780198744191.013.40.
  * [33] Breiman L, Friedman JH.  Estimating Optimal Transformations for Multiple Regression and Correlation.  Journal of the American Statistical Association. 1985;80(391):580-98.  Available from: https://www.tandfonline.com/doi/abs/10.1080/01621459.1985.10478157.
  * [34] Lorentz GG.  Bernstein Polynomials.  AMS Chelsea Publishing Series. Chelsea Publishing Company; 1986.  Available from: https://books.google.co.uk/books?id=nXiRAKzpdjUC.
  * [35] Hosaka M.  In: Bézier Curves and Control Points. Berlin, Heidelberg: Springer Berlin Heidelberg; 1992. p. 117-39.  Available from: https://doi.org/10.1007/978-3-642-76598-8_9.
  * [36] Piegl L, Tiller W.  In: B-spline Curves and Surfaces. Berlin, Heidelberg: Springer Berlin Heidelberg; 1997. p. 81-116.  Available from: https://doi.org/10.1007/978-3-642-59223-2_3.
  * [37] Mitz R, Shkolnisky Y.  ROIPCA: An Online PCA algorithm based on rank-one updates.  CoRR. 2019;abs/1911.11049.  Available from: http://arxiv.org/abs/1911.11049.
  * [38] Bunch JR, Nielsen CP, Sorensen DC.  Rank-one modification of the symmetric eigenproblem.  Physica Scripta. 1978.
  * [39] Mitz R, Sharon N, Shkolnisky Y.  Symmetric Rank-One Updates from Partial Spectrum with an Application to Out-of-Sample Extension.  SIAM J Matrix Anal Appl. 2017;40:973-97.
  * [40] Melman A.  Numerical solution of a secular equation.  Numerische Mathematik. 1995.
  * [41] Bergmeir C, Benítez JM.  On the use of cross-validation for time series predictor evaluation.  Information Sciences. 2012;191:192-213.  Data Mining for Software Trustworthiness.  Available from: https://www.sciencedirect.com/science/article/pii/S0020025511006773.
  * [42] Nicolaou N, Constandinou T.  A nonlinear causality estimator based on non-parametric multiplicative regression.  Frontiers in Neuroinformatics. 2016 Jun;10(JUNE).
  * [43] Sharpe WF.  The Sharpe Ratio.  The Journal of Portfolio Management. 1994;21(1):49-58.
  * [44] Bailey D, Lopez de Prado M.  The Deflated Sharpe Ratio: Correcting for Selection Bias, Backtest Overfitting, and Non-Normality.  The Journal of Portfolio Management. 2014 09;40:94-107.
  * [45] Kenneth R. French - Data Library;.  Accessed: 2023-02-02.  http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html.
  * [46] Schölkopf B, Williamson R, Smola A, Shawe-Taylor J, Platt J.  Support Vector Method for Novelty Detection.  In: Proceedings of the 12th International Conference on Neural Information Processing Systems. NIPS’99. Cambridge, MA, USA: MIT Press; 1999. p. 582–588.
  * [47] Godfrey MD.  Elements of Numerical Analysis, by Peter Henrici (John Wiley and Sons, New York, 1964), pp. 328.  Naval Research Logistics Quarterly. 1965 June;12(2):195-5.  Available from: https://ideas.repec.org/a/wly/navlog/v12y1965i2p195-195.html.

##  Appendix A A convexifying approach for solving the Secular equation

In this section, we give a brief description of the method we use to solve the
secular equation described in Section 3.1, which amounts to finding the roots
of

$f(\kappa)=1+\rho\sum_{j=1}^{d}\frac{z_{j}^{2}}{\lambda_{j}-\kappa}.$ (60)

The main difficulty in directly solving Equation 60 lies in the multitude of
discontinuities in the domain of $f$. An efficient approach to find the
$i^{\textrm{th}}$ root of the equation, would be to consider the
reparameterizations $\kappa=\lambda_{i}+\rho\xi$ and
$\delta_{j}=\frac{\lambda_{j}-\lambda_{i}}{\rho}$, which leads to the
formulation of Equation 60 as

$f_{i}(\xi)=1+\sum_{j=1}^{d}\frac{z_{j}^{2}}{\delta_{j}-\xi}.$ (61)

By assuming that each eigenvalue is distinct (a reasonable assumption when
dealing with real data), we notice that

$\delta_{1}<\delta_{2}<\dots<\delta_{i}=0<\dots<\delta_{n}.$ (62)

Therefore, we look for solutions in the set $\left(0,\delta_{i+1}\right)$, for
$i<n$, and in $\mathbb{R}_{+}$ for $i=n$. We can further transform $f_{i}(x)$
into a convex function, by considering $F_{i}(\gamma)=f_{i}(\gamma^{-1/2})$,
which can be shown to transform $F_{i}$ into a convex function in the interval
$(\frac{1}{\delta_{i+1}},\infty)$ for $i<n$ and $(0,\infty)$ for $i=n$.

Finally, by using the fact that, for a decreasing convex function, the Newton
method131313We remind the reader that the Newton method searches for a
solution to $f(x)=0$, starting from point $x_{0}$, and then considers the
sequence  $x_{n+1}=x_{n}-\frac{f(x_{n})}{f^{\prime}(x_{n})},$ until a solution
is obtained. monotonically converges to a solution [47], it can be established
that, if $\gamma^{*}$ is the solution to $F_{i}(\gamma)=0$, then the Newton
method converges monotonically to $\gamma^{*}$ with any starting point in
$\left[\frac{1}{\delta_{i+1}^{2}},\gamma^{*}\right]$. For a more detailed
explanation and an alternative based on the secant method, we refer the reader
to [40].

##  Appendix B Brief summary of baseline methods

The following methods are used as baselines against OFTER. We provide a short
description of each one.

### Ordinary Least Squares (OLS)

Ordinary Least Squares (OLS) is by far the most well known regression method,
where the target variable $y$ is modeled as a linear combination of the
features

$y_{i}=\beta_{0}+\sum_{j=1}^{d}\beta_{j}X_{i,j}+\epsilon_{i}.$ (63)

The coefficient vector $\beta=(\beta_{0},\dots,\beta_{d})$, is chosen as the
one maximizing $\sum\epsilon_{i}^{2}$, and can be found using the closed form
formula

$\beta=(X^{T}X)^{-1}X^{T}y.$ (64)

Since we include lagged features in the feature matrix, this method is very
similar to a Vector Auto-regressive Regression (VAR) model.

### Principal Component Regression

To combat the high correlation present in the features, Principal Component
Regression (PCR) first identifies the principal components of $X$, obtaining
$\tilde{X}$, and then a number of components $p<d$ is retained. OLS is then
used, but on the first $p$ columns of $\tilde{X}$ instead of $X$. Similar to
OFTER, we use $\delta=0.9$ for choosing the number of components.

### Least Absolute Shrinkage and Selection Operator

An alternative to reducing the dimensionality of $X$ before OLS, is to include
a regularization term in the objective function Equation 63. Least Absolute
Shrinkage and Selection Operator (LASSO) is a widely popular method, where the
modified objective function takes the following form

$(y-X\beta)^{T}(y-X\beta)+\lambda\sum_{j}|\beta_{j}|.$ (65)

Note that in this case, there is no closed-form solution, and thus a numerical
optimization method is utilized. We choose $\lambda$ using a cross-validation
procedure.

### Random Forests

A decision tree is a simple method inspired by decision analysis, where the
visualization of the model resembles a tree consisting of nodes and branches.
The initial feature space constitutes the root node. The observations in the
root node are partitioned into ’successor nodes’, based on a splitting rule
chosen so as to optimize some criterion (usually either entropy or Gini
Impurity). The process is repeated recursively on each of the successor nodes
until some stopping criterion is met. For each new data point, its features
are then used to determine its corresponding leaf node. Forecasts can then be
obtained by averaging the target value of all observations that fall within
each leaf node. This is a non-parametric method, but with a very high
variance.

A Random Forest is built by fitting multiple decision trees. Each decision
tree is fitted using only a random sample of observations and a random sample
of features. The forecasts of the decision trees are then combined into a
single forecast (usually by averaging). Random Forests are based on the idea
of bagging aimed at reducing variance and avoiding overfitting. For the
benchmark task, we use 200 trees.

### Support Vector Machines (SVM):

SVM is another popular machine learning model that is motivated by learning
theory. In particular, a vector $\beta$ and a value $\beta_{0}$ are found that
minimize

$\frac{1}{2}||\beta||^{2},\ \ \ \textrm{s.t}\ \ \
|(y_{i}-\beta^{T}x_{i}-\beta_{0})|\leq\epsilon.$ (66)

A prediction for a new point $x^{\prime}$ is then given by
$\beta_{0}+\beta^{T}x^{\prime}$.

### Recurrent Neural Networks:

Recurrent neural networks (RNN) form a class of neural networks specifically
aimed at data sets with temporal structures. This is generally achieved by
introducing additional skip connections that pass the output of a past layers
to the current one. One issue caused is that gradients have the tendency to
vanish or explode.

Long-Short Term Memory is a variant of Recurrent Neural Networks, specifically
designed to avoid the issue of exploding or vanishing gradients. This is
achieved by introducing a cell state, which contains a summary of the previous
cells and is controlled by a number of different gates. Denote by $x_{t}$ the
feature vector of a data point, $c_{t}$ the input of the cell state, and
$h_{t}$ the output of the module. The formula of the basic LSTM takes the form

$\begin{cases}f_{t}=\sigma\left(W_{f}x_{t}+U_{f}h_{t-1}+b_{f}\right),\\\
i_{t}=\sigma\left(W_{i}x_{t}+U_{i}h_{t-1}+b_{i}\right),\\\
o_{t}=\sigma\left(W_{o}x_{t}+U_{o}h_{t-1}+b_{o}\right),\\\
\hat{c}_{t}=\textrm{tanh}\left(W_{c}x_{t}+U_{c}h_{t-1}+b_{c}\right),\\\
c_{t}=f_{t}\odot
c_{t-1}+i_{t}\odot\hat{c}_{t}.h_{t}=c_{t}\odot\textrm{tanh}(c_{t}),\\\
\end{cases}$ (67)

where $W_{f},U_{f},b_{f},W_{i},U_{i},b_{i},W_{o},U_{o},b_{o},W_{c},U_{c}$ and
$b_{c}$ are trainable parameters.

We consider two more variants of LSTM. Bi-LSTM which utilizes two parallel and
reverse networks of LSTMs so that information can flow in both directions. We
also use Gated Recursion Unit (GRU) which has a similar yet simpler
architecture than LSTM.

For each of these two layers, we use a total of two layers and 300 epochs.

##  Appendix C Robustness analysis

In this short section, we demonstrate how changing five of the main parameters
used in OFTER-DR-FT affects its performance on the synthetic data set. We fix
five of the main parameters given in OFTER as follows: $\delta=0.9$, $L=800$,
$C_{\textrm{min}}=0.05$, $C_{\textrm{Original}}=0.05$, and
$p_{\textrm{ADF}}=0.05$. By keeping the other four parameters fixed, we vary
each of these parameters in turn and note the variation in its performance. In
Table 3, we show the correlation between the target and the forecasts when
changing each of the noted parameters to the given value.

In general, the performance of the algorithm does not significantly depend on
the choice of parameters, as forecasting each feature has consistent
performance. However there are exceptions; for example, $C_{\textrm{min}}$
appears to have a significant effect on the performance of the algorithm for
M2.3 and M2.4.

| $\delta$ | $p_{\textrm{ADF}}$ | $L$ | $C_{\textrm{min}}$ | $C_{\textrm{original}}$
---|---|---|---|---|---
Value | 0.8 | 0.95 | 0.01 | 0.1 | 500 | 1500 | 0.01 | 0.1 | 0.01 | 0.1
M1.1 | 0.470 | 0.459 | 0.467 | 0.467 | 0.441 | 0.477 | 0.464 | 0.464 | 0.465 | 0.461
M1.2 | 0.519 | 0.511 | 0.509 | 0.509 | 0.503 | 0.521 | 0.509 | 0.514 | 0.504 | 0.518
M1.3 | 0.077 | 0.071 | 0.085 | 0.085 | 0.085 | 0.113 | 0.119 | 0.116 | 0.085 | 0.117
M1.4 | 0.554 | 0.558 | 0.552 | 0.552 | 0.555 | 0.556 | 0.553 | 0.556 | 0.545 | 0.559
M1.5 | 0.094 | 0.107 | 0.101 | 0.101 | 0.147 | 0.120 | 0.102 | 0.110 | 0.101 | 0.113
M2.1 | 0.599 | 0.623 | 0.599 | 0.599 | 0.596 | 0.603 | 0.675 | 0.600 | 0.599 | 0.60
M2.2 | 0.535 | 0.535 | 0.535 | 0.535 | 0.534 | 0.539 | 0.897 | 0.535 | 0.536 | 0.536
M2.3 | 0.864 | 0.799 | 0.815 | 0.815 | 0.776 | 0.85 | 0.759 | 0.401 | 0.804 | 0.817
M2.4 | 0.962 | 0.964 | 0.963 | 0.963 | 0.963 | 0.963 | 0.962 | 0.963 | 0.963 | 0.963
M2.5 | 0.679 | 0.681 | 0.681 | 0.681 | 0.678 | 0.682 | 0.681 | 0.678 | 0.680 | 0.679
M3.1 | 0.758 | 0.758 | 0.758 | 0.758 | 0.751 | 0.762 | 0.758 | 0.758 | 0.758 | 0.757
M3.2 | 0.852 | 0.851 | 0.851 | 0.851 | 0.85 | 0.852 | 0.851 | 0.850 | 0.849 | 0.850
M3.3 | 0.669 | 0.664 | 0.665 | 0.665 | 0.662 | 0.668 | 0.665 | 0.666 | 0.664 | 0.667
M3.4 | 0.902 | 0.903 | 0.903 | 0.903 | 0.902 | 0.905 | 0.903 | 0.902 | 0.903 | 0.902
M3.5 | 0.826 | 0.827 | 0.827 | 0.827 | 0.825 | 0.828 | 0.827 | 0.826 | 0.825 | 0.825
Table 3: Robustness analysis of OFTER. The table displays the correlation
between the target and forecasts for each of the 15 features of the synthetic
data set, by varying a subset of five parameters of the algorithm.