$p$-Markov Gaussian Processes for Scalable and Expressive Online Bayesian
Nonparametric Time Series Forecasting

Yves-Laurent Kom Samo<EMAIL_ADDRESS>Stephen J. Roberts
<EMAIL_ADDRESS>
University of Oxford University of Oxford

###### Abstract

In this paper we introduce a novel online time series forecasting model we
refer to as the $p$M-GP filter. We show that our model is equivalent to
Gaussian process regression, with the advantage that both online forecasting
and online learning of the hyper-parameters have a constant (rather than
cubic) time complexity and a constant (rather than squared) memory requirement
in the number of observations, without resorting to approximations. Moreover,
the proposed model is expressive in that the family of covariance functions of
the implied latent process, namely the spectral Matérn kernels, have recently
been proven to be capable of approximating arbitrarily well any translation-
invariant covariance function. The benefit of our approach compared to
competing models is demonstrated using experiments on several real-life
datasets.

##  1 INTRODUCTION

Interest in analysing and forecasting time series is thousands of years old.
This field of research has witnessed two breakthroughs over the second-half of
the twentieth century in the development of the Box-Jenkins methodology (Box
et al. (1970)) and the derivation of the Kalman filter (Kalman (1960)). While
the ARIMA model underpinning the Box-Jenkins methodology has since become a
standard in graduate curricula worldwide, the Kalman filter has found a great
number of applications in addition to time series forecasting, including, but
not limited to, object tracking, navigation and computer vision.

###  1.1 Popular Approaches

The Box-Jenkins methodology for analysing time series consists of three steps.
Firstly, trends and seasonalities are removed from the training dataset by
differentiating the time series iteratively111That is the following time
series are constructed from the original time series $y_{k}$:
$\Delta^{(1)}y_{k}:=y_{k}-y_{k-1}$, …,
$\Delta^{(i+1)}y_{k}:=\Delta^{(i)}y_{k}-\Delta^{(i)}y_{k-1}$ etc… until an
ergodic-stationary time series is found. Secondly, the ergodic-stationary time
series is assumed to be a realisation of an ARMA process whose parameters are
estimated, typically by maximizing the likelihood. Finally, the calibrated
model is checked by analysing the sample autocorrelation of residuals.

Both the ARIMA model and the Kalman filter can be represented as discrete time
linear Gaussian state space models (see (Durbin, Koopman, 2012, chap. 3) or
(Commandeur, Koopman, 2007, chap. 10)). Discrete time state space models aim
at inferring the state $x_{k}\in\mathbb{R}^{p}$ of a latent dynamic system,
that is assumed to evolve according to Markovian dynamics, from some noisy
observations $y_{k}\in\mathbb{R}^{q}$. The models are characterised by a
probability distribution for the initial latent state of the system
$p(x_{0})$, a Markovian transition dynamics $p(x_{k}|x_{k-1})$, and a
measurement model $p(y_{k}|x_{k})$. It is also assumed that the observations
are independent from each other and from other latent states conditional on
their associated latent states (i.e. $\forall i\neq j~{}y_{i}\perp
y_{j}|x_{i}$ and $y_{i}\perp x_{j}|x_{i}$). The forecasting problem then
consists of recursively determining $p(x_{k+n}|y_{1},\dots,y_{k})$ for some
$n>0$.

In linear Gaussian state space models (LGSSM), the three distributions above
are taken to be Gaussian, and the transition and measurement distributions
have covariance matrices that do not depend on the state process. Moreover, it
is assumed that

$\displaystyle\text{E}\left(x_{k}|x_{k-1}\right)=F_{k}x_{k-1}+B_{k}u_{k},E\left(y_{k}|x_{k}\right)=H_{k}x_{k},$

where $u_{k}$, $B_{k}$, $F_{k}$ and $H_{k}$ are given (or estimated off-line).
Under this class of models the forecasting problem can be solved exactly and
very efficiently using the Kalman filter equations.

LGSSMs can be extended in several ways. Alternatives have been proposed,
including the popular extended Kalman filter, that relax the linearity
assumption in LGSSMs by postulating that

$\displaystyle\text{E}\left(x_{k}|x_{k-1}\right)=f(x_{k-1},u_{k})\text{,
}E\left(y_{k}|x_{k}\right)=h(x_{k}),$

for some (possibly nonlinear) functions $f$ and $h$, while retaining the
Gaussian assumptions on the initial, transition and measurement distributions.
Other approaches extend LGSSMs by relaxing the foregoing Gaussian assumptions.
For instance, discrete measurement probability distributions such as the
Poisson distribution or the negative binomial distribution would allow for
count (as opposed to real-valued) observations (West et al. (1985)), while
leptokurtic probability distributions such as the student-t distribution or
the Laplace distribution, used as measurement distribution, might improve the
robustness of forecasts to outliers. For a more comprehensive review of Box-
Jenkins models and state space methods for time series analysis, we refer the
reader to Commandeur, Koopman (2007), Durbin, Koopman (2012) and Hamilton
(1994).

More generally, time series forecasting can be regarded as a regression
problem, where one is interested in predicting future values and associated
confidence bounds from historical observations. Gaussian process regression or
GPR (Rasmussen, Williams, 2005, chap. 2) provides a flexible Bayesian
nonparametric framework for that purpose. The GPR approach to time series
modelling consists of regarding the time series as arising from the sum of a
continuous time Gaussian process $(z_{t})_{t\geq 0}$ with mean function $m$
and covariance function $k_{\theta}$, and an independent Gaussian white noise:

$\displaystyle(z_{t})_{t\geq
0}\sim\mathcal{GP}\left(m(.),k_{\theta}(.,.)\right),~{}(\epsilon_{t})_{t\geq
0}\sim\mathcal{GWN}(\sigma^{2}),$ $\displaystyle\forall
t_{i},~{}y_{t_{i}}=z_{t_{i}}+\epsilon_{t_{i}}.$

As a result, $(y_{t})_{t\geq 0}$ is also a Gaussian process (GP), and the
noise variance $\sigma^{2}$ and all other hyper-parameters can easily be
inferred from some historical data $\\{y_{t_{0}},\dots,y_{t_{k}}\\}$ by
maximizing the corresponding multivariate Gaussian likelihood. Predictive
distributions may then be obtained in closed-form.

More recently, online passive-aggressive algorithms have been proposed by
Crammer et al. (2006), that allow for online linear and basis function
regression with strong worst-case loss bounds.

###  1.2 Limitations of Popular Approaches

The critical assumption underpinning the Box-Jenkins methodology is that any
time series, after enough iterative differentiations, will become ergodic-
stationarity, or more precisely will pass standard ergodic-stationarity
statistical tests. The fundamental problem with this approach is that in
practice, one will have access to finite samples that might not be
sufficiently large or informative for the time series to pass any ergodic-
stationarity test. It is in this spirit that (Durbin, Koopman, 2012, §3.10.1)
noted that ‘in the economic and social fields, real series are never
stationary however much differencing is done’. The ergodic-stationarity
assumption is not a major limitation per se, but relying on the assumption
that sufficient data have been collected to fully characterise the latent
process (including testing for ergodic-stationarity) is a major limitation of
the Box-Jenkins methodology. In that regards, Kwiatkowski et al. (1992) noted
that ‘most economic time series are not very informative about whether or not
there is a unit root’. The state space approach on the other hand does not
require that the sample be informative enough to test for stationarity, which
is why Commandeur, Koopman (2007) and Durbin, Koopman (2012) argued it should
be preferred to the Box-Jenkins methodology.

Additionally, the scalability of state space models is very appealing.
However, the dynamics of the latent time series, characterised by $H_{k}$ ($h$
in the nonlinear case) is usually assumed to be known. Although this
assumption is fairly mild in many engineering applications where the dynamics
are derived from the laws of physics, when analysing time series arising from
complex systems where the dynamics are not known, $H_{k}$ (resp. $h$) fully
characterize the covariance function of the latent time series $(x_{t})_{t\geq
0}$. Thus, hand-picking $H_{k}$ or $h$ would be problematic as the covariance
function of the latent process should be flexibly learned from the data to
appropriately uncover and exploit patterns. Unfortunately, no state space
dynamics has been proposed in the literature, to the best of our knowledge,
that guarantees that the implied covariance functions of the latent processes
$(x_{t})_{t\geq 0}$ can approximate arbitrarily well any stationary covariance
function.

Families of covariance functions that are dense in the class of all continuous
bounded covariance functions have recently been proposed by Kom Samo, Roberts
(2015a), and can readily be used within the GPR framework. Unfortunately, with
a cubic time complexity and a squared memory requirement, GPR scales poorly
with the number of training samples. Sparse approximations have been proposed
that scale linearly with the number of observations (Quinonero-Candela et al.
(2005); Lazaro-Gredilla et al. (2010)). However, those methods are not good
enough for streaming applications as they require loss of information through
windowing to be of practical interest.

Passive-aggressive (PA) algorithms (Crammer et al. (2006)) have constant time
complexity and memory requirement for online linear regression problems.
However, their kernelized versions have memory requirement and computational
time that depend on the number of ‘support vectors’, and might increase
unboundedly with the number of observations (Wang, Vucetic (2010)). Thus, PA
algorithms in their original form are not appropriate for flexibly forecasting
time series in a high throughput setting. Wang, Vucetic (2010) have introduced
PA algorithms with an additional memory or CPU ‘budget’ constraint, but the
authors restricted themselves to classification problems. Moreover, unlike
probabilistic approaches, PA algorithms do not provide uncertainty around
predictions. A Bayesian PA-like online learning approach has been proposed by
Shi, Zhu (2014). However, like all other kernel-PA algorithms, the kernel is
assumed to be given and online learning of its parameters is not addressed by
the authors.

Our contibutions: In this paper we build on the work of Kom Samo, Roberts
(2015a) to propose a state space model for analysing and forecasting time
series that we prove is equivalent to GPR, under a family of kernels for the
latent GP that is dense in the family of all stationary kernels, and that
allows decoupling the flexibility of the kernel from the degree of smoothness
of the latent GP. More importantly, we derive a novel scheme for online
learning of hyper-parameters as a solution to a convex optimization problem,
which happens to be a passive-aggressive algorithm (Crammer et al. (2006))
with stochastic gradient descent updates (Bottou (1998)). Overall, our model
has constant time complexity and constant memory requirement, both for online
learning of the hyper-parameters and for online prediction.

The rest of the paper is structured as follows. In section 2 we provide the
intuition behind our approach and the related mathematical background. In
section 3 we formally introduce our model. We illustrate that our approach
outperforms competing approaches with some experiments in section 4. Finally,
we discuss possible extensions and conclude in section 5.

##  2 INTUITION AND BACKGROUND

In this paper we consider modelling real-valued time series. Our working
assumption is that each sample $\\{y_{t_{0}},\dots,y_{t_{N}}\\}$ of a time
series arises as noisy observations of a ‘smooth’222By smooth we mean that the
process is at least mean square continuous, and higher order derivatives might
exist. latent stochastic process $(z_{t})_{t\geq 0}$ pertubated by an
independent white noise $(\epsilon_{t})_{t\geq 0}$:

$\displaystyle\forall t,~{}y_{t}=z_{t}+\epsilon_{t},$ (1)
$\displaystyle(\epsilon_{t})_{t\geq
0}\sim\mathcal{WN}(\sigma^{2}),~{}(\epsilon_{t})_{t\geq 0}\perp(z_{t})_{t\geq
0}.$ (2)

In financial econometrics applications for instance, the process
$(z_{t})_{t\geq 0}$ may represent the logarithm of the equilibrium or fair
price of an asset, while $(\epsilon_{t})_{t\geq 0}$ accounts for short term
shocks, microstructure noise and so-called ‘fat-fingers’. More generally,
$(z_{t})_{t\geq 0}$ denotes the value of the time series of interest, after
taking out hazards such as side effects, measurement noise, recording errors
and so forth, accounted for by $(\epsilon_{t})_{t\geq 0}$.

Moreover, we will assume that the latent process $(z_{t})_{t\geq 0}$ is a
trend-stationary Gaussian process, i.e. a Gaussian process with a translation-
invariant covariance function and an arbitrary mean function:

$\displaystyle(z_{t})_{t\geq 0}\sim\mathcal{GP}(m(.),k(.)).$ (3)

It is our working assumption that the class of trend-stationary Gaussian
processes is large enough for our purpose. We refer the reader to Appendix A
for a discussion on why this is a relatively mild assumption. From a Bayesian
perspective, we are assuming a priori that the process evolves about a
deterministic trend, and that our prior uncertainty on how the process might
deviate from the trend is invariant by translation. In applications, the trend
$m$ will typically be taken to lie in a parametric family of functions.

Furthermore, we assume that the white noise process $(\epsilon_{t})_{t\geq 0}$
is Gaussian:

$\displaystyle(\epsilon_{t})_{t\geq 0}\sim\mathcal{GWN}(\sigma^{2}).$ (4)

We will discuss possible extensions to more robust noise models in section 5.
Our setup thus far is the same as GPR with a translation-invariant covariance
function. Our next objective is to investigate whether we could exploit the
scalability of the state space approach by providing an ‘equivalent’ state
space representation. More precisely, by ‘equivalent’ we mean one in which the
observable is the noisy time series value $y_{t}\in\mathbb{R}$, and for every
collection of observation times $\\{t_{0},\dots,t_{N}\\}$, the joint
distribution over $\left(y_{t_{0}},\dots,y_{t_{k}}\right)$ as per the state
space representation is the same as that implied by Equations (1) to (4).
Thus, we need to construct an appropriate Markov process $(x_{t})_{t\geq
0},x_{t}\in\mathbb{R}^{p}$ that is related to $(z_{t})_{t\geq 0}$ and from
which the state space dynamics will be deduced. We will then derive the
measurement dynamics from Equation (1). As $(x_{t})_{t\geq 0}$ should be a
Markov process, intuitively $x_{t}$ should contain all information about
present and past values of the latent time series $(z_{t})_{t\geq 0}$ if we
want the state space representation to be equivalent to Equations (1) to (4).
Thus, it seems natural to consider that $(z_{t})_{t\geq 0}$ is a coordinate
process of the vector state process $(x_{t})_{t\geq 0}$.

Case 1: $\forall t,~{}x_{t}=z_{t}$

If we take the state process to be the latent time series itself,
$(z_{t})_{t\geq 0}$ needs to be Markovian. However, this would imply that the
latent time series will either not be mean square differentiable, or, as the
following lemma states, it will almost surely be equal to the trend plus a
constant, which might be too restrictive.

###### Lemma 1

If a real-valued trend-stationary Gaussian process with differentiable mean
function is mean square differentiable and Markovian, then it has a constant
covariance function (or equivalently it is almost surely equal to its mean
function plus a constant).

