Forecasting and Granger Modelling with Non-linear Dynamical Dependencies

TL;DR

This paper introduces a novel kernel-based method for forecasting multivariate time series that captures non-linear dependencies and uncovers hidden dynamic relationships, outperforming traditional linear models.

Contribution

It develops a new approach for learning vector-valued functions with multiple matrix-valued kernels directly from data, enabling better modeling of non-linearities and dynamic dependencies.

Findings

Superior predictive performance on non-linear series

Effective recovery of hidden dynamic relationships

Alternative to existing graphical Granger methods

Abstract

Traditional linear methods for forecasting multivariate time series are not able to satisfactorily model the non-linear dependencies that may exist in non-Gaussian series. We build on the theory of learning vector-valued functions in the reproducing kernel Hilbert space and develop a method for learning prediction functions that accommodate such non-linearities. The method not only learns the predictive function but also the matrix-valued kernel underlying the function search space directly from the data. Our approach is based on learning multiple matrix-valued kernels, each of those composed of a set of input kernels and a set of output kernels learned in the cone of positive semi-definite matrices. In addition to superior predictive performance in the presence of strong non-linearities, our method also recovers the hidden dynamic relationships between the series and thus is a new alternative to existing graphical Granger techniques.

Tables

Table 1. Table 1 : Synthetic experiments: MSE (std) for 1-step ahead forecasts (hold-out sample average)

Train size	300		700		1000
Mean	0.925	(0.047)	0.923	(0.047)	0.923	(0.047)
LAR	0.890	(0.045)	0.890	(0.044)	0.890	(0.044)
LVAR	0.894	(0.045)	0.836	(0.041)	0.763	(0.035)
LVARL1	0.787	(0.037)	0.737	(0.031)	0.722	(0.030)
NVAR	0.835	(0.041)	0.735	(0.032)	0.719	(0.030)
\hdashlineNVARL1	0.754	(0.034)	0.706	(0.030)	0.679	(0.028)
NVARL12	0.808	(0.040)	0.710	(0.031)	0.684	(0.029)
Train size	1500		2000		3000
Mean	0.923	(0.047)	0.922	(0.047)	0.922	(0.047)
LAR	0.888	(0.045)	0.889	(0.045)	0.888	(0.045)
LVAR	0.751	(0.034)	0.741	(0.033)	0.687	(0.028)
LVARL1	0.710	(0.029)	0.701	(0.028)	0.693	(0.028)
NVAR	0.699	(0.028)	0.682	(0.027)	0.662	(0.026)
\hdashlineNVARL1	0.654	(0.026)	0.640	(0.025)	0.626	(0.025)
NVARL12	0.659	(0.027)	0.685	(0.028)	0.657	(0.027)

Equations

$$\hat{\mathbf{f}}=\underset{\mathbf{f}\in \mathcal{F}}{\mathrm{argmin}}R(\mathbf{f})$$

$$⟨\mathbf{y},\mathbf{g}(\mathbf{x})⟩=⟨{\mathbf{H}}_{\mathbf{x}}\mathbf{y},\mathbf{g}{⟩}_{\mathcal{F}}\forall (\mathbf{y},\mathbf{x},\mathbf{g})\in (\mathcal{Y},\mathcal{X},\mathcal{F}),$$

$$\hat{\mathbf{f}}=\sum _{t=1}^T{\mathbf{H}}_{{\mathbf{x}}_t}{\mathbf{c}}_t,$$

$$\sum_{t=1}^{T}\big{(}\mathbf{H}(\mathbf{x}_{s},\mathbf{x}_{t})+\lambda\delta_{st}\big{)}\,\mathbf{c}_{t}=\mathbf{y}_{s},\quad\forall s\in\mathbb{N}_{n}\enspace,$$

$${\mathcal{X}}^{(j)}=\{{\mathbf{x}}^{(j)}:x_i^{(j)}=x_{i^{(j)}}\forall i,\mathbf{x}\in \mathcal{X}\},{\cup }_j{\mathcal{X}}^{(j)}=\mathcal{X}$$

$$\mathbf{H}=\sum _{j=1}^m{\mathbf{H}}^{(j)}=\sum _{j=1}^m\sum _{i=1}^{s_j}{\mathbf{H}}^{(ji)}=\sum _{j=1}^m\sum _{i=1}^{s_j}{k}^{(ji)}{\mathbf{L}}^{(ji)}$$

$$J(\mathbf{f},\mathcal{L}):=R(\mathbf{f})+\tau Q(\mathcal{L}),\tau >0,$$

$$\tilde{R}(\mathbf{C}):=||\mathbf{Y}-\sum _{ji}{\mathbf{K}}^{(ji)}\mathbf{C}{\mathbf{L}}^{(ji)}|{|}_F^2+\lambda \sum _{ji}⟨{\mathbf{C}}^T{\mathbf{K}}^{(ji)}\mathbf{C},{\mathbf{L}}^{(ji)}{⟩}_F.$$