Proof  See Appendix B.
Intuitively, the memorylessness of the Markov assumption conflicts with the
memoryfulness of the differentiability assumption, so that the choice
$x_{t}=z_{t}$ might only be appropriate when assuming that the smooth latent
time series $(z_{t})_{t\geq 0}$ is mean square continuous but not
differentiable. An example covariance function yielding a continuous trend-
stationary Markov Gaussian process is the Matérn-$\frac{1}{2}$ kernel
$k_{\text{exp}}$ (see Rasmussen, Williams, 2005, Appendix §B.2). The
equivalent state space representation is easily derived from standard Gaussian
identities as:

$\begin{cases}z_{t_{0}}\sim\mathcal{N}\left(0,k_{\text{exp}}(0)\right),\\\
z_{t_{k}}=m(t_{k})+\alpha_{k}z_{t_{k-1}}+\sqrt{k_{\text{exp}}(0)\left(1-\alpha_{k}^{2}\right)}\xi_{t_{k}}\\\
y_{t_{k}}=z_{t_{k}}+\epsilon_{t_{k}}\end{cases}$

where $\alpha_{k}=\frac{k_{\text{exp}}(t_{k},t_{k-1})}{k_{\text{exp}}(0)}$,
and $(\xi_{t})_{t\geq 0}\sim\mathcal{GWN}(1)$.

Case 2: $\forall t,~{}x_{t}=\left(z_{t},z_{t}^{(1)},\dots,z_{t}^{(p)}\right)$

In many applications, we might be interested in postulating that the latent
time series $(z_{t})_{t\geq 0}$ is smoother, for instance $p$ times mean
square differentiable. In this case, we derive our intuition for what might be
good candidates to complete the state process into an appropriate Markov
vector process from the theory of analytic functions underpinning
implementations of scientific functions in most programming languages.
Conceptually, the aforementioned theory implies that under mild conditions,
the value of an infinitely smooth function and its derivatives
$f(t_{k}),f^{(1)}(t_{k}),\dots,f^{(i)}(t_{k}),\dots$ at a given point $t_{k}$
are sufficient to fully characterise the entire function elsewhere on the
domain. Although we only require $(z_{t})_{t\geq 0}$ to be $p$ times
differentiable, it is our hope that by augmenting the state process with the
$p$ consecutive derivatives of $(z_{t})_{t\geq 0}$, we could ensure that (i)
each state variable $x_{t}$ contains all information about past values of the
latent time series $(z_{u})_{u\leq t}$ or equivalently the state process
$(x_{t})_{t\geq 0}$ is Markovian, (ii) the latent time series $(z_{t})_{t\geq
0}$ is as smooth as desired, and (iii) the covariance function of the latent
time series is not strongly restricted.

When $(z_{t})_{t\geq 0}$ is a Gaussian process that is $p$ times continuously
differentiable in the mean square sense, we denote as p-derivative Gaussian
process ($p$DGP) the vector Gaussian process $(x_{t})_{t\geq 0}$ where
$x_{t}=\left(z_{t},z^{(1)}_{t},\dots,z^{(p)}_{t}\right)$ and $(z^{(i)})_{t\geq
0}$ is the $i$-th order mean square derivative of $(z_{t})_{t\geq 0}$. We say
that $(z_{t})_{t\geq 0}$ is a $p$-Markov Gaussian process ($p$M-GP) when it is
$p$ times continuously differentiable and the corresponding $p$DGP is
Markovian.333A $0$M-GP ($p=0$) is a Gaussian Markov process. Stationary
$p$M-GPs have been extensively studied in the seminal paper by Doob
(1944).444In Doob’s paper a $p$M-GP is denoted a t.h.G.M.(p-1) process. In
particular, the fundamental theorem below (see Doob, 1944, Theorem 4.9 ii)
provides necessary and sufficient conditions on our covariance function $k$
for the augmented state process to be Markovian when the trend is constant.

###### Theorem 2

A real-valued stationary Gaussian process $(z_{t})_{t\geq 0}$ with integrable
covariance function $k$ is $p$-Markov if and only if $k$ admits a spectral
density of the form

$\forall s>0,~{}S(s):=\int k(\tau)e^{-2\pi
is\tau}d\tau=\frac{c}{|A(is)|^{2}},$ (5)

where $c>0$ and $A$ is a polynomial of degree $p+1$ with real-valued
coefficients and no purely imaginary root.

###### Corollary 3

A stationary Gaussian process with covariance function the Matérn kernel

$k_{\text{ma}}(\tau;k_{0},l,\nu)=k_{0}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}|\tau|}{l}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}|\tau|}{l}\right),$

where $k_{0},l>0$, $\Gamma$ is the Gamma function, $K_{\nu}$ is the modified
Bessel function of second kind, and with $\nu=p+\frac{1}{2},~{}p\in\mathbb{N}$
is $p$-Markov.

Proof  The spectral density of a Matérn-$(p+\frac{1}{2})$ covariance function
is of the form

$\frac{c}{(\alpha+4\pi^{2}s^{2})^{p+1}}=\frac{c}{|(\sqrt{\alpha}+2\pi
is)^{p+1}|^{2}}$

with $c,\alpha>0$ (see Rasmussen, Williams, 2005).
Before deriving the state space representation of trend-stationary GPs with
Matérn-$(p+\frac{1}{2})$ covariance functions, we recall that a $p$DGP is a
vector Gaussian process and that in the trend-stationary case

$\forall~{}0\leq i,j\leq
p,~{}\text{cov}\left(z_{u}^{(i)},z_{v}^{(j)}\right)=(-1)^{j}k^{(i+j)}(u-v),$

where we use the superscript (i) to denote $i$-th order derivative or the
orginal function (or process) when $i=0$ (see Adler, Taylor, 2011, §2.6). If
we further denote

$K_{u,v}:=\left\\{\text{cov}\left(z_{u}^{(i)},z_{v}^{(j)}\right)\right\\}_{0\leq
i,j\leq p}$

the covariance matrix between $x_{u}$ and $x_{v}$, as

$K_{u|v}:=K_{u,u}-K_{u,v}K_{v,v}^{-1}K_{v,u}$

the auto-covariance matrix of $x_{u}$ conditional on $x_{v}$, and as $L_{u|v}$
any squared matrix satisfying555The Cholesky factorisation and the SVD provide
such a decomposition.

$K_{u|v}=L_{u|v}L_{u|v}^{T},$

then we obtain the following equivalent LGSSM representation (see Appendix C
for details):

$\begin{cases}x_{t_{0}}\sim\mathcal{N}(0,K_{t_{0},t_{0}})\\\
x_{t_{k}}=F_{t_{k}}x_{t_{k-1}}+L_{t_{k}|t_{k-1}}\xi_{t_{k}}\\\
y_{t_{k}}=m(t_{k})+H^{T}x_{t_{k}}+\epsilon_{t_{k}}\end{cases}$

with $F_{t_{k}}=K_{t_{k},t_{k-1}}K_{t_{k-1},t_{k-1}}^{-1},~{}H=(1,0,\dots,0),$
and where $(\xi_{t})_{t\geq 0}$ is now a $(p+1)$-dimensional standard Gaussian
white noise.

The idea of using derivative information to construct a scalable alternative
to GPs for MCMC inference has recently been developed in Kom Samo, Roberts
(2015b). The idea of representing GPR as a state space model with a Markovian
dynamics involving derivatives has been explored by Särkkä et al. (2013),
where the authors proposed approximating the spectral density of a given
covariance function by a function in the form of Equation (5). Unfortunately,
deriving such approximations might be tedious in certain cases. More
importantly, postulating that the covariance function has spectral density of
the form Equation (5) is still not flexible enough as it creates a link
between the degree of differentiability of the latent process and the number
of modes of the spectral density. To see why, we note that critical points (or
local extrema) of $\frac{c}{|A(is)|^{2}}$ and $|A(is)|^{2}$ are the same, and
that $|A(is)|^{2}$ being a polynomial of degree $2p+2$ in $s$, it can have at
most $2p+1$ critical points. This means that if for instance we would like the
latent function to be at most once differentiable ($p=1$), the approach of
Särkkä et al. (2013) would restrict the spectral density to have at most $3$
modes, rather than flexibly learning the number of modes from the data. This
also holds when the spectral density is approximated by a rational function in
the extension suggested by Särkkä et al. (2013), as the degree of the
numerator needs to be smaller than the degree of the denominator, and the
later is linked to the degree of differentiability of the latent GP.

Families of kernels, namely spectral Matérn kernels (Kom Samo, Roberts
(2015a)), do exist that are dense in the family of all stationary covariance
functions, and allow postulating or learning the degree of differentiability
of the corresponding GP. In the next section, we build on the state space
representation of $p$-Markov Matérn Gaussian processes to provide an exact
state space representation of our time series model (Equations (1) to (4))
when the latent time series $(z_{t})_{t\geq 0}$ is assumed to have a spectral
Matérn-($p+\frac{1}{2}$) covariance function (Kom Samo, Roberts (2015a)).
Moreover, we derive a novel passive-aggressive algorithm for online learning
of the hyper-parameters of the state space model, and that has constant memory
requirement and constant time complexity.

##  3 OUR MODEL

We start by introducing a state space time series model we refer to as the
$p$-Markov Gaussian process filter ($p$M-GP filter) that we prove is
equivalent to Gaussian process regression under a spectral Matérn kernel.

###  3.1 The $p$-Markov Gaussian Process Filter

###### Definition 4

($p$M-GP filter) We denote $p$-Markov Gaussian process filter ($p$M-GP filter)
a state space model of the form:

$\begin{cases}\perp\left\\{{}^{i}_{c}x_{t_{0}},{}^{i}_{s}x_{t_{0}}\right\\}_{i=0}^{n},\perp\left\\{({}^{i}_{c}\xi_{t})_{t\geq
0},({}^{i}_{s}\xi_{t})_{t\geq 0}\right\\}_{i=0}^{n}\\\ \forall
i,~{}{}^{i}_{c}x_{t_{k-1}}\perp{}^{i}_{c}\xi_{t_{k}},~{}{}^{i}_{s}x_{t_{k-1}}\perp{}^{i}_{s}\xi_{t_{k}}\\\
\forall
i,~{}{}^{i}_{c}x_{t_{k}}\perp\epsilon_{t_{k}},~{}{}^{i}_{s}x_{t_{k}}\perp\epsilon_{t_{k}}\\\
\forall
i,~{}{}^{i}_{c}x_{t_{0}},{}^{i}_{s}x_{t_{0}}\sim\mathcal{N}(0,{}^{i}K_{t_{0},t_{0}})\\\
\forall
i,~{}{}^{i}_{c}x_{t_{k}}={}^{i}F_{t_{k}}{}^{i}_{c}x_{t_{k-1}}+{}^{i}L_{t_{k}|t_{k-1}}{}^{i}_{c}\xi_{t_{k}}\\\
\forall
i,~{}{}^{i}_{s}x_{t_{k}}={}^{i}F_{t_{k}}{}^{i}_{s}x_{t_{k-1}}+{}^{i}L_{t_{k}|t_{k-1}}{}^{i}_{s}\xi_{t_{k}}\\\
y_{t_{k}}=H^{T}\sum_{i=0}^{n}\left(\cos(\omega_{i}t_{k}){}^{i}_{c}x_{t_{k}}+\sin(\omega_{i}t_{k}){}^{i}_{s}x_{t_{k}}\right)\\\
~{}~{}~{}~{}~{}+m(t_{k})+\epsilon_{t_{k}}\end{cases}$

where $\perp$ denotes mutual independence, $({}^{i}_{c}\xi_{t})_{t\geq 0}$ and
$({}^{i}_{s}\xi_{t})_{t\geq 0}$ are $(p+1)$-dimensional standard Gaussian
white noises, $(\epsilon_{t})_{t\geq 0}$ is a scalar Gaussian white noise with
variance $\sigma^{2}$, ${}^{i}K_{u,v}$ is the cross-covariance matrix of a
$p$DGP with Matérn kernel $k_{\text{ma}}(\tau;k_{0i},l_{i},p+\frac{1}{2})$,
and where $H,{}^{i}F_{t_{k}},{}^{i}L_{t_{k}|t_{k-1}}$ are as in the previous
section (replacing $K_{u,v}$ by ${}^{i}K_{u,v}$).

The following proposition establishes that the $p$M-GP filter is equivalent to
GPR under a spectral Matérn kernel, and a possibly non-constant mean function.

###### Proposition 5

Let $(\hat{z}_{t})_{t\geq 0}$ be a trend-stationary Gaussian process with mean
function $m$ and spectral Matérn covariance function

$k_{\text{sma}}(\tau)=\sum_{i=0}^{n}k_{\text{ma}}\left(\tau;k_{0i},l_{i},p+\frac{1}{2}\right)\cos(\omega_{i}\tau).$
(6)

Let $(\hat{\epsilon}_{t})_{t\geq 0}$ be a Gaussian white noise with variance
$\sigma^{2}$ that is independent from $(\hat{z}_{t})_{t\geq 0}$, and
$(\hat{y}_{t})_{t\geq 0}$ the process defined as

$\forall t\geq 0,~{}\hat{y}_{t}=\hat{z}_{t}+\hat{\epsilon}_{t}.$

Finally, let $(y_{t})_{t\geq 0}$ be the observation process of the $p$M-GP
filter with the same parameters
$m,n,p,\sigma,\\{k_{0i},l_{i},\omega_{i}\\}_{i=0}^{n}$. Then, the processes
$(y_{t})_{t\geq 0}$ and $(\hat{y}_{t})_{t\geq 0}$ have the same law, or
equivalently:

$\forall
t_{0}<\dots<t_{N},~{}(y_{t_{0}},\dots,y_{t_{N}})\sim(\hat{y}_{t_{0}},\dots,\hat{y}_{t_{N}}).$

Proof  See Appendix D.
The above proposition has profound implications. Firstly, the $p$M-GP filter
is the first equivalent state space representation of the Gaussian process
regression model (on a unidimensional input space) under a family of
covariance functions that is dense in the family of all stationary covariance
functions. This representation allows combining the great flexibility of GPR,
with the unmatched scalability of state space models, without resorting to
approximations. Secondly, the model allows controlling the degree of
differentiability of the latent time series independently from the flexibility
of the covariance function of the latent time series. Finally, unlike models
such as ARIMA that assume equally spaced observations, as any GPR model, the
$p$M-GP filter is inherently asynchronous and naturally copes with missing
data.

###  3.2 Solution to the Forecasting Problem

We note that the $p$M-GP filter can be rewritten as

$\begin{cases}\bm{x}_{t_{0}}\sim\mathcal{N}(0,\bm{K}_{t_{0},t_{0}})\\\
\bm{x}_{t_{k}}=\bm{F}_{t_{k}}\bm{x}_{t_{k-1}}+\bm{L}_{t_{k}|t_{k-1}}\bm{\xi}_{t_{k}}\\\
y_{t_{k}}=m(t_{k})+\bm{H}_{t_{k}}^{T}\bm{x}_{t_{k}}+\epsilon_{t_{k}}.\end{cases}$
(7)

Here, $\bm{x}_{t_{k}}$ (resp. $\bm{\xi}_{t_{k}}$) is the $2(p+1)(n+1)$ vector
obtained by stacking up the the vectors
$\\{\dots,{}^{i}_{c}x_{t_{k}},{}^{i}_{s}x_{t_{k}},\dots\\}$ (resp.
$\\{\dots,{}^{i}_{c}\xi_{t_{k}},{}^{i}_{s}\xi_{t_{k}},\dots\\}$). Similarly,
$\bm{K}_{t_{0},t_{0}}$ (resp.
$\bm{F}_{t_{k}},\bm{L}_{t_{k}|t_{k-1}},\bm{K}_{t_{k}|t_{k-1}}$) is the
$2(p+1)(n+1)\times 2(p+1)(n+1)$ block diagonal matrix whose $2i$-th and
$(2i+1)$-th diagonal blocks are both ${}^{i}K_{t_{0},t_{0}}$ (resp.
${}^{i}F_{t_{k}},{}^{i}L_{t_{k}|t_{k-1}},{}^{i}K_{t_{k}|t_{k-1}}$). Finally,
$\bm{H}_{t_{k}}$ is the $2(p+1)(n+1)$ vector with coordinates $0$ except at
indices multiple of $p+1$, and $\forall i\in[0\dots
n],~{}\bm{H}_{t_{k}}[2i(p+1)]=\cos(\omega_{i}t_{k})$ and
$\bm{H}_{t_{k}}[(2i+1)(p+1)]=\sin(\omega_{i}t_{k})$.

Thus the $p$M-GP filter is a LGSSM, and the forecasting problem can be solved
exactly, iteratively and in closed-form, and with memory requirement and time
complexity both constant in the total number of observations. Using standard
Kalman filter and Gaussian processes techniques (see Appendix E) we obtain:

$\begin{cases}\forall
t>t_{k-1},~{}\bm{x}_{t}|y_{t_{0}:t_{k-1}}&\sim\mathcal{N}(\bm{m}_{t}^{-},\bm{P}_{t}^{-})\\\
\forall
t>t_{k-1},~{}y_{t}|y_{t_{0}:t_{k-1}}&\sim\mathcal{N}(m(t)+\bm{H}_{t}^{T}\bm{m}_{t}^{-},v_{t}^{-})\\\
\forall
t>t_{k-1},~{}z_{t}|y_{t_{0}:t_{k-1}}&\sim\mathcal{N}(m(t)+\bm{H}_{t}^{T}\bm{m}_{t}^{-},v_{t}^{-}-\sigma^{2})\\\
\forall t\geq
t_{k-1},~{}\bm{x}_{t}|y_{t_{0}:t}&\sim\mathcal{N}(\bm{m}_{t},\bm{P}_{t})\\\
\forall t\geq
t_{k-1},~{}z_{t}|y_{t_{0}:t}&\sim\mathcal{N}(m(t)+\bm{H}_{t}^{T}\bm{m}_{t},v_{t})\end{cases}$
(8)

with $y_{t_{0}:t_{k-1}}=\\{y_{t_{0}},\dots,y_{t_{k-1}}\\}$,
$z_{t}=m(t)+\bm{H}_{t}^{T}\bm{x}_{t}$ and

$\displaystyle\underline{\text{Prediction step:}}$
$\displaystyle\bm{m}_{t_{0}}^{-}=0,~{}\bm{P}_{t_{0}}^{-}=\bm{K}_{t_{0},t_{0}}$
(9) $\displaystyle\forall k\geq
1,t>t_{k-1}\begin{cases}\bm{m}_{t}^{-}&=\bm{F}_{t}\bm{m}_{t_{k-1}}\\\
\bm{P}_{t}^{-}&=\bm{F}_{t}\bm{P}_{t_{k-1}}\bm{F}_{t}^{T}+\bm{K}_{t|t_{k-1}}\end{cases}$
$\displaystyle\underline{\text{Update step:}}$ $\displaystyle\forall
t\begin{cases}v_{t}^{-}&=\bm{H}_{t}^{T}\bm{P}_{t}^{-}\bm{H}_{t}+\sigma^{2}\\\
\bm{e}_{t}^{-}&=y_{t}-m(t)-\bm{H}_{t}^{T}\bm{m}_{t}^{-}\\\
\bm{G}_{t}&=\frac{1}{v_{t}^{-}}\bm{P}_{t}^{-}\bm{H}_{t}\\\
\bm{m}_{t}&=\bm{m}_{t}^{-}+\bm{e}_{t}^{-}\bm{G}_{t}\\\
\bm{P}_{t}&=\bm{P}_{t}^{-}-v_{t}^{-}\bm{G}_{t}\bm{G}_{t}^{T}\\\
v_{t}&=\bm{H}_{t}^{T}\bm{P}_{t}\bm{H}_{t}\end{cases}.$ (10)

###  3.3 Online Learning of Hyper-Parameters

Noting that $y_{t_{0}}\sim\mathcal{N}(m(t_{0}),v_{t_{0}}^{-})$, and that

$\log p(y_{t_{0}:t_{N}})=\log p(y_{t_{0}})+\sum_{k=1}^{T}\log
p(y_{t_{k}}|y_{t_{0}:t_{k-1}}),$

it follows that maximum likelihood inference of the hyper-parameters of $m$,
$\sigma$, and $\\{k_{0i},l_{i},\omega_{i}\\}_{i=0}^{n}$ may be achieved with
linear time complexity $\mathcal{O}(T)$ and with constant memory requirement
using Equations (8, 9, 10). However, to be of practical interest in streaming
applications, this time complexity requires loss of information through
windowing, which might be detrimental to forecasting performance. We herein
propose an alternative online approach that has constant time complexity.

So far we have assumed that the hyper-parameters are constant over time.
However, if we extend the model to iteration-specific hyper-parameters, the
solution to the forecasting problem (Equations (8, 9, 10)) will provably
remain unchanged. The intuition behind our approach to online learning of the
hyper-parameters is borrowed from maximum likelihood inference when the hyper-
parameters are constant, and online passive-aggressive algorithms for
regression, classification and semi-supervised learning (Crammer et al.
(2006); Wang, Vucetic (2010); Chang et al. (2010)). Let $\bm{\theta}_{t_{k}}$
be the vector of hyper-parameters (or log hyper-parameters for positive hyper-
parameters) prevailing at time $t_{k}$. We recall that when hyper-parameters
are constant, we may write

$\displaystyle\log p_{{\bm{\theta}}}\left(y_{t_{0}:t_{T}}\right)=\log
p_{{\bm{\theta}}}\left(y_{t_{0}}\right)+\sum_{k=1}^{T}\log
p_{{\bm{\theta}}}\left(y_{t_{k}}|y_{t_{0}:t_{k-1}}\right).$

Let us define

$\mathcal{L}^{l}_{t_{k}}(\bm{\theta}):=\begin{cases}\log
p_{\bm{\theta}}\left(y_{t_{0}}\right)&\text{if }k=0\\\ \log
p_{\bm{\theta}}\left(y_{t_{k}}|y_{t_{0}:t_{k-1}}\right)&\text{if }k\neq
0\end{cases}.$

At time $t_{k}$, $\mathcal{L}^{l}_{t_{k}}(\bm{\theta})$ represents the (log)
likelihood that $\bm{\theta}$ explains the new observation $y_{t_{k}}$ given
all previous observations (and hyper-parameters) and thus can be regarded as a
local utility function that we should aim to maximize. However, we should also
be mindful of the ‘utility’ of $\bm{\theta}$ in accounting for all previously
observed data $(y_{t_{0}},\dots,y_{t_{k-1}})$, which can be achieved by
ensuring that the new update is not be too far away from
$\bm{\theta}_{t_{k-1}}$. Deriving online learning updates as solutions to a
constrained optimization problem that provides a trade-off between retaining
information from previous iterations and locally reducing a loss function has
been advocated for a long time in the online learning literature (Littlestone
(1989); Kivinen, Warmuth (1997); Helmbold et al. (1999); Crammer et al.
(2006); Wang, Vucetic (2010); Chang et al. (2010)). The approach we adopt is
largely inspired by the PA-II algorithm of Crammer et al. (2006). We first
consider the problem:

$\displaystyle\begin{cases}\underset{\bm{\theta}}{\text{min}}~{}||\bm{\theta}-\bm{\theta}_{t_{k-1}}||^{2}+c_{k}\xi^{2}\\\
\text{ s.t.
}\max\left(-\epsilon-\mathcal{L}^{l}_{t_{k}}(\bm{\theta}),0\right)\leq\xi\end{cases},$
(11)

where $\epsilon\geq 0$ and $c_{k}>0$. $c_{k}$ can be regarded as an
aggressiveness parameter: the larger $c_{k}$ the more weight is given to
increasing the local utility compared to retaining information from previous
iterations. $\epsilon$ on the other hand is a margin or tolerance parameter
that allows for some degree of local sup-optimality, thereby helping to avoid
overfitting. Although this problem provides a suitable trade-off, it is
tedious to solve analytically, and it might not have a closed-form solution.
However, noting that the objective function of problem (11) aims at minimizing
$||\bm{\theta}-\bm{\theta}_{t_{k-1}}||$, we may replace
$\mathcal{L}^{l}_{t_{k}}(\bm{\theta})$ with its first order Taylor expansion
at $\bm{\theta}_{t_{k-1}}$. We adopt the resulting problem for online learning
of the hyper-parameters:

$\displaystyle\begin{cases}\bm{\theta}_{t_{k}}=\underset{\bm{\theta}}{\text{argmin}}~{}||\bm{\theta}-\bm{\theta}_{t_{k-1}}||^{2}+c_{k}\xi^{2}\\\
\text{ s.t.
}\max\left(-\epsilon-\hat{\mathcal{L}}^{l}_{t_{k}}(\bm{\theta}),0\right)\leq\xi\end{cases},$
(12)

with
$\hat{\mathcal{L}}^{l}_{t_{k}}(\bm{\theta})=\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})+\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})^{T}(\bm{\theta}-\bm{\theta}_{t_{k-1}})$.
The above optimization problem is convex and has closed-form solution (see
Appendix F for the derivation, and Appendix G for the derivation of
$\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta})$):

$\displaystyle\bm{\theta}_{t_{k}}=\bm{\theta}_{t_{k-1}}+c_{k}\frac{\max\left(-\epsilon-\mathcal{L}^{l}_{t_{k}}\left(\bm{\theta}_{t_{k-1}}\right),0\right)}{1+c_{k}\|\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\|^{2}}\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}}).$
(13)

Noting that the parameters do not change when
$\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})>-\epsilon$, it follows that
our learning scheme can be regarded as a passive-aggressive algorithm with
stochastic gradient descent update (Bottou (1998)). Moreover, as

$\displaystyle\hat{\mathcal{L}}^{l}_{t_{k}}(\bm{\theta}_{t_{k}})-\hat{\mathcal{L}}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})$
$\displaystyle=\max\left(-\epsilon-\mathcal{L}^{l}_{t_{k}}\left(\bm{\theta}_{t_{k-1}}\right),0\right)\frac{c_{k}\|\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\|^{2}}{1+c_{k}\|\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\|^{2}}\geq
0,$

each change in the hyper-parameters increases the approximate conditional log
likelihood $\hat{\mathcal{L}}^{l}_{t_{k}}$. Moreover, the one-side
$\epsilon$-insensitive loss term
$\max\left(-\epsilon-\mathcal{L}^{l}_{t_{k}}\left(\bm{\theta}_{t_{k-1}}\right),0\right)$
guards against overfitting when an update is required. Algorithm 1 summarizes
online forecasting and hyper-parameters learning under the $p$M-GP filter.

Algorithm 1 Forecasting with the $p$M-GP filter.

In: $\\{y_{t_{k}}\\}_{k\geq 0}$, average sampling frequency $F_{s}$.

Out: Predictive probability density functions
$\\{\dots,p(z_{t}|y_{t_{0}}:y_{t_{k}}),\dots\\},t\geq t_{k}$.

Set $\omega_{i}=\frac{1+i}{1+n}\pi F_{s}$, set the other parameters to $0$.

for all $(t_{k},y_{t_{k}})$ do

Evaluate $\bm{\theta}_{t_{k}}$ using Eq. (13).

Run the prediction step Eqs. (9) with $t=t_{k}$.

Run the update step Eqs. (10) with $t=t_{k}$.

For $t\geq t_{k}$ get $p(z_{t}|y_{t_{0}}:y_{t_{k}})$ from Eqs. (8, 9).

end for

In order to control the aggressiveness of our learning algorithm in a manner
that is consistent across datasets, we rewrite the aggressiveness parameter in
the form
$c_{k}:=c\frac{\|\bm{\theta}_{t_{k-1}}\|^{2}}{\left(\epsilon+\hat{\mathcal{L}}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\right)^{2}},c>0$,
that arises by normalizing the term
$||\bm{\theta}-\bm{\theta}_{t_{k-1}}||^{2}$ (resp. $\xi^{2}$) in the objective
function of problem (12) by $\|\bm{\theta}_{t_{k-1}}\|^{2}$ (resp.
$\left(\epsilon+\hat{\mathcal{L}}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\right)^{2}$).

##  4 EXPERIMENTS

Figure 1: CO2 experiment: (left) running normalized mean absolute error (NMAE), (middle) effect of the $p$M-GP aggressiveness parameter $c$ on NMAE, and (right) effect of the number of spectral components $K$ on NMAE. Table 1: Mean $\pm$ 1 standard deviation of normalized absolute errors in the experiments of section 4. | CO2 | Airline
---|---|---
$p$M-GP |  0.52 $\pm$0.70 |  0.52 $\pm$ 0.52
PAs AR(2) | 0.77 $\pm$ 0.57 | 0.94 $\pm$ 0.84
BLR AR(2) | 1.31 $\pm$ 0.58 | 0.77 $\pm$ 0.62
PAs AR(10) | 1.02 $\pm$ 0.57 | 0.92 $\pm$ 0.81
BLR AR(10) | 1.77 $\pm$ 1.00 | 0.77 $\pm$ 0.53

In this section we start by demonstrating that the $p$M-GP filter provides
considerably more accurate forecasts than competing fully-online approaches on
standard real-life time series datasets. We then experimentally illustrate the
sensitivity of the forecasting errors to the normalized aggressiveness
parameter $c$, the trend, and the number of spectral components $K$.

Benchmarking: We compare our model to competing fully-online alternatives on
the CO2 dataset of Rasmussen, Williams (2005), and the airline passengers
dataset of Box et al. (1970). We select as competing benchmarks two
autoregressive (AR) models, namely AR(2) and AR(10), and we use four different
algorithms to learn the autoregressive coefficients online, namely the PA,
PA-I and PA-II algorithms of Crammer et al. (2006), and Bayesian online linear
regression with i.i.d. standard normal priors on the weights (BLR). For
consistency with the $p$M-GP filter, we initialize the autoregressive weights
in all four algorithms to $0$. We choose $c=100.0,\epsilon=0.0$ for $p$M-GP
models and PA algorithms. We choose a linear trend and twice differentiability
($p=2$) for the $p$M-GP filter. BLR is run with a noise standard deviation of
$5\%$ of the sample standard deviation of the corresponding time series. For
each model we perform one-step ahead forecast. Mean absolute
errors666Excluding the first forecast, which is the same for all models.
normalized by the standard deviation of increments (NMAE) are reported in Tab.
1, where we refer to all three PA algorithms as PAs because they perform
identically up to two decimal points. The evolution of the running NMAE as a
function of time is illustrated in Fig. 1 (left) for the CO2 dataset. We note
that our approach provides more accurate forecasts than competing
alternatives.