$$\tilde{J}(\mathbf{C},\mathcal{L}):=\tilde{R}(\mathbf{C},\mathcal{L})+\tau Q(\mathcal{L}).$$

$$Q_1(\mathcal{L})=\sum _{ji}||{\mathbf{L}}^{(ji)}|{|}_1=\sum _{ji}\sum _s^m|L_{ss}^{(ji)}|.$$

$$\mathit{Q}_{1/2}(\mathcal{L})=\sum_{js}\sqrt{\sum_{i}\big{(}L^{(ji)}_{ss}\big{)}^{2}}\enspace.$$

$$\mathbf{A}=\big{(}diag(\mathbf{L}^{11}),\ldots,diag(\mathbf{L}^{ms_{m}})\big{)}^{T}$$

$$\underset{\mathbf{A},\mathbf{C}}{\mathrm{argmin}}\ddot{J}(\mathbf{C},\mathbf{A})$$

$$\ddot{R}(\mathbf{C},\mathbf{A})$$

$${\ddot{Q}}_1(\mathbf{A})=||\mathbf{A}|{|}_1=\sum _{rs}|A_{rs}|$$

$$\ddot{Q}_{1/2}(\mathbf{A})=\sum_{js}\sqrt{\sum_{i}\big{(}A_{(ji)s}\big{)}^{2}}$$

$$P(\mathbf{c},\mathbf{a}):=||\mathbf{y}-\sum _d{a}_d{\mathbf{K}}^d\mathbf{c}|{|}_2^2+\lambda \sum _{ji}{a}_d{\mathbf{c}}^T{\mathbf{K}}^d\mathbf{c}+\sum _d{a}_d.$$

$$\tilde{P}(\mathcal{Z},\mathbf{a}):=||\mathbf{y}-\sum _d{\mathbf{\Phi }}^d{\mathbf{z}}^d|{|}_2^2+\sum _d\left(\frac{\lambda ||{\mathbf{z}}^d|{|}_2^2}{a_d}+a_d\right).$$

$$\underset{\mathbf{a}}{\min }\tilde{P}(\mathcal{Z},\mathbf{a})=||\mathbf{y}-\sum _d{\mathbf{\Phi }}^d{\mathbf{z}}^d|{|}_2^2+2\sqrt{\lambda }\sum _d||{\mathbf{z}}^d|{|}_2.$$

$$(\sum _d{a}_d{\mathbf{K}}^d+\lambda {\mathbf{I}}_n)\mathbf{c}=\mathbf{y}.$$

$$P(\mathbf{c},\mathbf{a}):=||\mathbf{y}-\sum _{ji}{a}_{(ji)}{\mathbf{K}}^{(ji)}\mathbf{c}|{|}_2^2+\lambda \sum _{ji}{a}_{(ji)}{\mathbf{c}}^T{\mathbf{K}}^{(ji)}\mathbf{c}+\sum _j\sqrt{\sum _i{a}_{(ji)}^2}$$

Full text

1 Introduction

Traditional methods for forecasting stationary multivariate time series from their own past are derived from the classical linear ARMA modelling. In these, the prediction of the next point in the future of the series is constructed as a linear function of the past observations. The use of linear functions as the predictors is in part based on the Wold representation theorem (e.g. [6]) and in part, probably more importantly, on the fact that the linear predictor is the best predictor (in the mean-square-error sense) in case the time series is Gaussian.

The Gaussian assumption is therefore often adopted in the analysis of time series to justify the simple linear modelling. However, it is indeed a simplifying assumption since for non-Gaussian series the best predictor may very well be a non-linear function of the past observations. A number of parametric non-linear models has been proposed in the literature, each adapted to capture specific sources of non-linearity (for example multiple forms of regime-switching models, e.g. [20]).

In this paper we adopt an approach that does not rely on such prior assumptions for the function form. We propose to learn the predictor as a general vector-valued function $\mathbf{f}$ that takes as input the past observations of the multivariate series and outputs the forecast of the unknown next value (vector).

We have two principal requirements on the function $\mathbf{f}$. The first is the standard prediction accuracy requirement. That is, the function $\mathbf{f}$ shall be such that we can expect its outputs to be close (in the squared error sense) to the true future observations of the process. The second requirement is that the function $\mathbf{f}$ shall have a structure that will enable the analysis of the relationships amongst the subprocesses of the multivariate series. Namely, we wish to understand how parts of the series help in forecasting other parts of the multivariate series, a concept known in the time-series literature as graphical Granger modelling [11, 9].

To learn such a function $\mathbf{f}$ we employ the framework of regularised learning of vector-valued functions in the reproducing kernel Hilbert space (RKHS) [17]. Learning methods based on the RKHS theory have previously been considered for time series modelling (e.g. [10, 19, 15]). Though, as Pillonetto et al. note in their survey [18], their adoption for the dynamical system analysis is not a commonplace.