Sensitivity to Parameters: The sensitivity of the accuracy of the $p$M-GP
filter to the trend, the aggressiveness parameter $c$ and the number of
spectral components is illustrated in Fig. 1 (middle and right) for the CO2
dataset. Overall, for our data size (607 points), it can be seen that the
error decreases as a function of $c$. This is in line with the empirical
observation of Crammer et al. (2006) (Fig. 5) that, for small datasets large
$c$ perform better, and $c=0.001$ begins to outperform $c=100.0$ for the PA
classification problems the author considered when the data size is in the
thousands. Moreover, we note from Fig. 1 (right) that for small datasets, a
large $K$ should not necessarily be preferred to a smaller one, as more
samples will typically be required to learn the model hyper-parameters.

##  5 DISCUSSION

We propose an exact state space representation of Gaussian process regression
for time series forecasting, under a family of kernels that is dense in the
family of all stationary kernels, and that allows decoupling the flexibility
of the covariance structure from the differentiability requirement of the
latent time series. When hyper-parameters are known, exact GP predictive
inference can be performed in constant time and with constant memory
requirement. Critically, we propose a novel passive-aggressive algorithm for
online learning of the model hyper-parameters. The overall approach we refer
to as the $p$M-GP filter has constant complexity and memory requirement, and
provides more accurate forecasts than fully-online competing alternatives on
the standard CO2 and airline passengers datasets. The approach may easily be
extended to structured time series prediction, thus allowing for online
learning of correlations between time series. Techniques that extend the
Kalman filter to leptokurtic measurement noise models (e.g. Agamennoni et al.
(2011)) may also be used out-of-the box to robustify our approach.

## References

  * Adler, Taylor (2011) Adler R. J., Taylor J. E. Topological Complexity of Smooth Random Functions: École D’Été de Probabilités de Saint-Flour XXXIX-2009. 2011.  (Lecture Notes in Mathematics / École d’Été de Probabilités de Saint-Flour).
  * Agamennoni et al. (2011) Agamennoni G., Nieto J. I., Nebot E. M. An Outlier-Robust Kalman Filter  // Robotics and Automation (ICRA), 2011 IEEE International Conference on. 2011. 1551–1558.
  * Bottou (1998) Bottou Léon.  Online Learning and Stochastic Approximations  // Online Learning in Neural Networks. 1998. 17, 9. 25.
  * Box et al. (1970) Box G. E. P., Jenkins G. M., Reinsel G. C., M. Ljung G. Time Series Analysis: Forecasting and Control. 1970.
  * Boyd, Vandenberghe (2004) Boyd S., Vandenberghe L. Convex Optimization. 2004.
  * Chang et al. (2010) Chang C.-C., Lee Y.J., Pao H.K. A Passive-Aggressive Algorithm for Semi-Supervised Learning  // International Conference on Technologies and Applications of Artificial Intelligence (TAAI). 2010. 335–341.
  * Commandeur, Koopman (2007) Commandeur J. J. F., Koopman S. J. State Space Time Series Analysis. 2007.
  * Crammer et al. (2006) Crammer K., Dekel O., Keshet J., Shalev-Shwartz S., Singer Y. Online Passive-Aggressive Algorithms  // Journal of Machine Learning Research. 2006. 7. 551–585.
  * Doob (1944) Doob J. L. The Elementary Gaussian Processes  // The Annals of Mathematical Statistics. 1944. 15, 3. 229–282.
  * Durbin, Koopman (2012) Durbin J., Koopman S. J. Time Series Analysis by State Space Methods. 2012.
  * Hamilton (1994) Hamilton J. D. Time series analysis. 1994.
  * Helmbold et al. (1999) Helmbold D. P., Kivinen J., Warmuth M. K. Relative Loss Bounds for Single Neurons  // Neural Networks, IEEE Transactions on. 1999. 10, 6. 1291–1304.
  * Kalman (1960) Kalman R. E. A New Approach to Linear Filtering and Prediction Problems  // Journal of Fluids Engineering. 1960. 82, 1. 35–45.
  * Kivinen, Warmuth (1997) Kivinen J., Warmuth M. K. Exponentiated Gradient Versus Gradient Descent for Linear Predictors  // Information and Computation. 1997. 132, 1. 1–63.
  * Kom Samo, Roberts (2015a) Kom Samo Y.-L., Roberts S. J. Generalized Spectral Kernels  // pre-print. 2015a.  arXiv:1506.02236.
  * Kom Samo, Roberts (2015b) Kom Samo Y.-L., Roberts S. J. String Gaussian Processes  // pre-print. 2015b.  arXiv:1507.06977.
  * Kwiatkowski et al. (1992) Kwiatkowski D., Phillips P. C. B., Schmidt P., Shin Y. Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root?  // Journal of econometrics. 1992. 54, 1. 159–178.
  * Lazaro-Gredilla et al. (2010) Lazaro-Gredilla M., Quinonero-Candela J., Rasmussen C. E., Figueiras-Vida A. R. Sparse Spectrum Gaussian Process Regression  // The Journal of Machine Learning Research. 2010. 11. 1866–1881.
  * Littlestone (1989) Littlestone N. Mistake Bounds and Logarithmic Linear-Threshold Learning Algorithms. 1989\.
  * Papoulis, Pillai (2002) Papoulis A., Pillai S. U. Probability, Random Variables, and Stochastic Processes. 2002.
  * Quinonero-Candela et al. (2005) Quinonero-Candela J., Rasmussen C. E., Herbrich R. A Unifying View of Sparse Approximate Gaussian Process Regression  // Journal of Machine Learning Research. 2005. 6. 1939–1959.
  * Rasmussen, Williams (2005) Rasmussen C. E., Williams C. K. I. Gaussian Processes for Machine Learning. 2005.
  * Särkkä et al. (2013) Särkkä S., Solin A., Hartikainen J. Spatiotemporal Learning via Infinite-Dimensional Bayesian Filtering and Smoothing: A Look at Gaussian Process Regression Through Kalman Filtering.  // IEEE Signal Process. Mag. 2013. 30, 4. 51–61.
  * Shi, Zhu (2014) Shi T., Zhu J. Online Bayesian Passive-Aggressive Learning  // International Conference on Machine Learning (ICML). 2014. 378–386.
  * Wang, Vucetic (2010) Wang Z., Vucetic S. Online Passive-Aggressive Algorithms on a Budget  // International Conference on Artificial Intelligence and Statistics. 2010. 908–915.
  * West et al. (1985) West M., Harrison P. J., Migon H. S. Dynamic Generalized Linear Models and Bayesian Forecasting  // Journal of the American Statistical Association. 1985. 80, 389. 73–83.

Appendix

Unless stated otherwise, the stochastic processes we consider throughout this
appendix are indexed in $\mathbb{R}^{+}$. To ease notations, we use the
superscript (i) to denote the $i$-th order derivative of a function or
stochastic process when it exists, or the original function or stochastic
process when $i=0$. In the case of stochastic processes, the derivative is to
be understood in the mean square sense. Moreover, stationarity of stochastic
processes is always meant in the weak sense. Furthermore, observation times
are always assumed to be distinct and sorted by index:
$t_{0}<\dots<t_{k}<\dots$.

## Appendix Appendix A

Objective: In this section we discuss the appropriateness of assuming that the
latent time series is a trend-stationary Gaussian process. More precisely, we
argue that given a single path of a time series on some bounded interval
$[0,T]$, neither the trend-stationarity assumption nor the Gaussianity
assumption can be invalidated.

Trend-Stationarity: Firstly, we note that unless further assumptions are made
about the trend other than it being smooth, trend-stationarity cannot be
invalidated using discrete observations of a single path as we can always find
an infinitely smooth trend or mean function $\hat{m}$ (using polynomials for
instance), that coincides with all discrete observations, making the observed
path highly likely to result from any trend-stationary stochastic process with
mean function $\hat{m}$.

In practice however, the trend might be assumed to lie in a parametric family.
Even in that case, the stationarity of the residual time series can hardly be
invalidated. In effect, most stationarity statistical tests have as null
hypothesis that the sample comes from a nonstationary time series. Evidence is
then gathered from the sample through the test statistics to determine whether
the null hypothesis can be rejected with some confidence, or equivalently if
there is enough evidence in the sample to conclude stationarity. Hence,
fundamentally, such an approach cannot be used to conclude nonstationarity, as
the lack of evidence of stationarity in a given sample is not an evidence of
nonstationarity of the underlying process. It might well be that the sample
does not contain enough information to fully characterise the underlying
stochastic process so that, had we collected more data, we might have been
able to conclude stationarity. It is in the same spirit that Kwiatkowski et
al. (1992) noted that ‘most economic time series are not very informative
about whether or not there is a unit root’. As we do not assume that we have
enough data to characterize the latent process, our trend-stationarity
assumption cannot be invalidated experimentally.

To illustrate our point, we drew a path from a stationary Gaussian process
with mean $0$ and squared exponential covariance function
$k(u,v)=e^{-\frac{1}{2}(u-v)^{2}}$ on $[0,10]$ discretized with mesh size
$0.001$ (see Figure 2).

Figure 2: Draw from a stationary GP on $[0,10]$.

We ran three standard stationarity statistical tests on the sample, namely the
Augmented Dickey-Fuller test (ADF), the ADF-GLS test and the Phillips-Perron
test.

Table 2: Results of stationarity tests on the time series in Figure 2. Test | Statistics | p-Value | Lags
---|---|---|---
ADF | -2.12 | 0.24 | 22
ADF-GLS | -1.02 | 0.28 | 22
Phillips-Perron | -1.64 | 0.46 | 22

As can be seen in Table 2, all three tests failed to find evidence for (to
conclude) stationarity with $80\%$ confidence, despite the sample coming from
a stationary stochastic process. As all three stationarity tests rely on
ergodicity for the first two moments, one might be tempted to think that this
could be an indication that the underlying stochastic process is not ergodic.
However, a sufficient condition for a mean zero stationary Gaussian process to
be ergodic for the first two moments is that its covariance function vanishes
as the lag increases (see Papoulis, Pillai, 2002, §13.1)

$\displaystyle\forall
t,~{}k(\tau):=\text{cov}(z_{t},z_{t+\tau})\underset{\tau\to+\infty}{\longrightarrow}0,$

and this condition is satisfied by the squared exponential kernel. The real
issue at play here is that, given a finite sample of a time series, it is
hardly possible to say with confidence that it comes from a nonstationary time
series.

Gaussianity: Similarly, with one single realisation, it is hardly possible to
determine whether the sample comes from a multivariate Gaussian, and thus the
Gaussian process assumption cannot be invalidated either. Testing whether a
real-valued vector is a draw from a multivariate Gaussian without any
assumption on its mean and covariance matrix is as hopeless as testing whether
a real-value scalar is a draw from a univariate Gaussian random variable
without any assumption on its mean or variance.

## Appendix Appendix B

Objective: In this section we prove that if a real-valued trend-stationary
Gaussian process with a differentiable mean function is mean square
differentiable and Markovian, then it has a constant covariance function (or
equivalently it is almost surely equal to its mean function plus a constant).

Proof  Let us consider a real-valued trend-stationary, mean square
differentiable and Markovian Gaussian process $(z_{t})_{t\geq 0}$. Let us
denote $(u,v)\to k(u-v)$ its covariance function. It follows that

$\forall
t,h,~{}\text{cov}\left(z_{t+h},z_{t-h}|z_{t}\right)=k(2h)-\frac{k(h)^{2}}{k(0)}=0,$
(14)

where the first equality results from standard Gaussian identities and the
second equality results from $(z_{t})_{t\geq 0}$ being Markovian. As
$(z_{t})_{t\geq 0}$ is also assumed to be mean square differentiable and to
have a differentiable mean function, the centred Gaussian process
$(z_{t}-\text{E}(z_{t}))_{t\geq 0}$ is mean square differentiable, or
equivalently $k$ is twice differentiable everywhere. Hence, $h\to
k(2h)-\frac{k(h)^{2}}{k(0)}$ is also twice differentiable at $0$ and has
second order derivative $2k^{(2)}(0)$ at $0$. It then follows from Equation
(14) that $k^{(2)}(0)=0$. As $h\to-k^{(2)}(h)$ is the covariance function of
$(z^{(1)}_{t})_{t\geq 0}$, we have:

$\displaystyle
k^{(2)}(h)^{2}=\text{cov}\left(z^{(1)}_{t},z^{(1)}_{t+h}\right)^{2}$
$\displaystyle\leq\text{var}\left(z^{(1)}_{t}\right)\text{var}\left(z^{(1)}_{t+h}\right)$
$\displaystyle\leq k^{(2)}\left(0\right)^{2}=0.$

That is, $\forall h,k^{(2)}\left(h\right)=0$, or equivalently $k$ is a linear
function of $h$. Moreover, as $(z_{t})_{t\geq 0}$ is trend-stationary $k$ is
also bounded777More precisely, the positive-definiteness of the covariance
matrix between $z_{t}$ and $z_{t+h}$ implies $\forall h,|k(h)|\leq k(0)$..
Hence, $k$ is a constant function. This means that for any times $u$ and $v$,
$z_{u}-\text{E}(z_{u})$ and $z_{v}-\text{E}(z_{v})$ have the same mean, the
same variance and correlation $1$, which implies

$\forall u,v\geq
0,~{}z_{u}-\text{E}(z_{u})\overset{\text{a.s.}}{=}z_{v}-\text{E}(z_{v}).$

## Appendix Appendix C

Objective: In this section we derive the state space representation of the
Gaussian process regression model under a trend-stationary latent GP
$(z_{t})_{t\geq 0}$, with mean function $m$, Matérn-($p+\frac{1}{2}$) kernel,
and with a Gaussian white noise $(\epsilon_{t})_{t\geq 0}$.

Derivation: If we denote $f_{t}=z_{t}-m(t)$, $(x_{t})_{t\geq 0}$ with
$x_{t}=\left(f_{t},f_{t}^{(1)},\dots,f_{t}^{(p)}\right)$ the $p$DGP
corresponding to $(f_{t})_{t\geq 0}$, and $H=(1,0,\dots,0)$, it is easy to see
that if we use $x_{t}$ as state variable, the measurement equation of the
state space model should be

$\forall t,~{}y_{t}=m(t)+H^{T}x_{t}+\epsilon_{t}.$