A critical step in kernel-based methods for learning vector-valued functions is the specification of the operator-valued kernel that exploits well the relationships between the inputs and the outputs. A convenient and well-studied class of operator-valued kernels (e.g. in [7, 8, 12]) are those decomposable into a product of a scalar kernel on the input space (input kernel) and a linear operator on the output space (output kernel).

The kernel uniquely determines the function space within which the function $\mathbf{f}$ is learned. It thus has significant influence on both our objectives described above. Instead of having to choose the input and the output kernels a priori, we introduce a method for learning the input and output kernels from the data together with learning the vector-valued function $\mathbf{f}$.

Our method combines in a novel way the multiple-kernel learning (MKL) approach [14] with learning the output kernels within the space of positive semidefinite linear operators on the output space [12]. MKL methods for operator-valued kernels have recently been developed in [13] and [19]. The first learns a convex combination of a set of operator-valued kernels fixed in advance, the second combines a fixed set of input kernels with a single learned output kernel. To the best of our knowledge, ours is the first method in which the operator-valued kernel is learned by combining a set of input kernels with a set of multiple learned output kernels.

In accordance with our second objective stated above, we impose specific structural constraints on the function search space so that the learned function supports the graphical Granger analysis. We achieve this by working with matrix-valued kernels operating over input partitions restricted to single input scalar series (similar input partitioning has recently been used in [19]).

We impose diagonal structure on the output kernels to control the model complexity. Though this has a cost in the inability to model contemporaneous relationships, it addresses the strong over-parametrisation in a principled manner. It also greatly simplifies the final structure of the problem, which, in result, suitably decomposes into a set of smaller independent problems solvable in parallel.

We develop two forms of sparsity-promoting regularisation approaches for learning the output kernels. These are based on the $\ell_{1}$ and $\ell_{1}/\ell_{2}$ norms respectively and are motivated by the search for Granger-causality relationships. As to our knowledge, the latter has not been previously used in the context of MKL.

Finally, we confirm on experiments the benefits our methods can bring to forecasting non-Gaussian series in terms of improved predictive accuracy and the ability to recover hidden dynamic dependency structure within the time series systems. This makes them valid alternatives to the state-of-the-art graphical Granger techniques.

Notation

We use bold upper case and lower case letters for matrices and vectors respectively, and the plain letters with subscripts for their elements. For any matrix or vector the superscript T denotes its transpose. Vectors are by convention column-wise so that $\mathbf{x}=(x_{1},\ldots,x_{n})^{T}$ is the $n$-dimensional vector $\mathbf{x}$. $\mathbb{R},\mathbb{R}^{n},\mathbb{R}^{m\times n}$ are the sets of real scalars, $n$-dimensional vectors, and $m\times n$ dimensional matrices. $\mathbb{R}^{m\times n}_{+}$ is the set of non-negative matrices, $\mathbb{S}^{m}_{+}$ the set of positive semi-definite $m\times m$ matrices and $\mathbb{D}^{m}_{+}$ the set of non-negative diagonal matrices. $\mathbb{N}_{m}$ is the set of positive integers $\{1,\ldots,m\}$. For any vectors $\mathbf{x,y}\in\mathbb{R}^{n}$, $\langle\mathbf{x},\mathbf{y}\rangle,||\mathbf{x}||_{1},||\mathbf{x}||_{2}$ are the standard inner product, $\ell_{1}$ and $\ell_{2}$ norms in the real Hilbert spaces. For any square matrix $\mathbf{A}$, $\operatorname*{Tr}(\mathbf{A})$ denotes the trace. For any two matrices $\mathbf{A,B}\in\mathbb{R}^{m\times n}$, $\langle\mathbf{A},\mathbf{B}\rangle_{F}:=\operatorname*{Tr}(\mathbf{A}^{T}\mathbf{B})$ is the Frobenius inner product and $||\mathbf{A}||_{F}:=\sqrt{\langle\mathbf{A},\mathbf{A}\rangle_{F}}$ the Frobenius norm. $\langle.,.\rangle_{\mathcal{F}}$ and $||.||_{\mathcal{F}}$ are the inner product and norm in the Hilbert space $\mathcal{F}$.

2 Problem Formulation

Given a realisation of a discrete stationary multivariate time series process \big{\{}\mathbf{y}_{t}\in\mathcal{Y}\subseteq\mathbb{R}^{m}:t\in\mathbb{N}_{n}\}, our goal is to learn a vector-valued function $\mathbf{f}:\mathcal{Y}^{p}\to\mathcal{Y}$ that takes as input the $p$ past observations of the process and predicts its future vector value (one step ahead). The function $\mathbf{f}$ shall be such that

we can expect the prediction to be near (in the Euclidean distance sense) the unobserved future value 2. 2.

its structure allows to analyse if parts (subprocesses) of the series are useful for forecasting other subprocesses within the series or if some subprocesses can be forecast independently of the rest; in short, it allows Granger-causality analysis [11, 9].

For notational simplicity, from now on we indicate the output of the function $\mathbf{f}$ as $\mathbf{y}\in\mathcal{Y}\subseteq\mathbb{R}^{m}$ and the input as $\mathbf{x}\in\mathcal{X}\subseteq\mathbb{R}^{mp}$ (bearing in mind that $\mathcal{X}=\mathcal{Y}^{p}$ is in fact the $p$-th order Cartesian product of $\mathcal{Y}$ and that the inputs $\mathbf{x}$ and outputs $\mathbf{y}$ are the past and future observations of the same $m$-dimensional series). We also align the time indexes so that our data sample consists of input-output data pairs \big{\{}(\mathbf{y}_{t},\mathbf{x}_{t}):t\in\mathbb{N}_{n}\big{\}}.

Following the standard function learning theory, we will learn $\mathbf{f}\in\mathcal{F}$ by minimising the regularised empirical squared-error risk (with a regularization parameter $\lambda>0$)

[TABLE]

Here $\mathcal{F}$ is the reproducing kernel Hilbert space (RKHS) of $\mathbb{R}^{m}$-valued functions endowed with the norm $||.||_{\mathcal{F}}$ and the inner product $\langle.,.\rangle_{\mathcal{F}}$. The RKHS is uniquely associated with a symmetric positive-semidefinite matrix-valued kernel $\mathbf{H}:\mathcal{X}\times\mathcal{X}\to\mathbb{R}^{m\times m}$ with the reproducing property

[TABLE]

where the map $\mathbf{H}_{\mathbf{x}}:\mathcal{X}\to\mathbb{R}^{m\times m}$ is the kernel section of $\mathbf{H}$ centred at $\mathbf{x}$ such that $\mathbf{H}_{\mathbf{x}_{i}}(\mathbf{x}_{j})=\mathbf{H}(\mathbf{x}_{i},\mathbf{x}_{j})$ for all $(\mathbf{x}_{i},\mathbf{x}_{j})\in(\mathcal{X},\mathcal{X})$. From the classical result in [17], the unique solution $\mathbf{\widehat{f}}$ of the variational problem (1) admits a finite dimensional representation

[TABLE]

where the coefficients $\mathbf{c}_{t}\in\mathcal{Y}$ are the solutions of the system of linear equations

[TABLE]

where $\delta_{st}=1$ if $s=t$ and is zero otherwise.

2.1 Granger-causality Analysis

To study the dynamical relationships in time series processes, Granger [11] proposed a practical definition of causality based on the accuracy of least-squares predictor functions. In brief, for two time series processes $\mathbf{y}$ and $\mathbf{z}$, $\mathbf{y}$ is said to Granger-cause $\mathbf{z}$ ($\mathbf{y}\to\mathbf{z}$) if given all the other relevant information we can predict the future of $\mathbf{z}$ better (in the mean-square-error sense) using the history of $\mathbf{y}$ than without it.

Though the concept seems rather straightforward, there are (at least) three points worth considering. First, the notion is purely technical based on the predictive accuracy of functions with differing input sets; it does not seek to understand the underlying forces driving the relationships. Second, in practice the conditioning set of information needs to be reduced to all the available information instead of all the relevant information. Third, it only considers relationships between pairs of (sub-)processes and not the interactions amongst a set of series.

Eichler [9] extended the concept to multivariate analysis through graphical models. The discussion in the paper focuses on the notion of Granger non-causality rather than causality and describes the specific Markov properties (conditional non-causality) encoded in the graphs of Granger-causal relationships. In this sense, the absence of a variable in a set of inputs is more informative of the Granger (non-)causality than its presence. In result, graphical Granger methods are typically based on (structured) sparse modelling [4].

3 Function Space and Kernel Specification

The function space $\mathcal{F}$ within which $\mathbf{f}$ is learned is fully determined by the reproducing kernel $\mathbf{H}$. Its specification is therefore critical for achieving the two objectives for the function $\mathbf{f}$ defined in Sect. 2. We focus on the class of matrix-valued kernels decomposable into the product of input kernels, capturing the similarities in the inputs, and output kernels, encoding the relationships between the outputs.

To analyse the dynamical dependencies between the series, we need to be able to discern within the inputs of the learned function $\mathbf{f}$ the individual scalar series. Therefore we partition the elements of the input vectors according to the source scalar time series. In result, instead of a single kernel operating over the full vectors, we work with multiple partition-kernels, each of them operating over a single input series. We further propose to learn the partition-kernels by combining the MKL techniques with output kernel learning within the cone of positive semi-definite matrices.