Denoting $K_{u,v}$ the cross-covariance matrix between $x_{u}$ and $x_{v}$,
$K_{u|v}=L_{u|v}L_{u|v}^{T}$ the auto-covariance matrix of $x_{u}$ conditional
on $x_{v}$, and $t_{0}$ the initial time, we get the initial state
distribution:

$x_{t_{0}}\sim\mathcal{N}(0,K_{t_{0},t_{0}}).$

For observations times $t_{0},\dots,t_{T}$, using Bayes’ rule and Corollary 3,
we get

$\displaystyle p\left(x_{t_{0}},\dots,x_{t_{T}}\right)$
$\displaystyle=p(x_{t_{0}})\prod_{k=1}^{T}p\left(x_{t_{k}}|x_{t_{0}:t_{k-1}}\right)$
$\displaystyle=p(x_{t_{0}})\prod_{k=1}^{T}p\left(x_{t_{k}}|x_{t_{k-1}}\right).$

Moreover, we note that

$x_{t_{k}}|x_{t_{k-1}}\sim\mathcal{N}\left(F_{t_{k}}x_{t_{k-1}},K_{t_{k}|t_{k-1}}\right)$

with $F_{t_{k}}=K_{t_{k},t_{k-1}}K_{t_{k-1},t_{k-1}}^{-1}$, and we recall that
if $X$ is a vector of $(p+1)$ i.i.d. standard normal, $M$ a deterministic
vector and $L$ a square matrix both with appropriate dimensions, then $M+LX$
is a Gaussian vector with mean $M$ and covariance matrix $LL^{T}$. It then
follows that the full state space representation reads

$\begin{cases}x_{t_{0}}\sim\mathcal{N}(0,K_{t_{0},t_{0}})\\\
x_{t_{k}}=F_{t_{k}}x_{t_{k-1}}+L_{t_{k}|t_{k-1}}\xi_{t_{k}}\\\
y_{t_{k}}=m(t_{k})+H^{T}x_{t_{k}}+\epsilon_{t_{k}}\end{cases}$

where $(\xi_{t})_{t\geq 0}$ is a $(p+1)$-dimensional standard Gaussian white
noise.

## Appendix Appendix D

Objective: In this section we prove that the $p$M-GP filter is equivalent to
Gaussian process regression under a spectral Matérn kernel (Proposition 5).

Proof  Following the notations of Definition 4 and Proposition 5, and using
the result of Appendix C, we first note that by construction the processes
$\\{\dots,({}^{i}_{c}x_{t})_{t\geq 0},({}^{i}_{s}x_{t})_{t\geq 0},\dots\\}$
are mutually independent and both $({}^{i}_{c}x_{t})_{t\geq 0}$ and
$({}^{i}_{s}x_{t})_{t\geq 0}$ are $p$DGP with mean $0$ and Matérn covariance
function $k_{\text{ma}}(\tau;k_{0i},l_{i},p+\frac{1}{2})$. Writting

${}^{i}_{*}x_{t}=\left({}^{i}_{*}z_{t},{}^{i}_{*}z_{t}^{(1)},\dots,{}^{i}_{*}z_{t}^{(p)}\right)$

where $*$ is either $c$ or $s$, it is easy to see that the observations
process of the $p$M-GP filter can be written down as

$y_{t}=z_{t}+\epsilon_{t},$

where

$z_{t}=m(t)+\sum_{i=0}^{n}{}^{i}_{c}z_{t}\cos(\omega_{i}t)+{}^{i}_{s}z_{t}\sin(\omega_{i}t).$
(15)

It is also easy to see that $(z_{t})_{t\geq 0}$ is a Gaussian process with
mean function $m$, and by mutual independence of
$\\{\dots,({}^{i}_{c}x_{t})_{t\geq 0},({}^{i}_{s}x_{t})_{t\geq 0},\dots\\}$ we
also have that

$\displaystyle\text{cov}(z_{u},z_{v})$
$\displaystyle=\sum_{i=0}^{n}\bigg{(}\text{cov}({}^{i}_{c}z_{u},{}^{i}_{c}z_{v})\cos(\omega_{i}u)\cos(\omega_{i}v)$
$\displaystyle~{}~{}~{}~{}~{}~{}+\text{cov}({}^{i}_{s}z_{u},{}^{i}_{s}z_{v})\sin(\omega_{i}u)\sin(\omega_{i}v)\bigg{)}$
$\displaystyle=\sum_{i=0}^{n}k_{\text{ma}}\left(\tau;k_{0i},l_{i},p+\frac{1}{2}\right)\cos\left(\omega_{i}(u-v)\right).$

This proves that $(z_{t})_{t\geq 0}$ has the same law as $(\hat{z}_{t})_{t\geq
0}$. Given that $(\epsilon_{t})_{t\geq 0}$ has the same law as
$(\hat{\epsilon}_{t})_{t\geq 0}$, it follows that $(y_{t})_{t\geq 0}$ has the
same law as $(\hat{y}_{t})_{t\geq 0}$, which concludes the proof.

## Appendix Appendix E

Objective: In this section we derive the solution to the forecasting problem
for the $p$M-GP filter (Equations (8, 9, 10)); in particular we provide an
iterative algorithm to compute the posterior distribution over future values
of the latent time series given historical noisy observations:
$p(z_{t}|y_{t_{k}}\dots y_{t_{0}})$ for $0<t_{0}<\dots<t_{k}<t$.

Derivation: The derivation is almost identical to the Bayesian derivation of
the Kalman filter, except for the additional trend term $m$ in the observation
equation.

We note from Equation (7) that the processes $(y_{t})_{t\geq 0}$ and
$(\bm{x}_{t})_{t\geq 0}$ are jointly Gaussian. Hence, $\bm{x}_{t_{0}}$ and
$y_{t_{0}}$ are jointly Gaussian, which implies $\bm{x}_{t_{0}}|y_{t_{0}}$ is
Gaussian too. Using Equation (7) and Proposition 5 we get:

$\displaystyle\forall
t>0,~{}\text{E}(\bm{x}_{t})=0,~{}\text{E}(y_{t})=m(t),~{}\text{E}(z_{t})=m(t)$
$\displaystyle\text{cov}\left(\bm{x}_{t},y_{t}\right)$
$\displaystyle:=\text{E}\left(\bm{x}_{t}y_{t}\right)-\text{E}\left(\bm{x}_{t}\right)\text{E}\left(y_{t}\right)$
$\displaystyle=m(t)\text{E}\left(\bm{x}_{t}\right)+\text{E}\left(\bm{x}_{t}\bm{H}_{t}^{T}\bm{x}_{t}\right)+\text{E}\left(\bm{x}_{t}\epsilon_{t}\right)$
$\displaystyle=\bm{K}_{t,t}\bm{H}_{t}$
$\displaystyle\text{cov}\left(\bm{x}_{t},\bm{x}_{t}\right):=\bm{K}_{t,t}$
$\displaystyle\text{var}(y_{t})=\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}+\sigma^{2}$
$\displaystyle\text{cov}\left(y_{t},z_{t}\right)=\text{cov}\left(z_{t},z_{t}\right)=\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}.$

Moreover, using standard Gaussian identities, we get:

$\displaystyle\text{E}\left(\bm{x}_{t}|y_{t}\right)=\frac{1}{\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}+\sigma^{2}}\bm{K}_{t,t}\bm{H}_{t}\left(y_{t}-m(t)\right)$
$\displaystyle\text{cov}\left(\bm{x}_{t}|y_{t}\right)=\bm{K}_{t,t}-\frac{1}{\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}+\sigma^{2}}\bm{K}_{t,t}\bm{H}_{t}\bm{H}_{t}^{T}\bm{K}_{t,t}^{T}$
$\displaystyle\text{E}\left(z_{t}|y_{t}\right)$
$\displaystyle=m(t)+\frac{1}{\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}+\sigma^{2}}\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}\left(y_{t}-m(t)\right)$
$\displaystyle=m(t)+\bm{H}_{t}^{T}\text{E}\left(\bm{x}_{t}|y_{t}\right)$
$\displaystyle\text{cov}\left(z_{t}|y_{t}\right)$
$\displaystyle=\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}-\frac{\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}\bm{H}_{t}^{T}\bm{K}_{t,t}^{T}\bm{H}_{t}}{\bm{H}_{t}^{T}\bm{K}_{t,t}\bm{H}_{t}+\sigma^{2}}$
$\displaystyle=\bm{H}_{t}^{T}\text{cov}\left(\bm{x}_{t}|y_{t}\right)\bm{H}_{t}.$

Hence, with

$\displaystyle\bm{m}_{t_{0}}^{-}$ $\displaystyle=0$
$\displaystyle\bm{P}_{t_{0}}^{-}$ $\displaystyle=\bm{K}_{t_{0},t_{0}}$

and

$\displaystyle\forall
t\begin{cases}v_{t}^{-}&:=\bm{H}_{t}^{T}\bm{P}_{t}^{-}\bm{H}_{t}+\sigma^{2}\\\
\bm{e}_{t}^{-}&:=y_{t}-m(t)-\bm{H}_{t}^{T}\bm{m}_{t}^{-}\\\
\bm{G}_{t}&:=\frac{1}{v_{t}^{-}}\bm{P}_{t}^{-}\bm{H}_{t}\\\
\bm{m}_{t}&:=\bm{m}_{t}^{-}+\bm{e}_{t}^{-}\bm{G}_{t}\\\
\bm{P}_{t}&:=\bm{P}_{t}^{-}-v_{t}^{-}\bm{G}_{t}\bm{G}_{t}^{T}\\\
v_{t}&:=\bm{H}_{t}^{T}\bm{P}_{t}\bm{H}_{t}\end{cases},$ (16)

we get

$\displaystyle\bm{x}_{t_{0}}|y_{t_{0}}$
$\displaystyle\sim\mathcal{N}(\bm{m}_{t_{0}},\bm{P}_{t_{0}})$ $\displaystyle
z_{t_{0}}|y_{t_{0}}$
$\displaystyle\sim\mathcal{N}(m(t)+\bm{H}_{t_{0}}^{T}\bm{m}_{t_{0}},v_{t_{0}}).$

In order to derive the remaining steps, we need the following lemma (see Kom
Samo, Roberts, 2015b, Lemma 6, and Appendix G for the proof).

###### Lemma 6

Let $X$ be a multivariate Gaussian with mean $\mu_{X}$ and covariance matrix
$\Sigma_{X}$. If conditional on $X$, $Y$ is a multivariate Gaussian with mean
$MX+A$ and covariance matrix $\Sigma_{Y}^{c}$ where $M$, $A$ and
$\Sigma_{Y}^{c}$ do not depend on $X$, then $(X,Y)$ is a jointly Gaussian
vector with mean

$\mu_{X;Y}=\begin{bmatrix}\mu_{X}\\\ M\mu_{X}+A\end{bmatrix}$

and covariance matrix

$\Sigma_{X;Y}=\begin{bmatrix}\Sigma_{X}&\Sigma_{X}M^{T}\\\
M\Sigma_{X}&\Sigma_{Y}^{c}+M\Sigma_{X}M^{T}\end{bmatrix}.$

For $k>0$ and $t>t_{k-1}$ we proceed by iteration. We assume that

$\displaystyle\begin{cases}\bm{x}_{t_{k-1}}|y_{t_{0}:t_{k-1}}&\sim\mathcal{N}(\bm{m}_{t_{k-1}},\bm{P}_{t_{k-1}})\\\
z_{t_{k-1}}|y_{t_{0}:t_{k-1}}&\sim\mathcal{N}(m(t_{k-1})+\bm{H}_{t_{k-1}}^{T}\bm{m}_{t_{k-1}},v_{t_{k-1}})\end{cases}$
(17)

using the definitions in Equations (16), and we would like to prove that

$\displaystyle\begin{cases}\bm{x}_{t}|y_{t_{0}:t_{k-1}}&\sim\mathcal{N}(\bm{m}_{t}^{-},\bm{P}_{t}^{-})\\\
y_{t}|y_{t_{0}:t_{k-1}}&\sim\mathcal{N}(m(t)+\bm{H}_{t}^{T}\bm{m}_{t}^{-},v_{t}^{-})\\\
z_{t}|y_{t_{0}:t_{k-1}}&\sim\mathcal{N}(m(t)+\bm{H}_{t}^{T}\bm{m}_{t}^{-},v_{t}^{-}-\sigma^{2}),\end{cases}$
(18)

with

$\displaystyle\begin{cases}\bm{m}_{t}^{-}&=\bm{F}_{t}\bm{m}_{t_{k-1}}\\\
\bm{P}_{t}^{-}&=\bm{F}_{t}\bm{P}_{t_{k-1}}\bm{F}_{t}^{T}+\bm{K}_{t|t_{k-1}}.\end{cases}$

To do so, we note that

$\displaystyle p\left(\bm{x}_{t},\bm{x}_{t_{k-1}}|y_{t_{0}:t_{k-1}}\right)$
$\displaystyle=p\left(\bm{x}_{t}|\bm{x}_{t_{k-1}},y_{t_{0}:t_{k-1}}\right)p\left(\bm{x}_{t_{k-1}}|y_{t_{0}:t_{k-1}}\right)$
$\displaystyle=p\left(\bm{x}_{t}|\bm{x}_{t_{k-1}}\right)p\left(\bm{x}_{t_{k-1}}|y_{t_{0}:t_{k-1}}\right),$

where the first equality results from Bayes’ rule and the second result from
the Markov property of $(\bm{x}_{t})_{t\geq 0}$. As

$\displaystyle\bm{x}_{t}|\bm{x}_{t_{k-1}},y_{t_{0}:t_{k-1}}$
$\displaystyle\sim\mathcal{N}\left(\bm{F}_{t}\bm{x}_{t_{k-1}},\bm{K}_{t|t_{k-1}}\right)$
$\displaystyle\bm{x}_{t_{k-1}}|y_{t_{0}:t_{k-1}}$
$\displaystyle\sim\mathcal{N}(\bm{m}_{t_{k-1}},\bm{P}_{t_{k-1}}),$

it follows from Lemma 6 that

$\displaystyle\bm{x}_{t}|y_{t_{0}:t_{k-1}}$
$\displaystyle\sim\mathcal{N}(\bm{F}_{t}\bm{m}_{t_{k-1}},\bm{K}_{t|t_{k-1}}+\bm{F}_{t}\bm{P}_{t_{k-1}}\bm{F}_{t}^{T})$
$\displaystyle\sim\mathcal{N}(\bm{m}_{t}^{-},\bm{P}_{t}^{-}).$

Moreover, as $(y_{t})_{t\geq 0}$ is a Gaussian process,
$y_{t}|y_{t_{0}:t_{k-1}}$ is Gaussian. What’s more, as
$y_{t}=m(t)+\bm{H}_{t}^{T}\bm{x}_{t}+\epsilon_{t}$ and $\epsilon_{t}\perp
y_{t_{0}:t_{k-1}}$, we have that:

$\displaystyle\text{E}(y_{t}|y_{t_{0}:t_{k-1}})$
$\displaystyle=m(t)+\bm{H}_{t}^{T}\text{E}\left(\bm{x}_{t}|y_{t_{0}:t_{k-1}}\right)$
$\displaystyle=m(t)+\bm{H}_{t}^{T}\bm{m}_{t}^{-}$