More formally, the kernel we propose to use is constructed as a sum of kernels $\mathbf{H}=\sum_{j}^{m}\mathbf{H}^{(j)}$, where $m$ is the number of the individual scalar-valued series in the multivariate process (dimensionality of the output space $\mathcal{Y}$). Each $\mathbf{H}^{(j)}:\mathcal{X}^{(j)}\times\mathcal{X}^{(j)}\to\mathbb{R}^{m\times m}$ is a matrix-valued kernel that determines its own RKHS of vector-valued functions. The domains $\mathcal{X}^{(j)}\subseteq\mathbb{R}^{p}$ are sets of vectors constructed by selecting from the inputs $\mathbf{x}$ only the $p$ coordinates $i^{(j)}\in\mathbb{N}_{mp}$ that correspond to the past of a single scalar time series $j$.

[TABLE]

Further, instead of choosing the individual matrix-valued functions $\mathbf{H}^{(j)}$, we propose to learn them. We construct each $\mathbf{H}^{(j)}$ again as a sum of kernels $\mathbf{H}^{(j)}=\sum_{i}^{s_{j}}\mathbf{H}^{(ji)}$ of possibly uneven number of summands $s_{j}$ of matrix-valued kernels $\mathbf{H}^{(ji)}:\mathcal{X}^{(j)}\times\mathcal{X}^{(j)}\to\mathbb{R}^{m\times m}$. For this lowest level $\mathbf{H}^{(ji)}$ we focus on the family of decomposable kernels $\mathbf{H}^{(ji)}=k^{(ji)}\,\mathbf{L}^{(ji)}$. Here, the input kernels $k^{(ji)}:\mathcal{X}^{(j)}\times\mathcal{X}^{(j)}\to\mathbb{R}$ capturing the similarity between the inputs are fixed in advance from a dictionary of valid scalar-valued kernels (e.g. Gaussian kernels with varying scales). The set \mathcal{L}=\big{\{}\mathbf{L}^{(ji)}:j=\mathbb{N}_{m},\,i=\mathbb{N}_{s_{j}},\sum_{j}^{m}s_{j}=l\big{\}} of output kernels $\mathbf{L}^{(ji)}:\mathcal{Y}\to\mathcal{Y}$ encoding the relations between the outputs is learned within the cone of symmetric positive semidefinite matrices $\mathbb{S}^{m}_{+}$.

[TABLE]

3.1 Kernel Learning and Function Estimation

Learning all the output kernels $\mathbf{L}^{(ji)}$ as full PSD matrices implies learning more than $m^{3}$ parameters. To improve the generalization capability, we reduce the complexity of the problem drastically by restricting the search space for $\mathbf{L}$’s to PSD diagonal matrices $\mathbb{D}^{m}_{+}$. This essentially corresponds to the assumption of no contemporaneous relationships between the series. We return to this point in Sect. 5.

As explained in Sect. 2.1, Granger (non-)causality learning typically searches for sparse models. We bring this into our methods by imposing a further sparsity inducing regularizer $\mathit{Q}:(\mathbb{R}^{m\times m})^{l}\to\mathbb{R}$ on the set of the output kernels $\mathcal{L}$. We motivate and elaborate suitable forms of $\mathit{Q}$ in Sect. 3.2.

The joint learning of the kernels and the function can now be formulated as the problem of finding the minimising solution $\mathbf{f}\in\mathcal{F}$ and $\mathbf{L}$’s $\in\mathbb{D}^{m}_{+}$ of the regularised functional

[TABLE]

where $\mathit{R}(\mathbf{f})$ is the regularised risk from (1). By calling on the properties of the RKHS, we reformulate this as a finite dimensional problem that can be addressed by conventional finite-dimensional optimisation approaches. We introduce the gram matrices $\mathbf{K}^{(ji)}\in\mathbb{S}^{n}_{+}$ such that $K^{(ji)}_{ts}=k^{(ji)}(\mathbf{x}^{(j)}_{t},\mathbf{x}^{(j)}_{s})$ for all $t,s\in\mathbb{N}_{n}$, the output data matrix $\mathbf{Y}\in\mathbb{R}^{n\times m}$ such that $\mathbf{Y}=(\mathbf{y}_{1},\ldots\mathbf{y}_{n})^{T}$, and the coefficient matrix $\mathbf{C}\in\mathbb{R}^{n\times m}$ such that $\mathbf{C}=(\mathbf{c}_{1},\ldots\mathbf{c}_{n})^{T}$.

Using these and (2) it is easy to show that the minimisation of the regularised risk $\mathit{R}(\mathbf{f})$ in (1) with respect to $\mathbf{f}\in\mathcal{F}$ is equivalent to the minimisation with respect to $\mathbf{C}\in\mathbb{R}^{n\times m}$ of the objective

[TABLE]

The finite dimensional equivalent of (5) is thus the joint minimisation of

[TABLE]

3.2 Sparse Regularization

The construction of the kernel $\mathbf{H}$ and the function space $\mathcal{F}$ described in Sect. 3 imposes on the function $\mathbf{f}$ the necessary structure that allows the Granger-causality analysis (as per our 2nd objective set-out in Sect. 2). As explained in Sect. 2.1, the other ingredient we need to identify the Granger non-causalities is sparsity within the structure of the learned function.

In our methods, the sparsity is introduced by the regularizer $\mathit{Q}$. By construction of the function space, we can examine the elements of the output kernels $\mathbf{L}^{(ij)}$ (their diagonals) to make statements about the Granger non-causality. We say the $j$-th scalar time series is non-causal for the $s$ series (given all the remaining series in the process) if $\mathbf{L}^{(ji)}_{ss}=0$ for all $i\in\mathbb{N}_{s_{j}}$.

Essentially, any of the numerous regularizers that exist for sparse or structured sparse learning [3] could be used as $Q$, possibly based on some prior knowledge about the underlying dependencies within the time-series process.

We elaborate here two cases that do not assume any special structure in the dependencies as the base scenarios. The first is the entry-wise $\ell_{1}$ norm across all the output kernels so that

[TABLE]

The second is the $\ell_{1}/\ell_{2}$ grouped norm

[TABLE]

After developing the learning strategy for these in Sect(s). 4.1 and 4.2, we provide some more intuition of their effects on the models and link to some other known graphical Granger techniques in Sect. 5.

4 Learning Strategy

First of all, we simplify the final formulation of the problem (7) in Sect. 3.1. Rather than working with a set of diagonal matrices $\mathbf{L}^{(ji)}$, we merge the diagonals into a single matrix $\mathbf{A}$. We then re-formulate the problem with respect to this single matrix in place of the set and show how this reformulation can be suitably decomposed into smaller independent sub-problems.

We develop fit-to-purpose approaches for our two regularisers in Sect(s). 4.1 and 4.2. The first - based on the decomposition of the kernel matrices into the corresponding empirical features and on the variational formulation of norms [3] - shows the equivalence of the problem with group lasso [22, 23]. The second proposes a simple alternating minimisation algorithm to obtain the two sets of parameters.

We introduce the non-negative matrix $\mathbf{A}\in\mathbb{R}^{l\times m}_{+}$ such that

[TABLE]

(each row in $\mathbf{A}$ corresponds to the diagonal of one output kernel; if $s_{j}=1$ for all $j$ we have $A_{js}=\mathbf{L}^{(j1)}_{ss}$). Using this change of variable, the optimisation problem (7) can be written equivalently as

[TABLE]

where

[TABLE]

and $\ddot{Q}(\mathbf{A})$ is the equivalent of $Q(\mathcal{L})$ so that

[TABLE]

and

[TABLE]

In equations (12) and (14) we somewhat abuse the notation by using $\sum_{ji}A_{(ji)s}$ to indicate the sum across the rows of the matrix $\mathbf{A}$.

From (12)-(14) we observe that, with both of our regularizers, problem (11) is conveniently separable along $s$ into the sum of $m$ smaller independent problems, one per scalar output series. These can be efficiently solved in parallel, which makes our method scalable to very large multivariate systems. The final complexity depends on the choice of the regulariser $\mathit{Q}$ and the appropriate algorithm. The overhead cost can be significantly reduced by precalculating the gram matrices $\mathbf{K}^{(ij)}$ in a single preprocessing step and sharing these in between the $m$ parallel tasks.

4.1 Learning with $\ell_{1}$ Norm

To unclutter notation we replace the bracketed double superscripts $(ij)$ by a single superscript $d=1,\ldots,l$. We also drop the regularization parameter $\tau$ (fix it to $\tau=1$) as it is easy to show that any other value can be absorbed into the rescaling of $\lambda$ and the $\mathbf{C}$ and $\mathbf{A}$ matrices. For each of the $s$ parallel tasks we indicate $\mathbf{A}_{:s}=\mathbf{a}$, $\mathbf{C}_{:s}=\mathbf{c}$ and $\mathbf{Y}_{:s}=\mathbf{y}$ so that the individual problems are the minimisations with respect to $\mathbf{a}\in\mathbb{R}^{l}_{+}$ and $\mathbf{c}\in\mathbb{R}^{n}$ of

[TABLE]

We decompose (for example by eigendecomposition) each of the gram matrices as $\mathbf{K}^{d}=\mathbf{\Phi}^{d}(\mathbf{\Phi}^{d})^{T}$, where $\mathbf{\Phi}^{d}\in\mathbb{R}^{n\times n}$ is the matrix of the empirical features, and we introduce the variables $\mathbf{z}^{d}=a_{d}(\mathbf{\Phi}^{d})^{T}\mathbf{c}\in\mathbb{R}^{n}$ and the set $\mathcal{Z}=\{\mathbf{z}^{d}:\mathbf{z}^{d}\in\mathbb{R}^{n},\,d\in\mathbb{N}_{l}\}$. Using these we rewrite111We extend the function $x^{2}/y:\mathbb{R}\times\mathbb{R}_{+}\to\mathbb{R}_{+}$ to the point $(0,0)$ by taking the convention $0/0=0$. equation (15)