and

$\displaystyle\text{var}(y_{t}|y_{t_{0}:t_{k-1}})$
$\displaystyle=\bm{H}_{t}^{T}\text{cov}\left(\bm{x}_{t},\bm{x}_{t}|y_{t_{0}:t_{k-1}}\right)\bm{H}_{t}+\text{E}(\epsilon_{t}^{2})$
$\displaystyle=\bm{H}_{t}^{T}\bm{P}_{t}^{-}\bm{H}_{t}+\sigma^{2}$
$\displaystyle=v_{t}^{-}.$

The distribution $z_{t}|y_{t_{0}:t_{k-1}}$ is obtained in a similar fashion.

More generally, as the processes $(\bm{x}_{t})_{t\geq 0}$ and $(y_{t})_{t\geq
0}$ are jointly Gaussian, the random variables $\bm{x}_{t}$ and $y_{t}$ are
also jointly Gaussian conditional on $y_{t_{0}:t_{k-1}}$ and

$\displaystyle\text{cov}(\bm{x}_{t},y_{t}|y_{t_{0}:t_{k-1}})$
$\displaystyle=\text{cov}\left(\bm{x}_{t},m(t)+\bm{H}_{t}^{T}\bm{x}_{t}+\epsilon_{t}|y_{t_{0}:t_{k-1}}\right)$
$\displaystyle=\text{cov}\left(\bm{x}_{t},\bm{x}_{t}|y_{t_{0}:t_{k-1}}\right)\bm{H}_{t}+\text{E}(\bm{x}_{t}\epsilon_{t}|y_{t_{0}:t_{k-1}})$
$\displaystyle=\bm{P}_{t}^{-}\bm{H}_{t}+\text{E}(\bm{x}_{t}|y_{t_{0}:t_{k-1}})\text{E}(\epsilon_{t})$
$\displaystyle=\bm{P}_{t}^{-}\bm{H}_{t}.$

Finally, under the assumptions of Equations (17), which we recall imply
Equations (18), we would like to prove that

$\displaystyle\begin{cases}\bm{x}_{t}|y_{t_{0}:t}&\sim\mathcal{N}(\bm{m}_{t},\bm{P}_{t})\\\
z_{t}|y_{t_{0}:t}&\sim\mathcal{N}(m(t)+\bm{H}_{t}^{T}\bm{m}_{t},v_{t})\end{cases}.$
(19)

We have previously established that $\bm{x}_{t},y_{t}|y_{t_{0}:t_{k-1}}$ is
Gaussian and we have derived the corresponding mean and covariance matrix.
Noting that by definition

$\bm{x}_{t}|(y_{t_{0}:t_{k-1}},y_{t}):=\bm{x}_{t}|y_{t_{0}:t},$

it follows from standard Gaussian identities that

$\displaystyle\text{E}\left(\bm{x}_{t}|y_{t_{0}:t}\right)$
$\displaystyle=\text{E}\left(\bm{x}_{t}|y_{t_{0}:t_{k-1}}\right)+\frac{y_{t}-\text{E}\left(y_{t}|y_{t_{0}:t_{k-1}}\right)}{\text{var}\left(y_{t}|y_{t_{0}:t_{k-1}}\right)}$
$\displaystyle~{}~{}~{}~{}\times\text{cov}\left(\bm{x}_{t},y_{t}|y_{t_{0}:t_{k-1}}\right)$
$\displaystyle=\bm{m}_{t}^{-}+\frac{1}{v_{t}^{-}}\bm{P}_{t}^{-}\bm{H}_{t}\left(y_{t}-m(t)-\bm{H}_{t}^{T}\bm{m}_{t}^{-}\right)$
$\displaystyle=\bm{m}_{t}^{-}+\bm{G}_{t}\bm{e}_{t}^{-}$
$\displaystyle=\bm{m}_{t}$

and

$\displaystyle\text{cov}\left(\bm{x}_{t},\bm{x}_{t}|y_{t_{0}:t}\right)$
$\displaystyle=\text{cov}\left(\bm{x}_{t},\bm{x}_{t}|y_{t_{0}:t_{k-1}}\right)$
$\displaystyle-\frac{\text{cov}\left(\bm{x}_{t},y_{t}|y_{t_{0}:t_{k-1}}\right)\text{cov}\left(\bm{x}_{t},y_{t}|y_{t_{0}:t_{k-1}}\right)^{T}}{\text{var}\left(y_{t}|y_{t_{0}:t_{k-1}}\right)}$
$\displaystyle=\bm{P}_{t}^{-}-\frac{1}{v_{t}^{-}}\bm{P}_{t}^{-}\bm{H}_{t}\bm{H}_{t}^{T}\bm{P}_{t}^{-T}$
$\displaystyle=\bm{P}_{t}^{-}-v_{t}^{-}\bm{G}_{t}\bm{G}_{t}^{T}$
$\displaystyle=\bm{P}_{t}.$ (20)

This proves the first part of Equations (19). As for the second part, it is a
direct consequence of

$z_{t}=m(t)+\bm{H}_{t}^{T}\bm{x}_{t}.$

## Appendix Appendix F

Objective: In this section we derive the solution to the constrained
optimization problem (12).

Proof  We start by recalling the problem of interest:

$\displaystyle\begin{cases}\bm{\theta}_{t_{k}}=\underset{\bm{\theta}}{\text{argmin}}~{}||\bm{\theta}-\bm{\theta}_{t_{k-1}}||^{2}+c_{k}\xi^{2}\\\
\text{ s.t.
}\max\left(-\epsilon-\hat{\mathcal{L}}^{l}_{t_{k}}(\bm{\theta}),0\right)\leq\xi\end{cases},$

with
$\hat{\mathcal{L}}^{l}_{t_{k}}(\bm{\theta})=\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})+\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})^{T}(\bm{\theta}-\bm{\theta}_{t_{k-1}})$,
$c_{k}>0$ and $\epsilon\geq 0$.

It is easy to note that when
$\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})>-\epsilon$, the solution to
the optimization problem is $(\bm{\theta}_{t_{k-1}},0)$, so that we may focus
on the case $\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\leq-\epsilon$.

For $\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\leq-\epsilon$, the problem
can then be rewritten as the convex optimization problem:

$\displaystyle\begin{cases}\bm{\theta}_{t_{k}}=\underset{\bm{\theta}}{\text{argmin}}~{}||\bm{\theta}-\bm{\theta}_{t_{k-1}}||^{2}+c_{k}\xi^{2}\\\
\text{s.t.}-\epsilon-\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})-\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})^{T}(\bm{\theta}-\bm{\theta}_{t_{k-1}})-\xi\leq
0.\end{cases}$

As the constraints are linear and the domain of the objective is not
restricted, Slater’s condition is met and strong duality holds (see Boyd,
Vandenberghe, 2004, §5.2.3). The minimizer is therefore obtained by setting
the gradient of the Lagrangian

$\displaystyle||\bm{\theta}-\bm{\theta}_{t_{k-1}}||^{2}+c_{k}\xi^{2}+\lambda\bigg{(}-\epsilon-\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})-\xi$
$\displaystyle-\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})^{T}(\bm{\theta}-\bm{\theta}_{t_{k-1}})\bigg{)},\lambda\geq
0$ (21)

to $0$. Setting the gradient with respect to $\bm{\theta}$ to $0$ we get

$\bm{\theta}_{t_{k}}=\bm{\theta}_{t_{k-1}}+\frac{\lambda^{*}}{2}\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}}).$
(22)

Setting the derivative with respect to $\xi$ to zero, we get

$\xi^{*}=\frac{\lambda^{*}}{2c_{k}}.$ (23)

Finally, plugging Equations (22) and (23) into Equation (Appendix F), we can
rewrite the Lagrangian as

$-\left(\frac{\|\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\|^{2}}{4}+\frac{1}{4c}\right)(\lambda^{*})^{2}-\left(\epsilon+\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\right)\lambda^{*},$

which reaches its maximum at

$\lambda^{*}=-2c_{k}\frac{\epsilon+\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})}{1+c_{k}\|\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\|^{2}}.$
(24)

Using Equations (22) and (24), together with the result established for the
case $\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})>-\epsilon$, we conclude
that

$\bm{\theta}_{t_{k}}=\bm{\theta}_{t_{k-1}}+c_{k}\frac{\max\left(-\epsilon-\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}}),0\right)}{1+c_{k}\|\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}})\|^{2}}\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta}_{t_{k-1}}),$

which ends the proof.

## Appendix Appendix G

Objective: In this section we derive
$\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta})$. We assume that the trend
function $m$ is parametric and has parameters $\bm{\beta}$.

Derivation: We recall that $\mathcal{L}^{l}_{t_{k}}(\bm{\theta})$ is the
logarithm of the probability density function of a Gaussian. To ease
derivations we denote

$\bar{m}_{t_{k}}(\bm{\theta}):=m(t_{k})+\bm{H}_{t_{k}}^{T}\bm{m}_{t_{k}}^{-}\text{
and }\bar{v}_{t_{k}}(\bm{\theta}):=v_{t_{k}}^{-}$

the mean and variance of the corresponding Gaussian, so that

$\mathcal{L}^{l}_{t_{k}}(\bm{\theta})=-\frac{\log(2\pi)}{2}-\frac{\log\left(\bar{v}_{t_{k}}(\bm{\theta})\right)}{2}-\frac{\left(y_{t_{k}}-\bar{m}_{t_{k}}(\bm{\theta})\right)^{2}}{2\bar{v}_{t_{k}}(\bm{\theta})}.$
(25)

It then follows that:

$\displaystyle\nabla\mathcal{L}^{l}_{t_{k}}(\bm{\theta})$
$\displaystyle=\left(-\frac{1}{2\bar{v}_{t_{k}}(\bm{\theta})}+\frac{\left(y_{t_{k}}-\bar{m}_{t_{k}}(\bm{\theta})\right)^{2}}{2\bar{v}_{t_{k}}(\bm{\theta})^{2}}\right)\nabla\bar{v}_{t_{k}}(\bm{\theta})$
$\displaystyle+\left(\frac{y_{t_{k}}-\bar{m}_{t_{k}}(\bm{\theta})}{\bar{v}_{t_{k}}(\bm{\theta})}\right)\nabla\bar{m}_{t_{k}}(\bm{\theta}),$
(26)

so that all we need to do is derive the gradients of $\bar{m}_{t_{k}}$ and
$\bar{v}_{t_{k}}$ (using Equations (9) and (10)).

To do so, we recall that $\bm{\theta}$ is made of the parameters of $m$ that
we denote $\bm{\beta}$, $\log\sigma$ and $\\{\log k_{0i},\log
l_{i},\log\omega_{i}\\}_{i=0}^{n}$.

Derivatives with respect to $\bm{\beta}$:

$\displaystyle\frac{\partial\bar{m}_{t_{k}}(\bm{\theta})}{\partial\bm{\beta}}=\frac{\partial
m}{\partial\bm{\beta}},~{}\frac{\partial\bar{v}_{t_{k}}(\bm{\theta})}{\partial\bm{\beta}}=0.$

Derivatives with respect to $\log\sigma$:

$\displaystyle\frac{\partial\bar{m}_{t_{k}}(\bm{\theta})}{\partial\log\sigma}=0,~{}\frac{\partial\bar{v}_{t_{k}}(\bm{\theta})}{\partial\log\sigma}=\frac{\partial\bar{v}_{t_{k}}(\bm{\theta})}{\partial\sigma^{2}}\frac{\text{d}\sigma^{2}}{\text{d}\log\sigma}=2\sigma^{2}.$

Derivatives with respect to $\log\omega_{i}$:

$\displaystyle\frac{\partial\bar{m}_{t_{k}}(\bm{\theta})}{\partial\log\omega_{i}}=\frac{\partial\bm{H}_{t_{k}}^{T}}{\partial\log\omega_{i}}\bm{m}_{t_{k}}^{-}$
$\displaystyle\frac{\partial\bar{v}_{t_{k}}(\bm{\theta})}{\partial\log\omega_{i}}=\frac{\partial\bm{H}_{t_{k}}^{T}}{\partial\log\omega_{i}}\bm{P}_{t_{k}}^{-}\bm{H}_{t_{k}}+\bm{H}_{t_{k}}^{T}\bm{P}_{t_{k}}^{-}\frac{\partial\bm{H}_{t_{k}}}{\partial\log\omega_{i}}$

where $\frac{\partial\bm{H}_{t_{k}}^{T}}{\partial\log\omega_{i}}$ is identical
to $\bm{H}_{t_{k}}^{T}$ except that all terms in $\omega_{j},~{}j\neq i$ are
set to $0$, the term in $\cos(\omega_{i}t_{k})$ becomes
$-t_{k}\omega_{i}\sin(\omega_{i}t_{k})$, and the term in
$\sin(\omega_{i}t_{k})$ becomes $t_{k}\omega_{i}\cos(\omega_{i}t_{k})$.

Derivatives with respect to $\log k_{0i}$ and $\log l_{i}$ : We recall that

$\displaystyle\bar{m}_{t_{k}}(\bm{\theta})$
$\displaystyle=m(t_{k})+\bm{H}_{t_{k}}^{T}\bm{F}_{t_{k}}\bm{m}_{t_{k-1}}$
$\displaystyle\bar{v}_{t_{k}}(\bm{\theta})$
$\displaystyle=\bm{H}_{t_{k}}^{T}\bm{P}_{t_{k}}^{-}\bm{H}_{t_{k}}+\sigma^{2}$

with
$\bm{P}_{t_{k}}^{-}=\bm{F}_{t_{k}}\bm{P}_{t_{k-1}}\bm{F}_{t_{k}}^{T}+\bm{K}_{t_{k}|t_{k-1}}$.
Hence, the only terms that depend on $\log k_{0i}$ and $\log l_{i}$ are
$\bm{F}_{t_{k}}$ and $\bm{K}_{t_{k}|t_{k-1}}$, and partial derivatives of
$\bar{m}_{t_{k}}$ and $\bar{v}_{t_{k}}$ with respect to $\log k_{0i}$ and
$\log l_{i}$ are easily derived from that of $\bm{F}_{t_{k}}$ and
$\bm{K}_{t_{k}|t_{k-1}}$:

$\displaystyle\frac{\partial\bar{m}_{t_{k}}(\bm{\theta})}{\partial\log
k_{0i}}=\bm{H}_{t_{k}}^{T}\frac{\partial\bm{F}_{t_{k}}}{\partial\log
k_{0i}}\bm{m}_{t_{k-1}}$
$\displaystyle\frac{\partial\bar{m}_{t_{k}}(\bm{\theta})}{\partial\log
l_{i}}=\bm{H}_{t_{k}}^{T}\frac{\partial\bm{F}_{t_{k}}}{\partial\log
l_{i}}\bm{m}_{t_{k-1}}$
$\displaystyle\frac{\partial\bar{v}_{t_{k}}(\bm{\theta})}{\partial\log
k_{0i}}=\bm{H}_{t_{k}}^{T}\frac{\partial\bm{P}_{t_{k}}^{-}}{\partial\log
k_{0i}}\bm{H}_{t_{k}}$
$\displaystyle\frac{\partial\bar{v}_{t_{k}}(\bm{\theta})}{\partial\log
l_{i}}=\bm{H}_{t_{k}}^{T}\frac{\partial\bm{P}_{t_{k}}^{-}}{\partial\log
l_{i}}\bm{H}_{t_{k}}$