[TABLE]

We first find the closed form of the minimising solution for $\mathbf{a}$ as $a_{d}=\sqrt{\lambda}||\mathbf{z}^{d}||_{2}$ for all $d$. Plugging this back to (16) we obtain

[TABLE]

Seen as a minimisation with respect to the set $\mathcal{Z}$ this is the classical group-lasso formulation with the empirical features $\mathbf{\Phi}^{d}$ as inputs. Accordingly, it can be solved by any standard method for group-lasso problems such as the proximal gradient descent method, e.g. [3], which we employ in our experiments. After solving for $\mathcal{Z}$ we can directly recover $\mathbf{a}$ from the above minimising identity and then obtain the parameters $\mathbf{c}$ from the set of linear equations

[TABLE]

The algorithm outlined above takes advantage of the convex group-lasso reformulation (17) and has the standard convergence and complexity properties of proximal gradient descent. The empirical features $\mathbf{\Phi}^{d}$ can be pre-calculated and shared amongst the $m$ tasks to reduce the overhead cost.

4.2 Learning with $\ell_{1}/\ell_{2}$ Norm

For the $\ell_{1}/\ell_{2}$ regularization, we need to return to the double indexation $(ji)$ to make clear how the groups are created. As above, for each of the $s$ parallel tasks we use the vectors $\mathbf{a},\mathbf{c}$ and $\mathbf{y}$. However, for vector $\mathbf{a}$ we will keep the $(ji)$ notation for its elements. The individual problems are the minimisations with respect to $\mathbf{a}\in\mathbb{R}^{l}_{+}$ and $\mathbf{c}\in\mathbb{R}^{n}$ of

[TABLE]

We propose to use the alternating minimisation with a proximal gradient step. At each iteration, we alternatively solve for $\mathbf{c}$ and $\mathbf{a}$. For fixed $\mathbf{a}$ we obtain $\mathbf{c}$ from the set of linear equations (18). With fixed $\mathbf{c}$, problem (19) is a group lasso for $\mathbf{a}$ with groups defined by the sub-index $j$ within the double $(ji)$ indexation of the elements of $\mathbf{a}$. Here, the proximal gradient step takes place to move along the descend direction for $\mathbf{a}$. Though convex in $\mathbf{a}$ and $\mathbf{c}$ individually, the problem (19) is jointly non-convex and therefore can converge to local minima.

5 Interpretation and Crossovers

To help the understanding of the inner workings of our methods and especially the effects of the two regularizers, we discuss here the crossovers to other existing methods for MKL and Granger modelling.

$\ell_{1}$ Norm

The link to group-lasso demonstrated in Sect. 4.1 is not in itself too surprising. The formulation in (15) can be recognised as a sparse multiple kernel learning problem which has been previously shown to relate to group-lasso (eg. [2], [21]). We derive this link in Sect. 4.1 using the empirical feature representation to i) provide better intuition for the structure of the learned function $\mathbf{\widehat{f}}$, ii) develop an efficient algorithm for solving problem (15).

The re-formulation in terms of the empirical features $\mathbf{\Phi}^{d}$ creates an intuitive bridge to the classical linear models. Each $\mathbf{\Phi}^{d}$ can be seen as a matrix of features generated from a subset $\mathcal{X}^{(j)}$ of the input coordinates relating to the past of a single scalar time series $j$. The group-lasso regularizer in equation (17) has a sparsifying effect at the level of these subsets zeroing out (or not) the whole groups of parameters $\mathbf{z}^{d}$. In the context of linear methods, this approach is known as the grouped graphical Granger modelling [16].

Within the non-linear approaches to time series modelling, Sindhwani et al. [19] recently derived a similar formulation. There the authors followed a strategy of multiple kernel learning from a dictionary of input kernels combined with a single learned output kernel (as opposed to our multiple output kernels). They obtain their IKL model, which is in its final formulation equivalent to problem (15), by fixing the output kernel to identity.

Though we initially formulate our problem quite differently, the diagonal constraint we impose on the output kernels essentially prevents the modelling of any contemporaneous relationships between the series (as does the identity output kernel matrix in IKL). What remains in our methods are the diagonal elements, which are non-constant and sparse, and which can be interpreted as the weights of the input kernels in the standard MKL setting.

$\ell_{1}/\ell_{2}$ Norm

The more complex $\ell_{1}/\ell_{2}$ regularisation discussed in Sect. 4.2 is to the best of our knowledge novel in the context of multiple kernel learning. It has again a strong motivation and clear interpretation in terms of the graphical Granger modelling. The norm has a sparsifying effect not only at the level of the individual kernels but at the level of the groups of kernels operating over the same input partitions $\mathcal{X}^{(j)}$. In this respect our move from the $\ell_{1}$ to the $\ell_{1}/\ell_{2}$ norm has a parallel in the same move in linear graphical Granger techniques. The $\ell_{1}$ norm Lasso-Granger method [1] imposes the sparsity on the individual elements of the parameter matrices in a linear model, while the $\ell_{1}/\ell_{2}$ of the grouped-Lasso-Granger [16] works with groups of the corresponding parameters of a single input series across the multiple lags $p$.

6 Experiments

To document the performance of our method, we have conducted a set of experiments on real and synthetic datasets. In these we simulate real-life forecasting exercise by splitting the data into a training and a hold-out set which is unseen by the algorithm when learning the function $\mathbf{\widehat{f}}$ and is only used for the final performance evaluation.

We compare our methods with the output kernel $\ell_{1}$ regularization (NVARL1) and and $\ell_{1}/\ell_{1}$ (NVARL12) with simple baselines (which nevertheless are often hard to beat in practical time series forecasting) as well as with the state-of-the-art techniques for forecasting and Granger modelling. Namely, we compare with simple mean and univariate linear autoregressive models (LAR), multivariate linear vector autoregressive model with $\ell_{2}$ penalty (LVARL2), the group-lasso Granger method [16] (LVARL1), and a sparse MKL without the $\mathcal{X}^{(j)}$ input partitioning (NVAR). Of these, the last two are the most relevant competitors. LVARL1, similarly to our methods, aims at recovering the Granger structure but is strongly constrained to linear modelling only. NVAR has no capability to capture the Granger relationships but, due to the lack of structural constraints, it is the most flexible of all the models.

We evaluate our results with respect to the two objectives for the function $\mathbf{\widehat{f}}$ defined in Sect. 2. We measure the accuracy of the one-step ahead forecasts by the mean square error (MSE) for the whole multivariate process averaged over 500 hold-out points. The structural objective allowing the analysis of dependencies between the sub-processes is wired into the method itself (see Sect(s). 3 and 2.1) and is therefore satisfied by construction. We produce adjacency matrices of the graphs of the learned dependencies, compare these with the ones produced by the linear Granger methods and comment on the observed results.

6.1 Technical Considerations

For each experiment we preprocessed the data by removing the training sample mean and rescaling with the training sample standard deviation. We fix the number of kernels for each input partition to six ($s_{j}=6$ for all $j$) and use the same kernel functions for all experiments: a linear, 2nd order and 3rd polynomial, and Gaussian kernels with width 0.5, 1 and 2. We normalise the kernels so that the training Gram matrices have trace equal to the size of the training sample.

We search for the hyper-parameter $\lambda$ by a 5-fold cross-validation within a 15-long logarithmic grid $\lambda\in\{10^{-3},\ldots,10^{4}\}\sqrt{n}l$, where $n$ is the training sample size and $l$ is the number of kernels or groups (depending on the method). In each grid search, we use the previous parameter values as warm starts. We do not perform an exhaustive search for the optimal lag for each of the scalar input series by some of the classical testing procedures (based on AIC, BIC etc.). We instead fix it to $p=5$ for all series in all experiments and rely on the regularization to control any excess complexity.

We implemented our own tools for all the tested methods based on variations of proximal gradient descent with ISTA line search [5]. The full Matlab code is available at https://bitbucket.org/dmmlgeneva/nonlinear-granger

6.2 Synthetic Experiments

We have simulated data from a five dimensional non-Gaussian time-series process generated through a linear filter of a 5-dimensional i.i.d. exponential white noise $\mathbf{e}_{t}$ with identity covariance matrix (re-centered to zero and re-scaled to unit variance). The matrix $\mathbf{\Psi}=[0.7,1.3,0,0,0;0,0.6,-1.5,0,0;0,-1.2,1.46,0,0;0,0,0,0.6,1.4;0,0,0,1.3,-0.5]$ in the filter $\mathbf{y}_{t}=\mathbf{e}_{t}+\mathbf{\Psi}\mathbf{e}_{t-1}$ is such that the process consists of two independent internally interrelated sub-processes, one composed of the first 3 scalar series, the other of the remaining two series. This structural information, though known to us, is unknown to the learning methods (not considered in the learning process).

We list in Table 6.2 the predictive performance of the tested methods in terms of the average hold-out MSE based on training samples of varying size. Our methods clearly outperform all the linear models. The functionally strongly constrained linear LVARL1 performs roughly on par with our methods for the small sample sizes. But for larger sample sizes, the higher flexibility of the function space in our methods yields significantly more accurate forecasts (as much as 10% MSE improvement).