with

$\displaystyle\frac{\partial\bm{P}_{t_{k}}^{-}}{\partial\log
k_{0i}}=\frac{\partial\bm{K}_{t_{k}|t_{k-1}}}{\partial\log
k_{0i}}+\frac{\partial\bm{F}_{t_{k}}}{\partial\log
k_{0i}}\bm{P}_{t_{k-1}}\bm{F}_{t_{k}}^{T}+$
$\displaystyle+\bm{F}_{t_{k}}\bm{P}_{t_{k-1}}\frac{\partial\bm{F}_{t_{k}}^{T}}{\partial\log
k_{0i}}$ $\displaystyle\frac{\partial\bm{P}_{t_{k}}^{-}}{\partial\log
l_{i}}=\frac{\partial\bm{K}_{t_{k}|t_{k-1}}}{\partial\log
l_{i}}+\frac{\partial\bm{F}_{t_{k}}}{\partial\log
l_{i}}\bm{P}_{t_{k-1}}\bm{F}_{t_{k}}^{T}+$
$\displaystyle+\bm{F}_{t_{k}}\bm{P}_{t_{k-1}}\frac{\partial\bm{F}_{t_{k}}^{T}}{\partial\log
l_{i}}.$


* [21] Robin Cooper  “Copredication, Quantification and Frames”  In _Logical Aspects of Computational Linguistics_ , Lecture Notes in Computer Science 6736  BerlinHeidelberg: Springer, 2011, pp. 64–79  DOI: 10.1007/978-3-642-22221-4_5
  * [22] Robin Cooper  “From Perception to Communication”, Oxford Studies in Semantics and Pragmatics 16  Oxford, UK: Oxford University Press, 2023  DOI: 10.1093/oso/9780192871312.001.0001
  * [23] Robin Cooper, Simon Dobnik, Shalom Lappin and Staffan Larsson  “Probabilistic Type Theory and Natural Language Semantics”  In _Linguistic Issues in Language Technology – LiLT_ 10.4, 2015, pp. 1–43  URL: https://aclanthology.org/2015.lilt-10.4
  * [24] Kensy Cooperrider  “Foreground gesture, background gesture”  In _Gesture_ 16.2, 2017, pp. 176–202  DOI: 10.1075/gest.16.2.02coo
  * [25] Jaap M. Does and Michiel Lambalgen  “A Logic of Vision”  In _Linguistics and Philosophy_ 23.1, 2000, pp. 1–92  DOI: 10.1023/A:1005584121495
  * [26] Susan D. Duncan  “Gesture, Verb Aspect, and the Nature of Iconic Imagery in Natural Discourse”  In _Gesture_ 2.2, 2002, pp. 183–206
  * [27] Cornelia Ebert  “The Non-At-Issue Contributions of Gestures”, Talk given at the Workshop on Demonstration and Demonstratives, April 11-12 2014, Stuttgart, 2014
  * [28] Cornelia Ebert  “Semantics of Gesture”  In _Annual Review of Linguistics_ , 2024  DOI: 10.1146/annurev-linguistics-022421-063057
  * [29] Cornelia Ebert and Christian Ebert  “Gestures, demonstratives, and the attributive/referential distinction”, Handout of a talk given at Semantics and Philosophy in Europe (SPE 7), Berlin, 25
  * [30] Paul Egré and Guy Politzer  “On the negation of indicative conditionals”  In _Proceedings of the 19th Amsterdam Colloquium_ , 2013, pp. 10–18
  * [31] Stefan Engelberg  “Verben, Ereignisse und das Lexikon”, Linguistische Arbeiten 414  Niemeyer, 2000
  * [32] Maria Esipova  “Composition and projection of co-speech gestures”  In _Proceedings of the 29th Semantics and Linguistic Theory Conference_ , SALT 29, 2019, pp. 117–137
  * [33] Charles J. Fillmore  “The case for case”  In _Universals of linguistic theory_ New York, NY: Holt, Rinehart & Winston, 1968, pp. 1–88
  * [34] Charles J. Fillmore and Colin Baker  “A Frames Approach to Semantic Analysis”  In _The Oxford Handbook of Linguistic Analysis_ Oxford: Oxford University Press, 2010, pp. 313–340
  * [35] Charles J. Fillmore, Collin F. Baker and Hiroaki Sato  “FrameNet as a “Net””  In _Proceedings of the Fourth International Conference on Language Resources and Evaluation_ , LREC’2004, 2004, pp. 1091–1094
  * [36] Jonathan Ginzburg  “The Interactive Stance: Meaning for Conversation”  Oxford, UK: Oxford University Press, 2012
  * [37] Jonathan Ginzburg and Robin Cooper  “Clarification, Ellipsis, and the Nature of Contextual Updates in Dialogue”  In _Linguistics and Philosophy_ 27.3, 2004, pp. 297–365  DOI: 10.1023/B:LING.0000023369.19306.90
  * [38] Jonathan Ginzburg and Andy Lücking  “I thought pointing is rude: A dialogue-semantic analysis of pointing at the addressee” Special Session: Gestures and Natural Language Semantics  In _Proceedings of Sinn und Bedeutung 25_ , SuB 25, 2021, pp. 276–291  DOI: 10.18148/sub/2021.v25i0.937
  * [39] Jonathan Ginzburg and Andy Lücking  “Requesting clarifications with speech and gestures”  In _Proceedings of the 1st Workshop on Multimodal Semantic Representations_ , MMSR  Groningen, Netherlands (Online): Association for Computational Linguistics, 2021, pp. 21–31  URL: https://aclanthology.org/2021.mmsr-1.3
  * [40] Jonathan Ginzburg et al.  “Characterizing the Response Space of Questions: Data and Theory”  In _Dialogue & Discourse_ 13.2, 2022, pp. 79–132  DOI: org/10.5210/dad.2022.203
  * [41] Gianluca Giorgolo  “Space and Time in Our Hands”, 2010
  * [42] Nelson Goodman  “Languages of Art. An Approach to a Theory of Symbols”  Idianapolis: Hackett Publishing Company, Inc., 1976
  * [43] Gabriel Greenberg  “The Semiotic Spectrum”, 2011
  * [44] Gabriel Greenberg  “Semantics of Pictorial Space”  In _Review of Philosophy and Psychology_ 12, 2021, pp. 847–887  DOI: 10.1007/s13164-020-00513-6
  * [45] Marianne Gullberg and Sotaro Kita  “Attention to Speech-Accompanying Gestures: Eye Movements and Information Uptake”  In _Journal of Nonverbal Behavior_ 33.4, 2009, pp. 251–277  DOI: 10.1007/s10919-009-0073-2
  * [46] Uri Hadar  “Two Types of Gesture and Their Role in Speech Production”  In _Journal of Language and Social Psychology_ 8.3-4, 1989, pp. 221–228  DOI: 10.1177/0261927X8983004
  * [47] Uri Hadar  “Co-verbal gestures: Between communication and speech production”  In _Body – Language – Communication. An International Handbook on Multimodality in Human Interaction_ 1.38, Handbooks of Linguistics and Communication Science  BerlinBoston: De Gruyter Mouton, 2013, pp. 804–821  DOI: https://doi.org/10.1515/9783110261318.804
  * [48] Uri Hadar and Robert K. Krauss  “Iconic gestures: The grammatical categories of lexical affiliates”  In _Journal of Neurolinguistics_ 12.1, 1999, pp. 1–12  DOI: 10.1016/S0911-6044(99)00001-9
  * [49] Florian Hahn and Hannes Rieser  “Explaining Speech Gesture Alignment in MM Dialogue Using Gesture Typology”  In _Proceedings of the 14th Workshop on the Semantics and Pragmatics of Dialogue: Aspects of Semantics and Pragmatics of Dialogue_ , SemDial 2010, 2010, pp. 99–111  Polish Society for Cognitive Science
  * [50] Stevan Harnad  “The Symbol Grounding Problem”  In _Physica D_ 42, 1990, pp. 335–346
  * [51] Patrick G.T. Healey, Nicola Plant, Christine Howes and Mary Lavelle  “When Words Fail: Collaborative Gestures During Clarification Dialogues”  In _Turn-Taking and Coordination in Human-Machine Interaction: Papers from the 2015 AAAI Spring Symposium_ , 2015, pp. 23–29
  * [52] Irene Heim and Angelika Kratzer  “Semantics in Generative Grammar”, Blackwell Textbooks in Linguistics 13  Malden, MAOxford, UK: Blackwell, 1998
  * [53] Judith Holler and Stephen C. Levinson  “Multimodal Language Processing in Human Communication” Opinion  In _Trends in Cognitive Sciences_ 23.8, 2019, pp. 639–652  DOI: 10.1016/j.tics.2019.05.006
  * [54] Yufang Hou, Katja Markert and Michael Strube  “Unrestricted Bridging Resolution”  In _Computational Linguistics_ 44.2 Cambridge, MA: MIT Press, 2018, pp. 237–284  DOI: 10.1162/COLI_a_00315
  * [55] Julian Hough, Iwan Kok, David Schlangen and Stefan Kopp  “Timing and Grounding in Motor Skill Coaching Interaction: Consequences for the Information State”  In _Proceedings of SEMDIAL 2015_ , goDIAL 2015, 2015, pp. 86–94
  * [56] Julie Hunter  “Relating gesture to speech: Reflections on the role of conditional presuppositions”  In _Linguistics and Philosophy_ 42, 2019, pp. 317–332  DOI: 10.1007/s10988-018-9244-0
  * [57] Matthias Irmer  “Inferring Implicatures and Discourse Relations from Frame Information” Special Issue: Implicature and Discourse Structure  In _Lingua_ 132, 2013, pp. 29–50  DOI: 10.1016/j.lingua.2012.11.009
  * [58] Ray Jackendoff  “On Beyond Zebra: The Relation of Linguistic and Visual Information”  In _Cognition_ 26, 1987, pp. 89–114
  * [59] Ray Jackendoff  “Parts and Boundaries”  In _Cognition_ 41, 1991, pp. 9–45
  * [60] Gunnar Johansson  “Visual Perception of Biological Motion and a Model for its Analysis”  In _Perception & Psychophysics_ 14.2, 1973, pp. 201–211  DOI: doi.org/10.3758/BF03212378
  * [61] Gunnar Johansson  “Spatio-Temporal Differentiation and Integration in Visual Motion Perception. An Experimental and Theoretical Analysis of Calculus-Like Functions in Visual Data Processing”  In _Psychol. Res._ 38, 1976, pp. 379–393
  * [62] Gunnar Johansson, Claes Hofsten and Gunnar Jansson  “Event Perception”  In _Annual Review of Psychology_ 31, 1980, pp. 27–63
  * [63] Hans Kamp and Antje Rossdeutscher  “DRS-Construction and Lexically Driven Inference”  In _Theoretical Linguistics_ 20.2-3, 1994, pp. 165–235  DOI: 10.1515/thli.1994.20.2-3.165
  * [64] Hans Kamp and Antje Rossdeutscher  “Remarks on Lexical Structure and DRS Construction”  In _Theoretical Linguistics_ 20.2-3, 1994, pp. 97–164  DOI: 10.1515/thli.1994.20.2-3.97
  * [65] Adam Kendon  “Gesture: Visible Action as Utterance”  Cambridge, MA: Cambridge University Press, 2004
  * [66] Christopher Kennedy  “Vagueness and grammar: The semantics of relative and absolute gradable adjectives”  In _Linguistics and Philosophy_ 30.1, 2007, pp. 1–45  DOI: 10.1007/s10988-006-9008-0
  * [67] Casey Kennington and David Schlangen  “Simple Learning and Compositional Application of Perceptually Grounded Word Meanings for Incremental Reference Resolution”  In _Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing_ , 2015, pp. 292–301
  * [68] Sotaro Kita and Aslı Özyürek  “What Does Cross-Linguistic Variation in Semantic Coordination of Speech and Gesture Reveal?: Evidence for an Interface Representation of Spatial Thinking and Speaking”  In _Journal of Memory and Language_ 48.1, 2003, pp. 16–32  DOI: 10.1016/S0749-596X(02)00505-3
  * [69] Stefan Kopp, Paul Tepper and Justine Cassell  “Towards Integrated Microplanning of Language and Iconic Gesture for Multimodal Output”  In _Proceedings of the 6th International Conference on Multimodal Interfaces_ , ICMI ’04  State College, PA, USA: ACM, 2004, pp. 97–104  DOI: 10.1145/1027933.1027952
  * [70] Alfred Kranstedt et al.  “Deictic Object Reference in Task-oriented Dialogue”  In _Situated Communication_ Berlin: Mouton de Gruyter, 2006, pp. 155–207
  * [71] Alfred Kranstedt et al.  “Measuring and Reconstructing Pointing in Visual Contexts”  In _Proceedings of the 10th Workshop on the Semantics and Pragmatics of Dialogue_ , Brandial’06  Potsdam: Universitätsverlag Potsdam, 2006, pp. 82–89
  * [72] Robert M. Krauss, Yihsiu Chen and Rebecca F. Gottesmann  “Lexical Gestures and Lexical Access: A Process Model”  In _Language and Gesture_ Cambridge, UK: Cambridge University Press, 2000, pp. 261–283
  * [73] Manfred Krifka  “The Origins of Telicity”  In _Events and Grammar_ , Studies in Linguistics and Philosophy 70  Dordrecht: Springer, 1998  DOI: org/10.1007/978-94-011-3969-4_9
  * [74] Manfred Krifka  “Some remarks on event structure, conceptual spaces and the semantics of verbs”  In _Theoretical Linguistics_ 38.3-4, 2012, pp. 223–236  DOI: 10.1515/tl-2012-0014
  * [75] Schuyler Laparle  “Tracking Discourse Topics in Co-speech Gesture”  In _Digital Human Modeling and Applications in Health, Safety, Ergonomics and Risk Management. Human Body, Motion and Behavior_ , Lecture Notes in Computer Science 12777  Cham: Springer International Publishing, 2021, pp. 233–249
  * [76] Shalom Lappin  “Curry Typing, Polymorphism, and Fine-Grained Intensionality”  In _The Handbook of Contemporary Semantic Theory_ John Wiley & Sons, Ltd, 2015, pp. 408–428  DOI: 10.1002/9781118882139.ch13
  * [77] Staffan Larsson  “Formal Semantics for Perceptual Classification”  In _Journal of Logic and Computation_ 25.2, 2015, pp. 335–369  DOI: 10.1093/logcom/ext059
  * [78] Staffan Larsson  “Extensions Are Indeterminate if Intensions Are Classifiers”  In _Proceedings of the 24th Workshop on the Semantics and Pragmatics of Dialogue – Full Papers_ Virtually at Brandeis, Waltham, New Jersey: SEMDIAL, 2020  URL: http://semdial.org/anthology/Z20-Larsson_semdial_0019.pdf
  * [79] Alex Lascarides and Matthew Stone  “A Formal Semantic Analysis of Gesture”  In _Journal of Semantics_ 26.4, 2009, pp. 393–449  DOI: 10.1093/jos/ffp004
  * [80] David Lewis  “Convention: A Philosophical Study”  Cambridge, MA: Harvard University Press, 1969
  * [81] David Lewis  “General Semantics”  In _Synthese_ 22.1/2, Semantics of Natural Language II, 1970, pp. 18–67
  * [82] Daniel Loehr  “Apects of Rhythm in Gesture in Speech”  In _Gesture_ 7.2, 2007, pp. 179–214
  * [83] Andy Lücking  “Ikonische Gesten. Grundzüge einer linguistischen Theorie”  BerlinBoston: De Gruyter, 2013
  * [84] Andy Lücking  “Multimodal Propositions? From Semiotic to Semantic Considerations in the Case of Gestural Deictics”  In _Poster Abstracts of the Proceedings of the 17th Workshop on the Semantics and Pragmatics of Dialogue_ , SemDial 2013, 2013, pp. 221–223
  * [85] Andy Lücking  “Modeling Co-verbal Gesture Perception in Type Theory with Records”  In _Proceedings of the 2016 Federated Conference on Computer Science and Information Systems_ 8, Annals of Computer Science and Information Systems  IEEE, 2016, pp. 383–392  DOI: 10.15439/2016F83
  * [86] Andy Lücking  “Witness-loaded and Witness-free Demonstratives”  In _Atypical Demonstratives. Syntax, Semantics and Pragmatics_ , Linguistische Arbeiten 568  BerlinBoston: De Gruyter, 2018, pp. 255–284
  * [87] Andy Lücking  “Gesture”  In _Head Driven Phrase Structure Grammar: The handbook_ , Empirically Oriented Theoretical Morphology and Syntax 9  Berlin: Language Science Press, 2021, pp. 1201–1250  DOI: 10.5281/zenodo.5543318
  * [88] Andy Lücking  “Aspects of Multimodal Communication”, 2022
  * [89] Andy Lücking and Jonathan Ginzburg  “Referential transparency as the proper treatment of quantification”  In _Semantics and Pragmatics_ 15.4, 2022, pp. 1–58  DOI: 10.3765/sp.15.4
  * [90] Andy Lücking and Jonathan Ginzburg  “Leading voices: Dialogue semantics, cognitive science, and the polyphonic structure of multimodal interaction”  In _Language and Cognition_ 15.1, 2023, pp. 148–172  DOI: 10.1017/langcog.2022.30
  * [91] Andy Lücking, Thies Pfeiffer and Hannes Rieser  “Pointing and Reference Reconsidered”  In _Journal of Pragmatics_ 77, 2015, pp. 56–79  DOI: 10.1016/j.pragma.2014.12.013
  * [92] Andy Lücking et al.  “The Bielefeld Speech and Gesture Alignment Corpus (SaGA)”  In _Multimodal Corpora: Advances in Capturing, Coding and Analyzing Multimodality_ , LREC 2010, 2010, pp. 92–98  7th International Conference for Language ResourcesEvaluation  DOI: 10.13140/2.1.4216.1922
  * [93] Andy Lücking et al.  “Data-based Analysis of Speech and Gesture: The Bielefeld Speech and Gesture Alignment Corpus (SaGA) and its Applications”  In _Journal on Multimodal User Interfaces_ 7.1-2, 2013, pp. 5–18  DOI: 10.1007/s12193-012-0106-8
  * [94] David Marr  “Vision”  The MIT Press, 2010  DOI: 10.7551/mitpress/9780262514620.001.0001
  * [95] David Marr and Lucia Vaina  “Representation and recognition of the movements of shapes”  In _Proceedings of the Royal Society of London. Series B. Biological Sciences_ 214.1197, 1982, pp. 501–524  DOI: 10.1098/rspb.1982.0024
  * [96] David McNeill  “So You Think Gestures are Nonverbal?”  In _Psychological Review_ 92.3, 1985, pp. 350–371
  * [97] David McNeill  “Hand and Mind – What Gestures Reveal about Thought”  Chicago: Chicago University Press, 1992
  * [98] David McNeill  “Gesture and Thought”  Chicago: University of Chicago Press, 2005
  * [99] David McNeill, Justine Cassell and Elena T. Levy  “Abstract deixis”  In _Semiotica_ 95.1-2, 1993, pp. 5–19  DOI: 10.1515/semi.1993.95.1-2.5
  * [100] David McNeill and Susan D. Duncan  “Growth points in thinking-for-speaking”  In _Language and Gesture_ Cambridge, MA: Cambridge University Press, 2000, pp. 141–161
  * [101] Alexander Mehler and Andy Lücking  “Pathways of Alignment between Gesture and Speech: Assessing Information Transmission in Multimodal Ensembles”  In _Proceedings of the International Workshop on Formal and Computational Approaches to Multimodal Communication under the auspices of ESSLLI 2012, Opole, Poland, 6-10 August_ , 2012
  * [102] Vedran Miletić et al.  “The beamer package. A LaTeX class for producing presentations and slides”, 2024  URL: https://ctan.org/pkg/beamer
  * [103] Irene Mittelberg and Linda B. Waugh  “Gestures and metonymy”  In _Body – Language – Communication. An International Handbook on Multimodality in Human Interaction_ 2.38, Handbooks of Linguistics and Communication Science  BerlinBoston: De Gruyter Mouton, 2014, pp. 1747–1766  DOI: 10.1515/9783110302028.1747
  * [104] Marcin Morzycki  “Degree modification of gradable nouns: size adjectives and adnominal degree morphemes”  In _Natural Language Semantics_ 17, 2009, pp. 175–203  DOI: 10.1007/s11050-009-9045-7
  * [105] Cornelia Müller  “Gestural modes of representation as techniques of depiction”  In _Body – Language – Communication. An International Handbook on Multimodality in Human Interaction_ 2.38, Handbooks of Linguistics and Communication Science  BerlinBoston: De Gruyter Mouton, 2014, pp. 1687–1702
  * [106] Kevin Mulligan, Peter Simons and Barry Smith  “Truth-Makers”  In _Philosophy and Phenomenological Research_ 44, 1984, pp. 287–321
  * [107] Reinhard A. Muskens  “Combining Montague Semantics and Discourse Representation”  In _Linguistics and Philosophy_ 19, 1996, pp. 143–186
  * [108] Cornelia Müller  “How recurrent gestures mean: Conventionalized contexts-of-use and embodied motivation”  In _Gesture_ 16.2 John Benjamins, 2017, pp. 277–304  DOI: 10.1075/gest.16.2.05mul
  * [109] John O’Keefe  “The Spatial Prepositions in English, Vector Grammar, and the Cognitive Map Theory”  In _Language and Space_ Cambridge, MA: MIT Press, 1996, pp. 277–316
  * [110] John O’Keefe  “Vectors Grammar, Places, and the Functional Role of the Spatial Prepositions in English”  In _Representing Direction in Language and Space_ , Explorations in Language and Space 1  Oxford, NY: Oxford University Press, 2003, pp. 69–85
  * [111] Aslı Özyürek  “Hearing and Seeing Meaning in Speech and Gesture: Insights from Brain and Behaviour”  In _Philosophical Transactions of the Royal Society B: Biological Sciences_ 369.1651, 2014  DOI: 10.1098/rstb.2013.0296
  * [112] Allan Paivio  “Mental Representations: A Dual Coding Approach”  New York: Oxford University Press, 1986
  * [113] Terence Parsons  “Events in the Semantics of English”, Current Studies in Linguistics Series 19  Cambridge: MIT Press, 1990
  * [114] Barbara H. Partee  “Possible Worlds Semantics and Linguistic Theory”  In _The Monist_ 60.3 Oxford University Press, 1977, pp. 303–326
  * [115] Miriam R. L. Petruck and Ellen K Dodge  “MetaNet: Repository, Identification System, and Applications”  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics: Tutorial Abstracts_ Berlin, Germany: Association for Computational Linguistics, 2016  URL: https://aclanthology.org/P16-5008
  * [116] Massimo Poesio et al.  “Computational Models of Anaphora”  In _Annual Review of Linguistics_ 9.1, 2023, pp. 561–587  DOI: 10.1146/annurev-linguistics-031120-111653
  * [117] Nicholas Rescher  “How Many Possible Worlds Are There?”  In _Philosophy and Phenomenological Research_ 59.2, 1999, pp. 403–420
  * [118] Hannes Rieser  “Pointing in Dialogue”  In _Proceedings of the Eighth Workshop on the Semantics and Pragmatics of Dialogue_ , Catalog ’04, 2004, pp. 93–100
  * [119] Hannes Rieser  “Aligned Iconic Gesture in Different Strata of MM Route-Description”  In _LonDial 2008: The 12th Workshop on the Semantics and Pragmatics of Dialogue (SEMDIAL)_ , 2008, pp. 167–174
  * [120] Hannes Rieser  “How to Disagree on a Church Window’s Shape Using Gesture”  In _Dimensionen der Analyse von Texten und Diskursen. Festschrift für János Sándor Petőfi_ Berlin: LIT Verlag, 2011, pp. 231–246
  * [121] Hannes Rieser  “When Hands Talk to Mouth. Gesture and Speech as Autonomous Communicating Processes”  In _Proceedings of the 19th Workshop on the Semantics and Pragmatics of Dialogue_ , SEMDIAL 2015: goDIAL, 2015, pp. 122–130
  * [122] Hannes Rieser and Massimo Poesio  “Interactive Gestures in Dialogue: a PTT Model”  In _Proceedings of the 10th Annual Meeting of the Special Interest Group in Discourse and Dialogue_ , SIGDIAL 2009, 2009, pp. 87–96
  * [123] Jan Peter Ruiter  “The Production of Gesture and Speech”  In _Language and Gesture_ Cambridge, UK: Cambridge University Press, 2000, pp. 284–311
  * [124] Emanuel A. Schegloff  “On some Gestures’ Relation to Talk”  In _Structures of Social Action. Studies in Conversational Analysis_ , Studies in Emotion and Social Interaction  Cambridge, MA: Cambridge University Press, 1984, pp. 266–296
  * [125] Philippe Schlenker  “Gesture projection and cosuppositions”  In _Linguistics and Philosophy_ 41.3, 2018, pp. 295–365  DOI: 10.1007/s10988-017-9225-8
  * [126] Philippe Schlenker  “Gestural semantics”  In _Natural Language & Linguistic Theory_ 37, 2019, pp. 735–784
  * [127] Philippe Schlenker and Emmanuel Chemla  “Gestural agreement”  In _Natural Language & Linguistic Theory_ 36, 2018, pp. 587–625  DOI: 10.1007/s11049-017-9378-8
  * [128] Atsushi Shimojima  “Semantic Properties of Diagrams and Their Cognitive Potentials”, Studies in the Theory and Applications of Diagrams  Stanford, CA: CSLI Publications, 2015
  * [129] Timo Sowa  “Understanding Coverbal Iconic Gestures in Shape Descriptions”  Berlin: Akademische Verlagsgesellschaft, 2006
  * [130] Luc Steels and Tony Belpaeme  “Coordinating perceptually grounded categories through language: A case study for colour”  In _Behavioral and Brain Sciences_ 28.4, 2005, pp. 469–529
  * [131] Jürgen Streeck  “Depicting by Gesture”  In _Gesture_ 8.3, 2008, pp. 285–301  DOI: 10.1075/gest.8.3.02str
  * [132] Leonard Talmy  “Force dynamics in language and cognition”  In _Cognitive Science_ 12.1, 1988, pp. 49–100  DOI: 10.1016/0364-0213(88)90008-0
  * [133] Shimon Ullman  “High-Level Vision”, A Bradford Book  Cambridge, MA: MIT Press, 1996
  * [134] Massimo Warglien, Peter Gärdenfors and Matthijs Westera  “Event structure, conceptual spaces and the semantics of verbs”  In _Theoretical Linguistics_ 38.3-4, 2012, pp. 159–193  DOI: 10.1515/tl-2012-0010
  * [135] Matthias Weisgerber  “Decomposing Path Shapes: About an Interplay of Manner of Motion and ‘The Path”’  In _Proceedings of Sinn und Bedeutung_ , SuB 10, 2006, pp. 405–419  DOI: 10.18148/sub/2006.v10i2.741
  * [136] Henk Zeevat  “A Compositional Approach to Discourse Representation Theory”  In _Linguistics and Philosophy_ 12.1, 1989, pp. 95–131
  * [137] Thomas Ede Zimmermann  “Model-Theoretic Semantics”  In _Semantics: An International Handbook of Natural Language Meaning_ 1.33, Handbücher zur Sprach- und Kommunikationswissenschaft / Handbooks of Linguistics and Communication Science (HSK)  BerlinBoston: De Gruyter Mouton, 2011, pp. 762–802  DOI: 10.1515/9783110226614.946
  * [138] Thomas Ede Zimmermann  “Intensions, Types and Models”  In _The Semantics of Plurals, Focus, Degrees, and Times: Essays in Honor of Roger Schwarzschild_ Cham: Springer International Publishing, 2019, pp. 41–56  DOI: 10.1007/978-3-030-04438-1_3
  * [139] Thomas Ede Zimmermann  “Extensions in Compositional Semantics”  In _Language, Logic, and Computation_ Cham: Springer International Publishing, 2022, pp. 148–172
  * [140] Joost Zwarts  “Vectors as Relative Positions: A Compositional Semantics of Modified PPs”  In _Journal of Semantics_ 14.1, 1997, pp. 57–86
  * [141] Joost Zwarts  “Vectors Across Spatial Domains: From Place to Size, Orientation, Shape, and Parts”  In _Representing Direction in Language and Space_ , Explorations in Language and Space 1  Oxford, NY: Oxford University Press, 2003, pp. 39–68
  * [142] Joost Zwarts  “Spatial Semantics: Modeling the Meaning of Prepositions” e12241 LNCO-0444.R2  In _Language and Linguistics Compass_ 11.5, 2017, pp. e12241–n/a  DOI: 10.1111/lnc3.12241
  * [143] Joost Zwarts and Peter Gärdenfors  “Locative and Directional Prepositions in Conceptual Spaces: The Role of Polar Convexity”  In _Journal of Logic, Language and Information_ 25.1, 2016, pp. 109–138
  * [144] Joost Zwarts and Yoad Winter  “Vector Space Semantics: A Model-Theoretic Analysis of Locative Prepositions”  In _Journal of Logic, Language, and Information_ 9.2, 2000, pp. 169–211