# Topological magnons in Kitaev magnets with finite Dzyaloshinskii-Moriya
interaction at high field

Kangkang Li<EMAIL_ADDRESS>Department of Physics, Zhejiang Normal
University, Jinhua, 321004, China

###### Abstract

There have been intensive studies on Kitaev materials for the sake of the
realization of exotic states such as quantum spin liquid and topological
orders. In realistic materials, the Kitaev interaction may coexist with the
Dzyaloshinskii-Moriya interaction (DMI), and it is a challenge to distinguish
their magnitudes separately. Here, we study the topological magnon excitations
and related thermal Hall conductivity of kagome magnet exhibiting Heisenberg,
Kitaev and DM interactions exposed to a magnetic field. In a strong magnetic
field perpendicular to the plane of the lattice ([111] direction) that bring
the system into the fully polarized paramagnetic phase, we find that the
magnon bands carry nontrivial Chern numbers in the full region of the phase
diagram. Furthermore, there are phase transitions relate two topological
phases with opposite Chern numbers, which lead to the sign changes of the
thermal Hall conductivity. In the phase with negative thermal conductivity,
the Kitaev interaction is relatively large and the width of the phase
increases with the strength of DMI. Hence the study here will contribute to
the understanding of related compounds.

##  I Introduction

In recent years, quantum materials with bond-dependent anisotropic Kitaev spin
interactions have been the subject of much experimental and theoretical
studies, because there are frustrations from the competing exchange couplings
coexisting with the geometric frustration of the underlying lattices, which
may lead to unusual magnetic orders as well as gapped and gapless spin liquids
with fractional excitationsRousochatzakis  _et al._ (2015); Lee  _et al._
(2014); Nasu  _et al._ (2014); Takayama  _et al._ (2015); Jahromi  _et al._
(2016, 2018); Kishimoto  _et al._ (2018); Kimchi and Vishwanath (2014).
However, the Kitaev interaction is often accompanied by the Dzyaloshinskii-
Moriya interaction (DMI)Dzyaloshinsky (1958); Moriya (1960) in realistic
materials, and it is a challenge to distinguish their magnitudes separatelyGao
_et al._ (2019); Zhang  _et al._ (2021).

The Kitaev interaction also shows its exotic properties in magnonic side, for
instance, it can realize topological magnon bands in various lattices such as
honeycomb latticeJoshi (2018); Chern  _et al._ (2021); McClarty  _et al._
(2018); Zhang  _et al._ (2021) and kagome latticeLi (2022). However, the
effect of the Kitaev interaction can be generally similar to that of the DMI,
which can also induce topological magnons and has been already studied in lots
of experimental and theoretical worksKatsura  _et al._ (2010); Onose  _et al._
(2010); Hirschberger  _et al._ (2015); Chisnell  _et al._ (2015); Owerre
(2016); Zhang  _et al._ (2013); Li  _et al._ (2017); Li and Hu (2017); Cao
_et al._ (2015). In some cases, the spin wave spectrum of the Heisenberg-
Kitaev model can even reduce to that of the Heisenberg model with DMI, and
then hosts the same topological propertyMcClarty  _et al._ (2018). Therefore,
the same magnon bands and topological property can be generated either by
Kitaev interaction, DMI, or their combination. Fortunately, it was shown that
one can distinguish whether a system is Kitaev interaction or DMI dominated by
further investigating the magnonic transport properties, for example, in
honeycomb Kitaev magnetGao  _et al._ (2019); Zhang  _et al._ (2021).

Figure 1: (a) The structure of the kagome lattice. There are three spins
reside in a primitive cell, which are denoted by red, green and blue sites.
Red, green and blue bonds between NN sites $(i,j)$ carry three distinct Kitaev
couplings $S_{i}^{x}S_{j}^{x}$, $S_{i}^{y}S_{j}^{y}$ and $S_{i}^{x}S_{j}^{x}$,
respectively. Meanwhile, there are isotropic Heisenberg couplings in all the
NN bonds. The arrows in a triangle denote the coupling directions whose DM
vectors are in the [111] direction. The field is also applied in the [111]
direction. The kagome lattice sits on the (111) plane, and we define a new 2D
frame $x^{\prime}y^{\prime}$ on the kagome lattice with basis vectors
$\bm{a}_{1}=(1,0)$ and $\bm{a}_{2}=(1,\sqrt{3})/2$. (b) The first Brillouin
zone of the kagome lattice. Figure 2: (a) The topological phase diagram of our
model as a function of $\psi$ and $D$. The topological phases are
characterized by the Chern numbers of magnon bands, and the phase diagram is
filled by two phases $(1,0,-1)$ and $(-1,0,1)$. The three red points on the
dashed line with $D=0.7$ reside in phase $(-1,0,1)$ with $\psi=1.57\pi$, at
phase boundary with $\psi=1.726\pi$ and in phase $(1,0,-1)$ with $\psi=2\pi$,
respectively. (b) The Berry curvatures of the first and third magnon bands for
$\psi=1.57\pi$ (top panel) and $\psi=2\pi$ (bottom panel) with $D=0.7$.

In this work, we study the topological magnon excitations and related thermal
Hall conductivity of kagome magnet exhibiting Heisenberg, Kitaev and DM
interactions exposed to a magnetic field. It is well known that the magnon
bands of kagome ferromagnet with DMI carry nonzero Chern numbersChisnell  _et
al._ (2015). On the other hand, the kagome magnet with Kitaev interaction,
whose ground state is canted ferromagnetic orderMorita  _et al._ (2018); Yang
_et al._ (2020), also hosts topological magnon excitationsLi (2022). Hence,
one may wonder what would happen when there are DMI and Kitaev interaction in
kagome magnet simultaneously. Since the Kitaev interaction does not support
ferromagnetic orderMorita  _et al._ (2018); Yang  _et al._ (2020), we apply a
magnetic field perpendicular to the plane of the lattice ([111] direction),
which is strong enough to bring the system into fully polarized paramagnetic
phaseMcClarty  _et al._ (2018). We find that the magnon bands carry nontrivial
Chern numbers in the full region of the phase diagram. Furthermore, there are
phase transitions relate two topological phases with opposite Chern numbers,
which lead to the sign changes of the thermal Hall conductivity. In the phase
with negative thermal conductivity, the Kitaev interaction is relatively large
and the width of the phase increases with the strength of DMI. Thus the study
here will contribute to the understanding of related compounds.

This paper is organized as follows. In section II, we introduce the model and
related methods. In section III, we discuss the topological magnons. In
section IV, we present the thermal Hall conductivity with sign change.
Finally, a summary is given in section V.

##  II Model and method

We consider interacting spins reside on the kagome lattice as shown in Fig.
1(a). The model is described by the spin Hamiltonian

$\displaystyle H=$ $\displaystyle J\sum_{\langle ij\rangle}{\bf
S}_{i}\cdot{\bf S}_{j}+K\sum_{\langle
ij\rangle}S_{i}^{\gamma_{ij}}S_{j}^{\gamma_{ij}}$
$\displaystyle+D\sum_{\langle
ij\rangle}\nu_{ij}\hat{{\bm{z}}^{\prime}}\cdot({\bf S}_{i}\times{\bf
S}_{j})-h\sum_{i}\hat{{\bm{z}}^{\prime}}\cdot{\bf S}_{i},$ (1)

where ${\bf S}_{i}$ and ${\bf S}_{j}$ are spin $S=1/2$ spins reside on the
nearest-neighbor (NN) lattice sites, and $J$ and $K$ denote the Heisenberg and
Kitaev exchange couplings, respectively. The Cartesian components
$\gamma_{ij}$ equals $x$, $y$ or $z$, depending on the bond type as shown in
Fig. 1(a). We parameterize the Heisenberg and Kitaev terms by
$J=\cos\psi,\quad K=\sin\psi,\quad\psi\in[0,2\pi)$, with the energy unit
$J^{2}+K^{2}=1$. The NN DMI is represented by the third term, where
$\hat{{\bm{z}}^{\prime}}$ is the unit vector along [111] direction and
$\nu_{ij}=\pm 1$ correspond to the anticlockwise and clockwise directions of
the NN couplings in a triangle, respectively. In kagome magnet the Kitaev
interaction does not support ferromagnetic orderMorita  _et al._ (2018); Yang
_et al._ (2020), hence we add the final term of Zeeman coupling to a applied
magnetic field along the [111] direction, which is strong enough to bring the
system into fully polarized paramagnetic phase. We study a reasonable
parameter range of $D\in[0,1]$, and set $h=10$ (in unit of $S$) all through
the work to ensure the stability of the system.

Further, the Holstein-Primakoff (HP) transformationHolstein and Primakoff
(1940) is employed to rewrite the Hamiltonian in terms of magnon creation and
annihilation operators $a_{i}^{\dagger}$ and $a_{i}$. Under the linear spin
wave approximation, we keep only the quadratic terms of the magnon operators.
Then the Fourier transformation is performed to rewrite the Hamiltonian in
momentum space with basis
$\Psi_{\bm{k}^{\prime}}^{\dagger}=(a_{1\bm{k}^{\prime}}^{\dagger},a_{2\bm{k}^{\prime}}^{\dagger},a_{3\bm{k}^{\prime}}^{\dagger},a_{1-\bm{k}^{\prime}},a_{2-\bm{k}^{\prime}},a_{3-\bm{k}^{\prime}})$,
and we get the magnon Hamiltonian matrix $h(\bm{k}^{\prime})$ (see Appendix
for details). Finally, the eigenvalues and eigenvectors are obtained by
diagonalizing the dynamic matrix $I_{-}h(\bm{k}^{\prime})$ with
$I_{-}=[(I,0),(0,-I)]$, $I$ as the $3\times 3$ identity matrix.

##  III Topological magnons

Figure 3: Top panel: magnon bands for (a) $\psi=1.57\pi$, (b) $\psi=1.726\pi$
and (c) $\psi=2\pi$ all with $D=0.7$, which are denoted in Fig. 2(a) as red
points. For the point (b) sits at the phase boundary, the bands are gapless
and there are Dirac points at the $K$ points. Note that the top band in (b)
has small dispersion. Bottom panel: the corresponding magnon bands in a strip
geometry. There are in-gap edge modes in (d) and (f), and nondispersive edge
modes connecting the Dirac points in (e).

We characterize the topological property of the model by the Chern numberHasan
and Kane (2010); Qi and Zhang (2011); Xiao  _et al._ (2010) of the $n$th
magnon band

$\displaystyle C_{n}=\frac{1}{2\pi}\int_{\rm BZ}{\rm d}k_{x}^{\prime}{\rm
d}k_{y}^{\prime}B_{k_{x}^{\prime}k_{y}^{\prime}}^{n},$ (2)

with the Berry curvature of the $n$th band

$\displaystyle B_{k_{x}^{\prime}k_{y}^{\prime}}^{n}\\!\\!={\rm
i}\\!\\!\sum_{n^{\prime}\neq n}\\!\\!\frac{\langle\phi_{n}|\frac{\partial
h(\bm{k}^{\prime})}{\partial
k_{x}^{\prime}}|\phi_{n^{\prime}}\rangle\langle\phi_{n^{\prime}}|\frac{\partial
h(\bm{k}^{\prime})}{\partial
k_{y}^{\prime}}|\phi_{n}\rangle\\!\\!-\\!\\!(k_{x}^{\prime}\leftrightarrow
k_{y}^{\prime})}{(E_{n}-E_{n^{\prime}})^{2}},$ (3)

where $E_{n}$ and $\phi_{n}$ are the eigenvalue and eigenvector of the $n$th
band respectively.

The Heisenberg kagome magnet with only DMI or Kitaev interaction both host
topological magnon excitationsChisnell  _et al._ (2015); Li (2022). The
similar effects on magnon topology of the DMI and Kitaev interaction stem from
their common origin, the spin-orbit coupling. However, due to their different
formalisms of spin operators, there should also be competing effects. Here we
find that the magnon bands of our model carry nontrivial Chern numbers in the
full region of the phase diagram with the tuning of $D$ and $\psi$, as shown
in Fig. 2(a). We characterize the topological phase by $(C_{1},C_{2},C_{3})$,
and the phase diagram is filled by two phases $(1,0,-1)$ and $(-1,0,1)$. Note
that we only show the parameter range $D\in[0.02,1]$ for the phase diagram,
because there are complicated phases and phase transitions in the range
$D\in[0,0.02]$ and the magnitude of the related thermal conductivity is nearly
zero, and these results are deviated from the theme of this paper.

One may note that, the two separated regions of the phase $(-1,0,1)$ have
nearly the same shape, and they may can be related to each other by
$\psi+\pi$. However, it is not exactly the case. The mapping
$\psi\rightarrow\psi+\pi$ will take $(J,K)$ to $(-J,-K)$ and the Zeeman
coupling to the field will only give a energy shift of $h$ for the whole
spectrum, then if we have the eigenvalue of the $n$th band
$E_{n,\psi,D}(\bm{k}^{\prime})=\widetilde{E}_{n}(\bm{k}^{\prime})+h$, we must
also have
$E_{n,\psi+\pi,-D}(\bm{k}^{\prime})=-\widetilde{E}_{n}(\bm{k}^{\prime})+h$.
Hence we have the relation
$(E_{n,\psi,D}(\bm{k}^{\prime})+E_{n,\psi+\pi,-D}(\bm{k}^{\prime}))/2=h$,
which means $E_{n,\psi,D}(\bm{k}^{\prime})$ and
$E_{n,\psi+\pi,-D}(\bm{k}^{\prime})$ are mirror reflections of each other
about the energy axis $E=h$. Thus the closing and re-opening of the band gap,
which denotes a topological phase transition, is simultaneous for the bands
$E_{n,\psi,D}(\bm{k}^{\prime})$ and $E_{n,\psi+\pi,-D}(\bm{k}^{\prime})$.
Therefore, if we have
$E_{n,\psi+\pi,-D}(\bm{k}^{\prime})=E_{n,\psi+\pi,D}(\bm{k}^{\prime})$, the
two separated regions of the phase $(-1,0,1)$ will be related to each other by
$\psi+\pi$ and have exactly the same shape. However, we have checked that the
sign change of $D$ will affect the band shape slightly and break the
simultaneity of the related phase transitions. Consequently, we only have
similar but not the same shape of the two regions of the phase $(-1,0,1)$. Now
the same Chern numbers of the two regions can also be understood easily. The
mirror reflection about a constant energy axis does not change the Chern
numbers of the bands, however, it will change the order of the bands. If the
bands $E_{n,\psi,D}(\bm{k}^{\prime})$ carry Chern numbers $(-1,0,1)$, then the
mirror-reflected bands $E_{n,\psi+\pi,-D}(\bm{k}^{\prime})$ will carry Chern
numbers $(1,0,-1)$. However, the sign change of $D$ will change the sign of
the Chern number of every band. Thus the bands
$E_{n,\psi+\pi,D}(\bm{k}^{\prime})$ will carry Chern numbers $(-1,0,1)$, which
are the same with that of the bands $E_{n,\psi,D}(\bm{k}^{\prime})$.

To gain a deeper insight of the topological property, we calculate the Berry
curvatures of the magnon bands for $\psi=1.57\pi$ and $\psi=2\pi$ both with
$D=0.7$, which are in the phases $(-1,0,1)$ and $(1,0,-1)$ respectively, as
shown in Fig. 2(b). As for $\psi=1.57\pi$ with $D=0.7$, the Berry curvature of
the first band is negative all over the Brillouin zone, and the integral of it
gives Chern number $-1$. While for the third band, its Berry curvature is
always positive and gives Chern number $1$. For $\psi=2\pi$ with $D=0.7$, the
situation is reversed. The sign of the Berry curvature of the second band is
momentum dependent for both cases, and they are not shown here.

Figure 4: (a) The transverse thermal conductivity as function of temperature
for $\psi=1.57\pi$, $\psi=1.726\pi$ and $\psi=2\pi$ all with $D=0.7$, which
are denoted in Fig. 2(a) as red points. The thermal conductivity are always
negative and positive for the phases $(-1,0,1)$ and $(1,0,-1)$ respectively,
and it is always zero for the phase boundary. (b) The transverse thermal
conductivity as function of $\psi$ with $D=0.7$ and $T=10$. There are sign
changes at the phase boundaries. Here we set $\hbar=k_{B}=1$.

In Fig. 3 we show the bulk magnon bands and their corresponding band
structures in a strip geometry for $\psi=1.57\pi$, $\psi=1.726\pi$ and
$\psi=2\pi$ all with $D=0.7$, which are denoted in Fig. 2(a) as red points.
For the point sits at the phase boundary, the band structure is nearly the
same with that of the Heisenberg kagome ferromagnet without DMI, the bands are
gapless and there are Dirac points at the $K$ pointsChisnell  _et al._ (2015),
except that here the top band has small dispersion. Due to the bulk-edge
correspondenceHatsugai (1993), there are in-gap edge modes in the gapped band
structures and nondispersive edge modes connecting the Dirac points.

##  IV Transverse thermal Hall conductivity with sign change

The momentum independent sign structure of the Berry curvatures will lead to
the sign change of the magnon thermal Hall conductivity, as we will show in
this section. The transverse thermal Hall conductivity of magnon can be
calculated asMatsumoto and Murakami (2011a, b); Matsumoto  _et al._ (2014)

$\displaystyle\kappa_{x^{\prime}y^{\prime}}=\frac{k_{B}^{2}T}{(2\pi)^{2}\hbar}\sum_{n}\int_{\rm
BZ}c_{2}(\rho_{n})B_{k_{x}^{\prime}k_{y}^{\prime}}^{n}{\rm
d}k_{x}^{\prime}{\rm d}k_{y}^{\prime},$ (4)

with the sum running over the three bands and the integral is over the first
Brillouin zone. $\rho_{n}=1/(\exp(E_{n}/k_{B}T)-1)$ is the Bose distribution
with $E_{n}$ as the $n$th eigenvalue. $c_{2}$ is given by

$\displaystyle
c_{2}(\rho_{n})=(1+\rho_{n})(\rm{ln}\frac{1+\rho_{n}}{\rho_{n}})^{2}-(\rm{ln}\rho_{n})^{2}-2Li_{2}(-\rho_{n}),$
(5)

where $\rm{Li}_{2}(x)$ is the polylogarithm function of order $2$.

As shown in Fig. 4(a) and 4(b), the thermal Hall conductivity are always
negative and positive for the phases $(-1,0,1)$ and $(1,0,-1)$ respectively,
with the sign changes happen at the phase boundaries, whose conductivity are
always zero. The sign structure of the thermal conductivity stems from the
sign structure of the Berry curvatures of the magnon bands. As for the phase
$(-1,0,1)$, the negative and positive Berry curvatures of the first and third
bands contribute to the Chern numbers -1 and 1 respectively. However, in the
integral of the thermal conductivity, there is a weight $c_{2}(\rho_{n})\geq
0$ which decreases monotonically with energyCao  _et al._ (2015), thus the
negative contribution from the first band is greater than the positive
contribution from the third band, which leads to the negative conductivity.
For the phase $(1,0,-1)$, the situation is reversed and thus its conductivity
is positive. We have checked that, the contributions from the second band are
nearly zero for both phases. From the physical angle, the lowest magnon band
will be thermally occupied maximally at low temperatures, and then dominate
the thermal conductivity.

We note that, the phase $(-1,0,1)$ with negative thermal conductivity resides
around $\psi=0.6\pi$ and $\psi=1.6\pi$ in the phase diagram, where the Kitaev
interaction is relatively large with $\left|K\right|\approx 0.95$. Moreover,
the width of the phase increases with the strength of DMI. These phenomena
will be helpful to distinguish the relative strength of the isotropic
Heisenberg interaction, the Kitaev interaction and the DMI.

##  V Summary

We have studied the topological magnon excitations and related thermal Hall
conductivity in a kagome magnet exhibiting Heisenberg, Kitaev and DM
interactions exposed to a magnetic field. We consider a strong enough filed to
bring the system into fully polarized paramagnetic phase. We find that the
magnon bands carry nontrivial Chern numbers in the full region of the phase
diagram. Furthermore, there are phase transitions relate two topological
phases with opposite Chern numbers, which lead to the sign changes of the
thermal Hall conductivity. In the phase with negative thermal conductivity,
the Kitaev interaction is relatively large and the width of the phase
increases with the strength of DMI. Since the effects of magnon-magnon
interactions are suppressed by the exchange scale divided by the applied field
strength, the paramagnetic phase here is suitable for the exploration of
related physics. Thus we believe that the study here will contribute to the
understanding of related compounds.

## Acknowledgements

We thank Changle Liu for helpful discussions. This work was supported by the
National Natural Science Foundation of China (Grant NO. 12104407) and the
Natural Science Foundation of Zhejiang Province (Grant NO. LQ20A040004).

## Appendix: magnon Hamiltonian matrix

We denote the directions of spins ${\bf S}_{i=1,2,3}$ by their polar angles
$\theta_{1,2,3}$ and azimuthal angles $\phi_{1,2,3}$ in the global frame. The
HP transformation for a spin in its local frame reads

$\displaystyle S_{x}^{0}=$ $\displaystyle\frac{\sqrt{2S}}{2}(a+a^{\dagger}),$
(6) $\displaystyle S_{y}^{0}=$ $\displaystyle\frac{\sqrt{2S}}{2{\rm
i}}(a-a^{\dagger}),$ (7) $\displaystyle S_{z}^{0}=$ $\displaystyle
S-a^{\dagger}a,$ (8)

where $a^{\dagger}$ and $a$ are the magnon creation and annihilation operators
respectively, which obey the boson commutation rules. Then by multiplying a
rotation matrix we get the HP transformation for spin ${\bf S}_{i}$ in the
global frame

$\displaystyle\left(\begin{array}[]{c}S_{ix}\\\\[8.0pt] S_{iy}\\\\[8.0pt]
S_{iz}\end{array}\right)=\left(\begin{array}[]{ccc}\cos\theta_{i}\cos\phi_{i}&-\sin\phi_{i}&\sin\theta_{i}\cos\phi_{i}\\\\[8.0pt]
\cos\theta_{i}\sin\phi_{i}&\cos\phi_{i}&\sin\theta_{i}\sin\phi_{i}\\\\[8.0pt]
-\sin\theta_{i}&0&\cos\theta_{i}\end{array}\right)\left(\begin{array}[]{c}S_{ix}^{0}\\\\[8.0pt]
S_{iy}^{0}\\\\[8.0pt] S_{iz}^{0}\end{array}\right).$ (18)

Substituting $S_{ix,y,z}$ into the Hamiltonian (1) and then do the Fourier
transformation, we get the quadratic Hamiltonian in momentum space

$\displaystyle
H=\frac{1}{2}\sum_{\bm{k}^{\prime}}\Psi_{\bm{k}^{\prime}}^{\dagger}h(\bm{k}^{\prime})\Psi_{\bm{k}^{\prime}},$
(19)

where
$\Psi_{\bm{k}^{\prime}}^{\dagger}=(a_{1\bm{k}^{\prime}}^{\dagger},a_{2\bm{k}^{\prime}}^{\dagger},a_{3\bm{k}^{\prime}}^{\dagger},a_{1-\bm{k}^{\prime}},a_{2-\bm{k}^{\prime}},a_{3-\bm{k}^{\prime}})$.
The magnon Hamiltonian matrix is

$\displaystyle
h(\bm{k}^{\prime})=\left(\begin{array}[]{cc}A_{\bm{k}^{\prime}}&B_{\bm{k}^{\prime}}^{\dagger}\\\\[4.0pt]
B_{\bm{k}^{\prime}}&A_{\bm{k}^{\prime}}^{*}\end{array}\right)S,$ (22)

with $A_{\bm{k}^{\prime}}$ and $B_{\bm{k}^{\prime}}$ $3\times 3$ matrices.
Their elements are as follows

$\displaystyle A_{11}=$ $\displaystyle-2J(c_{z}+e_{z})-2K(c_{8}+e_{7})$
$\displaystyle-2\frac{D}{\sqrt{3}}(c_{D}+e_{D})+\frac{h}{\sqrt{3}}h_{1},$ (23)
$\displaystyle A_{22}=$ $\displaystyle-2J(c_{z}+d_{z})-2K(c_{8}+d_{6})$
$\displaystyle-2\frac{D}{\sqrt{3}}(c_{D}+d_{D})+\frac{h}{\sqrt{3}}h_{2},$ (24)
$\displaystyle A_{33}=$ $\displaystyle-2J(d_{z}+e_{z})-2K(d_{6}+e_{7})$
$\displaystyle-2\frac{D}{\sqrt{3}}(d_{D}+e_{D})+\frac{h}{\sqrt{3}}h_{3},$ (25)

$\displaystyle A_{12}=$ $\displaystyle[J(c_{x}+c_{y}-{\rm i}c_{xy}+{\rm
i}c_{yx})+Kc_{3}+\frac{D}{\sqrt{3}}c_{B}^{*}]\cos\bm{k}^{\prime}\cdot\bm{\delta}_{1},$
(26) $\displaystyle A_{21}=$ $\displaystyle A_{12}^{*},$ (27) $\displaystyle
A_{13}=$ $\displaystyle[J(e_{x}+e_{y}+{\rm i}e_{xy}-{\rm i}e_{yx})$
$\displaystyle+K(e_{2}+e_{5}+{\rm i}e_{10}-{\rm
i}e_{12})+\frac{D}{\sqrt{3}}e_{B}]\cos\bm{k}^{\prime}\cdot\bm{\delta}_{3},$
(28) $\displaystyle A_{31}=$ $\displaystyle A_{13}^{*},$ (29) $\displaystyle
A_{23}=$ $\displaystyle[J(d_{x}+d_{y}-{\rm i}d_{xy}+{\rm i}d_{yx})$
$\displaystyle+K(d_{1}+d_{4}-{\rm i}d_{9}+{\rm
i}d_{11})+\frac{D}{\sqrt{3}}d_{B}^{*}]\cos\bm{k}^{\prime}\cdot\bm{\delta}_{2},$
(30) $\displaystyle A_{32}=$ $\displaystyle A_{23}^{*},$ (31)

$\displaystyle B_{12}=$ $\displaystyle B_{21}=[J(c_{x}-c_{y}-{\rm
i}c_{xy}-{\rm i}c_{yx})$
$\displaystyle+Kc_{3}+\frac{D}{\sqrt{3}}c_{A}]\cos\bm{k}^{\prime}\cdot\bm{\delta}_{1},$
(32) $\displaystyle B_{13}=$ $\displaystyle B_{31}=[J(e_{x}-e_{y}-{\rm
i}e_{xy}-{\rm i}e_{yx})$ $\displaystyle+K(e_{2}-e_{5}-{\rm i}e_{10}-{\rm
i}e_{12})+\frac{D}{\sqrt{3}}e_{A}]\cos\bm{k}^{\prime}\cdot\bm{\delta}_{3},$
(33) $\displaystyle B_{23}=$ $\displaystyle B_{32}=[J(d_{x}-d_{y}-{\rm
i}d_{xy}-{\rm i}d_{yx})$ $\displaystyle+K(d_{1}-d_{4}-{\rm i}d_{9}-{\rm
i}d_{11})+\frac{D}{\sqrt{3}}d_{A}]\cos\bm{k}^{\prime}\cdot\bm{\delta}_{2},$
(34) $\displaystyle B_{11}=$ $\displaystyle B_{22}=B_{33}=0,$ (35)

where

$\displaystyle c_{1}=$
$\displaystyle\cos\theta_{1}\cos\phi_{1}\cos\theta_{2}\cos\phi_{2},$ (36)
$\displaystyle c_{2}=$
$\displaystyle\cos\theta_{1}\sin\phi_{1}\cos\theta_{2}\sin\phi_{2},$ (37)
$\displaystyle c_{3}=$ $\displaystyle\sin\theta_{1}\sin\theta_{2},$ (38)
$\displaystyle c_{4}=$ $\displaystyle\sin\phi_{1}\sin\phi_{2},$ (39)
$\displaystyle c_{5}=$ $\displaystyle\cos\phi_{1}\cos\phi_{2},$ (40)
$\displaystyle c_{6}=$
$\displaystyle\sin\theta_{1}\cos\phi_{1}\sin\theta_{2}\cos\phi_{2},$ (41)
$\displaystyle c_{7}=$
$\displaystyle\sin\theta_{1}\sin\phi_{1}\sin\theta_{2}\sin\phi_{2},$ (42)
$\displaystyle c_{8}=$ $\displaystyle\cos\theta_{1}\cos\theta_{2},$ (43)
$\displaystyle c_{9}=$ $\displaystyle-\cos\theta_{1}\cos\phi_{1}\sin\phi_{2},$
(44) $\displaystyle c_{10}=$
$\displaystyle\cos\theta_{1}\sin\phi_{1}\cos\phi_{2},$ (45) $\displaystyle
c_{11}=$ $\displaystyle-\sin\phi_{1}\cos\theta_{2}\cos\phi_{2},$ (46)
$\displaystyle c_{12}=$ $\displaystyle\cos\phi_{1}\cos\theta_{2}\sin\phi_{2},$
(47)

$\displaystyle c_{x}=$ $\displaystyle c_{1}+c_{2}+c_{3},$ (48) $\displaystyle
c_{y}=$ $\displaystyle c_{4}+c_{5},$ (49) $\displaystyle c_{z}=$
$\displaystyle c_{6}+c_{7}+c_{8},$ (50) $\displaystyle c_{xy}=$ $\displaystyle
c_{9}+c_{10},$ (51) $\displaystyle c_{yx}=$ $\displaystyle c_{11}+c_{12},$
(52)

$\displaystyle c_{13}=$
$\displaystyle\sin\theta_{1}\cos\theta_{2}\sin\phi_{2},$ (53) $\displaystyle
c_{14}=$ $\displaystyle-\cos\theta_{1}\sin\phi_{1}\sin\theta_{2},$ (54)
$\displaystyle c_{15}=$
$\displaystyle\sin\theta_{1}\sin\phi_{1}\cos\theta_{2},$ (55) $\displaystyle
c_{16}=$ $\displaystyle-\cos\theta_{1}\sin\theta_{2}\sin\phi_{2},$ (56)
$\displaystyle c_{17}=$ $\displaystyle-\cos\phi_{1}\sin\theta_{2},$ (57)
$\displaystyle c_{18}=$ $\displaystyle\sin\theta_{1}\cos\phi_{2},$ (58)
$\displaystyle c_{19}=$
$\displaystyle\cos\theta_{1}\cos\phi_{1}\sin\theta_{2},$ (59) $\displaystyle
c_{20}=$ $\displaystyle-\sin\theta_{1}\cos\theta_{2}\cos\phi_{2},$ (60)
$\displaystyle c_{21}=$
$\displaystyle\cos\theta_{1}\sin\theta_{2}\cos\phi_{2},$ (61) $\displaystyle
c_{22}=$ $\displaystyle-\sin\theta_{1}\cos\phi_{1}\cos\theta_{2},$ (62)
$\displaystyle c_{23}=$ $\displaystyle\sin\theta_{1}\sin\phi_{2},$ (63)

$\displaystyle c_{24}=$ $\displaystyle-\sin\phi_{1}\sin\theta_{2},$ (64)
$\displaystyle c_{25}=$
$\displaystyle\cos\theta_{1}\cos\phi_{1}\cos\theta_{2}\sin\phi_{2},$ (65)
$\displaystyle c_{26}=$
$\displaystyle-\cos\theta_{1}\sin\phi_{1}\cos\theta_{2}\cos\phi_{2},$ (66)
$\displaystyle c_{27}=$ $\displaystyle\cos\phi_{1}\sin\phi_{2},$ (67)
$\displaystyle c_{28}=$ $\displaystyle-\sin\phi_{1}\cos\phi_{2},$ (68)
$\displaystyle c_{29}=$
$\displaystyle\sin\theta_{1}\cos\phi_{1}\sin\theta_{2}\sin\phi_{2},$ (69)
$\displaystyle c_{30}=$
$\displaystyle-\sin\theta_{1}\sin\phi_{1}\sin\theta_{2}\cos\phi_{2},$ (70)
$\displaystyle c_{31}=$ $\displaystyle\cos\theta_{1}\cos\phi_{1}\cos\phi_{2},$
(71) $\displaystyle c_{32}=$
$\displaystyle\cos\theta_{1}\sin\phi_{1}\sin\phi_{2},$ (72) $\displaystyle
c_{33}=$ $\displaystyle-\sin\phi_{1}\cos\theta_{2}\sin\phi_{2},$ (73)
$\displaystyle c_{34}=$
$\displaystyle-\cos\phi_{1}\cos\theta_{2}\cos\phi_{2},$ (74)

$\displaystyle c_{xxx}=$ $\displaystyle c_{13}+c_{14},$ (75) $\displaystyle
c_{xzz}=$ $\displaystyle c_{15}+c_{16},$ (76) $\displaystyle c_{yxx}=$
$\displaystyle c_{19}+c_{20},$ (77) $\displaystyle c_{yzz}=$ $\displaystyle
c_{21}+c_{22},$ (78) $\displaystyle c_{zxx}=$ $\displaystyle c_{25}+c_{26},$
(79) $\displaystyle c_{zyy}=$ $\displaystyle c_{27}+c_{28},$ (80)
$\displaystyle c_{zzz}=$ $\displaystyle c_{29}+c_{30},$ (81) $\displaystyle
c_{zxy}=$ $\displaystyle c_{31}+c_{32},$ (82) $\displaystyle c_{zyx}=$
$\displaystyle c_{33}+c_{34},$ (83)

$\displaystyle c_{A}=$ $\displaystyle c_{xxx}+c_{yxx}+c_{zxx}-c_{zyy}$
$\displaystyle-{\rm i}(c_{17}+c_{24}+c_{zyx}+c_{18}+c_{23}+c_{zxy}),$ (84)
$\displaystyle c_{B}=$ $\displaystyle c_{xxx}+c_{yxx}+c_{zxx}+c_{zyy}$
$\displaystyle-{\rm i}(c_{17}+c_{24}+c_{zyx}-c_{18}-c_{23}-c_{zxy}),$ (85)
$\displaystyle c_{D}=$ $\displaystyle c_{xzz}+c_{yzz}+c_{zzz},$ (86)

and change the corresponding subscripts in $\theta_{1,2}$ and $\phi_{1,2}$ to
$\theta_{2,3}$ and $\phi_{2,3}$, we get the corresponding expressions for
$d_{i}$ with

$\displaystyle i=$ $\displaystyle 1-34,x,y,z,xy,yx,xxx,xzz,yxx,yzz,$
$\displaystyle zxx,zyy,zzz,zxy,zyx,A,B,D.$ (87)

Similarly, change $\theta_{1,2}$ and $\phi_{1,2}$ to $\theta_{3,1}$ and
$\phi_{3,1}$, we get the corresponding expressions for $e_{i}$. And we have

$\displaystyle
h_{1}=\sin\theta_{1}\cos\phi_{1}+\sin\theta_{1}\sin\phi_{1}+\cos\theta_{1},$
(88) $\displaystyle
h_{2}=\sin\theta_{2}\cos\phi_{2}+\sin\theta_{2}\sin\phi_{2}+\cos\theta_{2},$
(89) $\displaystyle
h_{3}=\sin\theta_{3}\cos\phi_{3}+\sin\theta_{3}\sin\phi_{3}+\cos\theta_{3},$
(90)

Note that for the polarized phase here, we have
$\theta_{1}=\theta_{2}=\theta_{3}=\arctan\sqrt{2}$ and
$\phi_{1}=\phi_{2}=\phi_{3}=\pi/4$. The vectors $\bm{\delta}_{1,2,3}$ are the
NN vectors of the kagome lattice with $\bm{\delta}_{1}=(1/2,0)$,
$\bm{\delta}_{2}=(-1,\sqrt{3})/4$, $\bm{\delta}_{3}=(-1,-\sqrt{3})/4$, which
are defined in the new 2D frame $x^{\prime}y^{\prime}$. Note that to get the
eigenvalues and eigenvectors of bosonic quadratic Hamiltonian, we need to
diagonalize the matrix $I_{-}h(\bm{k}^{\prime})$ instead of
$h(\bm{k}^{\prime})$, where

$\displaystyle I_{-}=\left(\begin{array}[]{cc}I&0\\\\[4.0pt]
0&-I\end{array}\right),$ (93)

with $I$ the $3\times 3$ identity matrix .

## References

  * Rousochatzakis  _et al._ (2015) Ioannis Rousochatzakis, Johannes Reuther, Ronny Thomale, Stephan Rachel, and N. B. Perkins, “Phase diagram and quantum order by disorder in the kitaev ${K}_{1}-{K}_{2}$ honeycomb magnet,” Phys. Rev. X 5, 041035 (2015).
  * Lee  _et al._ (2014) Eric Kin-Ho Lee, Robert Schaffer, Subhro Bhattacharjee, and Yong Baek Kim, “Heisenberg-kitaev model on the hyperhoneycomb lattice,” Phys. Rev. B 89, 045117 (2014).
  * Nasu  _et al._ (2014) Joji Nasu, Masafumi Udagawa, and Yukitoshi Motome, “Vaporization of kitaev spin liquids,” Phys. Rev. Lett. 113, 197205 (2014).
  * Takayama  _et al._ (2015) T. Takayama, A. Kato, R. Dinnebier, J. Nuss, H. Kono, L. S. I. Veiga, G. Fabbris, D. Haskel, and H. Takagi, “Hyperhoneycomb iridate $\beta\text{$-$}{\mathrm{li}}_{2}{\mathrm{iro}}_{3}$ as a platform for kitaev magnetism,” Phys. Rev. Lett. 114, 077202 (2015).
  * Jahromi  _et al._ (2016) Saeed S. Jahromi, Mehdi Kargarian, S. Farhad Masoudi, and Abdollah Langari, “Topological spin liquids in the ruby lattice with anisotropic kitaev interactions,” Phys. Rev. B 94, 125145 (2016).
  * Jahromi  _et al._ (2018) Saeed S. Jahromi, Román Orús, Mehdi Kargarian, and Abdollah Langari, “Infinite projected entangled-pair state algorithm for ruby and triangle-honeycomb lattices,” Phys. Rev. B 97, 115161 (2018).
  * Kishimoto  _et al._ (2018) Masanori Kishimoto, Katsuhiro Morita, Yukihiro Matsubayashi, Shigetoshi Sota, Seiji Yunoki, and Takami Tohyama, “Ground state phase diagram of the kitaev-heisenberg model on a honeycomb-triangular lattice,” Phys. Rev. B 98, 054411 (2018).
  * Kimchi and Vishwanath (2014) Itamar Kimchi and Ashvin Vishwanath, “Kitaev-heisenberg models for iridates on the triangular, hyperkagome, kagome, fcc, and pyrochlore lattices,” Phys. Rev. B 89, 014414 (2014).
  * Dzyaloshinsky (1958) I. Dzyaloshinsky, “A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics,” Journal of Physics and Chemistry of Solids 4, 241–255 (1958).
  * Moriya (1960) Tôru Moriya, “Anisotropic superexchange interaction and weak ferromagnetism,” Phys. Rev. 120, 91–98 (1960).
  * Gao  _et al._ (2019) Yong Hao Gao, Ciarán Hickey, Tao Xiang, Simon Trebst, and Gang Chen, “Thermal hall signatures of non-kitaev spin liquids in honeycomb kitaev materials,” Phys. Rev. Research 1, 013014 (2019).
  * Zhang  _et al._ (2021) Li-Chuan Zhang, Fengfeng Zhu, Dongwook Go, Fabian R. Lux, Flaviano José dos Santos, Samir Lounis, Yixi Su, Stefan Blügel, and Yuriy Mokrousov, “Interplay of dzyaloshinskii-moriya and kitaev interactions for magnonic properties of heisenberg-kitaev honeycomb ferromagnets,” Phys. Rev. B 103, 134414 (2021).
  * Joshi (2018) Darshan G. Joshi, “Topological excitations in the ferromagnetic kitaev-heisenberg model,” Phys. Rev. B 98, 060405 (2018).
  * Chern  _et al._ (2021) Li Ern Chern, Emily Z. Zhang, and Yong Baek Kim, “Sign structure of thermal hall conductivity and topological magnons for in-plane field polarized kitaev magnets,” Phys. Rev. Lett. 126, 147201 (2021).
  * McClarty  _et al._ (2018) P. A. McClarty, X.-Y. Dong, M. Gohlke, J. G. Rau, F. Pollmann, R. Moessner, and K. Penc, “Topological magnons in kitaev magnets at high fields,” Phys. Rev. B 98, 060404 (2018).
  * Li (2022) Kangkang Li, “Thermal hall conductivity with sign change in the heisenberg-kitaev kagome magnet,” arXiv:2209.04133 (2022), 10.48550/arXiv.2209.04133.
  * Katsura  _et al._ (2010) Hosho Katsura, Naoto Nagaosa, and Patrick A. Lee, “Theory of the thermal hall effect in quantum magnets,” Phys. Rev. Lett. 104, 066403 (2010).
  * Onose  _et al._ (2010) Y. Onose, T. Ideue, H. Katsura, Y. Shiomi, N. Nagaosa, and Y. Tokura, “Observation of the magnon hall effect,” Science 329, 297–299 (2010), https://www.science.org/doi/pdf/10.1126/science.1188260 .
  * Hirschberger  _et al._ (2015) Max Hirschberger, Robin Chisnell, Young S. Lee, and N. P. Ong, “Thermal hall effect of spin excitations in a kagome magnet,” Phys. Rev. Lett. 115, 106603 (2015).
  * Chisnell  _et al._ (2015) R. Chisnell, J. S. Helton, D. E. Freedman, D. K. Singh, R. I. Bewley, D. G. Nocera, and Y. S. Lee, “Topological magnon bands in a kagome lattice ferromagnet,” Phys. Rev. Lett. 115, 147201 (2015).
  * Owerre (2016) S A Owerre, “A first theoretical realization of honeycomb topological magnon insulator,” Journal of Physics: Condensed Matter 28, 386001 (2016).
  * Zhang  _et al._ (2013) Lifa Zhang, Jie Ren, Jian-Sheng Wang, and Baowen Li, “Topological magnon insulator in insulating ferromagnet,” Phys. Rev. B 87, 144101 (2013).
  * Li  _et al._ (2017) Kangkang Li, Chenyuan Li, Jiangping Hu, Yuan Li, and Chen Fang, “Dirac and nodal line magnons in three-dimensional antiferromagnets,” Phys. Rev. Lett. 119, 247202 (2017).
  * Li and Hu (2017) Kang-Kang Li and Jiang-Ping Hu, “Weyl and nodal ring magnons in three-dimensional honeycomb lattices,” Chinese Physics Letters 34, 077501 (2017).
  * Cao  _et al._ (2015) Xiaodong Cao, Kai Chen, and Dahai He, “Magnon hall effect on the lieb lattice,” Journal of Physics: Condensed Matter 27, 166003 (2015).
  * Morita  _et al._ (2018) Katsuhiro Morita, Masanori Kishimoto, and Takami Tohyama, “Ground-state phase diagram of the kitaev-heisenberg model on a kagome lattice,” Phys. Rev. B 98, 134437 (2018).
  * Yang  _et al._ (2020) Yang Yang, Natalia B. Perkins, Fulya Koç, Chi-Huei Lin, and Ioannis Rousochatzakis, “Quantum-classical crossover in the spin-$\frac{1}{2}$ heisenberg-kitaev kagome magnet,” Phys. Rev. Research 2, 033217 (2020).
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff, “Field dependence of the intrinsic domain magnetization of a ferromagnet,” Phys. Rev. 58, 1098–1113 (1940).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, “Colloquium: Topological insulators,” Rev. Mod. Phys. 82, 3045–3067 (2010).
  * Qi and Zhang (2011) Xiao-Liang Qi and Shou-Cheng Zhang, “Topological insulators and superconductors,” Rev. Mod. Phys. 83, 1057–1110 (2011).
  * Xiao  _et al._ (2010) Di Xiao, Ming-Che Chang, and Qian Niu, “Berry phase effects on electronic properties,” Rev. Mod. Phys. 82, 1959–2007 (2010).
  * Hatsugai (1993) Yasuhiro Hatsugai, “Chern number and edge states in the integer quantum hall effect,” Phys. Rev. Lett. 71, 3697–3700 (1993).
  * Matsumoto and Murakami (2011a) Ryo Matsumoto and Shuichi Murakami, “Theoretical prediction of a rotating magnon wave packet in ferromagnets,” Phys. Rev. Lett. 106, 197202 (2011a).
  * Matsumoto and Murakami (2011b) Ryo Matsumoto and Shuichi Murakami, “Rotational motion of magnons and the thermal hall effect,” Phys. Rev. B 84, 184406 (2011b).
  * Matsumoto  _et al._ (2014) Ryo Matsumoto, Ryuichi Shindou, and Shuichi Murakami, “Thermal hall effect of magnons in magnets with dipolar interaction,” Phys. Rev. B 89, 054420 (2014).


# Autoregressive Denoising Diffusion Models for Multivariate Probabilistic
Time Series Forecasting

Kashif Rasul  Calvin Seward  Ingmar Schuster  Roland Vollgraf

###### Abstract

In this work, we propose TimeGrad, an autoregressive model for multivariate
probabilistic time series forecasting which samples from the data distribution
at each time step by estimating its gradient. To this end, we use diffusion
probabilistic models, a class of latent variable models closely connected to
score matching and energy-based methods. Our model learns gradients by
optimizing a variational bound on the data likelihood and at inference time
converts white noise into a sample of the distribution of interest through a
Markov chain using Langevin sampling. We demonstrate experimentally that the
proposed autoregressive denoising diffusion model is the new state-of-the-art
multivariate probabilistic forecasting method on real-world data sets with
thousands of correlated dimensions. We hope that this method is a useful tool
for practitioners and lays the foundation for future research in this area.

Time Series and Sequences, Generative Models

##  1 Introduction

Classical time series forecasting methods such as those in (Hyndman &
Athanasopoulos, 2018) typically provide univariate point forecasts, require
hand-tuned features to model seasonality, and are trained individually on each
time series. Deep learning based time series models (Benidis et al., 2020) are
popular alternatives due to their end-to-end training of a global model, ease
of incorporating exogenous covariates, and automatic feature extraction
abilities. The task of modeling uncertainties is of vital importance for
downstream problems that use these forecasts for (business) decision making.
More often the individual time series for a problem data set are statistically
dependent on each other. Ideally, deep learning models need to incorporate
this inductive bias in the form of multivariate (Tsay, 2014) probabilistic
methods to provide accurate forecasts.

To model the full predictive distribution, methods typically resort to
tractable distribution classes or some type of low-rank approximations,
regardless of the true data distribution. To model the distribution in a
general fashion, one needs probabilistic methods with tractable likelihoods.
Till now several deep learning methods have been proposed for this purpose
such as autoregressive (van den Oord et al., 2016c) or generative ones based
on normalizing flows (Papamakarios et al., 2019) which can learn flexible
models of high dimensional multivariate time series. Even if the full
likelihood is not be tractable, one can often optimize a tractable lower bound
to the likelihood. But still, these methods require a certain structure in the
functional approximators, for example on the determinant of the Jacobian (Dinh
et al., 2017) for normalizing flows. _Energy-based models_ (EBM) (Hinton,
2002; LeCun et al., 2006) on the other hand have a much less restrictive
functional form. They approximate the unnormalized log-probability so that
density estimation reduces to a non-linear regression problem. EBMs have been
shown to perform well in learning high dimensional data distributions at the
cost of being difficult to train (Song & Kingma, 2021).

In this work, we propose autoregressive EBMs to solve the multivariate
probabilistic time series forecasting problem via a model we call TimeGrad and
show that not only are we able to train such a model with all the inductive
biases of probabilistic time series forecasting, but this model performs
exceptionally well when compared to other modern methods. This autoregressive-
EBM combination retains the power of autoregressive models, such as good
performance in extrapolation into the future, with the flexibility of EBMs as
a general purpose high-dimensional distribution model, while remaining
computationally tractable.

The paper is organized as follows. In Section 2 we first set up the notation
and detail the EBM of (Ho et al., 2020) which forms the basis of our per time-
step distribution model. Section 3 introduces the multivariate probabilistic
time series problem and we detail the TimeGrad model. The experiments with
extensive results are detailed in Section 4. We cover related work in Section
5 and conclude with some discussion in Section 6.

##  2 Diffusion Probabilistic Model

Let $\mathbf{x}^{0}\sim q_{\mathcal{X}}(\mathbf{x}^{0})$ denote the
multivariate training vector from some input space
$\mathcal{X}=\mathbb{R}^{D}$ and let $p_{\theta}(\mathbf{x}^{0})$ denote the
probability density function (PDF) which aims to approximate
$q_{\mathcal{X}}(\mathbf{x}^{0})$ and allows for easy sampling. Diffusion
models (Sohl-Dickstein et al., 2015) are latent variable models of the form
$p_{\theta}(\mathbf{x}^{0}):=\int
p_{\theta}(\mathbf{x}^{0:N})\,\mathrm{d}\mathbf{x}^{1:N}$, where
$\mathbf{x}^{1},\ldots,\mathbf{x}^{N}$ are latents of dimension
$\mathbb{R}^{D}$. Unlike in variational autoencoders (Kingma & Welling, 2019)
the approximate posterior $q(\mathbf{x}^{1:N}|\mathbf{x}^{0})$,

$q(\mathbf{x}^{1:N}|\mathbf{x}^{0})=\Pi_{n=1}^{N}q(\mathbf{x}^{n}|\mathbf{x}^{n-1})$

is not trainable but fixed to a Markov chain (called the _forward_ process)
that gradually adds Gaussian noise to the signal:

$q(\mathbf{x}^{n}|\mathbf{x}^{n-1}):=\mathcal{N}(\mathbf{x}^{n};\sqrt{1-\beta_{n}}\mathbf{x}^{n-1},\beta_{n}\mathbf{I}).$

The forward process uses an increasing variance schedule
$\beta_{1},\ldots,\beta_{N}$ with $\beta_{n}\in(0,1)$. The joint distribution
$p_{\theta}(\mathbf{x}^{0:N})$ is called the _reverse_ process, and is defined
as a Markov chain with learned Gaussian transitions starting with
$p(\mathbf{x}^{N})=\mathcal{N}(\mathbf{x}^{N};\mathbf{0},\mathbf{I})$, where
each subsequent transition of

$p_{\theta}(\mathbf{x}^{0:N}):=p(\mathbf{x}^{N})\Pi_{n=N}^{1}p_{\theta}(\mathbf{x}^{n-1}|\mathbf{x}^{n})$

is given by a parametrization of our choosing denoted by

$p_{\theta}(\mathbf{x}^{n-1}|\mathbf{x}^{n}):=\mathcal{N}(\mathbf{x}^{n-1};\mu_{\theta}(\mathbf{x}^{n},n),\Sigma_{\theta}(\mathbf{x}^{n},n)\mathbf{I}),$
(1)

with shared parameters $\theta$. Both
$\mu_{\theta}:\mathbb{R}^{D}\times\mathbb{N}\to\mathbb{R}^{D}$ and
$\Sigma_{\theta}:\mathbb{R}^{D}\times\mathbb{N}\to\mathbb{R}^{+}$ take two
inputs, namely the variable $\mathbf{x}^{n}\in\mathbb{R}^{D}$ as well as the
noise index $n\in\mathbb{N}$. The goal of
$p_{\theta}(\mathbf{x}^{n-1}|\mathbf{x}^{n})$ is to eliminate the Gaussian
noise added in the diffusion process. The parameters $\theta$ are learned to
fit the data distribution $q_{\mathcal{X}}(\mathbf{x}^{0})$ by minimizing the
negative log-likelihood via a variational bound using Jensen’s inequality:

$\begin{split}\min_{\theta}\mathbb{E}_{q(\mathbf{x}^{0})}[-\log
p_{\theta}(\mathbf{x}^{0})]\leq\\\
\min_{\theta}\mathbb{E}_{q(\mathbf{x}^{0:N})}[-\log
p_{\theta}(\mathbf{x}^{0:N})+\log
q(\mathbf{x}^{1:N}|\mathbf{x}^{0})].\end{split}$

This upper bound can be shown to be equal to

$\min_{\theta}\mathbb{E}_{q(\mathbf{x}^{0:N})}\left[-\log
p(\mathbf{x}^{N})-\sum_{n=1}^{N}\log\frac{p_{\theta}(\mathbf{x}^{n-1}|\mathbf{x}^{n})}{q(\mathbf{x}^{n}|\mathbf{x}^{n-1})}\right].$
(2)

As shown by (Ho et al., 2020), a property of the forward process is that it
admits sampling $\mathbf{x}^{n}$ at any arbitrary noise level $n$ in closed
form, since if $\alpha_{n}:=1-\beta_{n}$ and
$\bar{\alpha}_{n}:=\Pi_{i=1}^{n}\alpha_{i}$ its cumulative product, we have:

$q(\mathbf{x}^{n}|\mathbf{x}^{0})=\mathcal{N}(\mathbf{x}^{n};\sqrt{\bar{\alpha}_{n}}\mathbf{x}^{0},(1-\bar{\alpha}_{n})\mathbf{I}).$
(3)

By using the fact that these processes are Markov chains, the objective in (2)
can be written as the KL-divergence between Gaussian distributions:

$-\log
p_{\theta}(\mathbf{x}^{0}|\mathbf{x}^{1})+D_{\mathrm{KL}}(q(\mathbf{x}^{N}|\mathbf{x}^{0})||p(\mathbf{x}^{N}))\\\
+\sum_{n=2}^{N}D_{\mathrm{KL}}(q(\mathbf{x}^{n-1}|\mathbf{x}^{n},\mathbf{x}^{0})||p_{\theta}(\mathbf{x}^{n-1}|\mathbf{x}^{n})),$
(4)

and (Ho et al., 2020) shows that by the property (3) the forward process
posterior in these KL divergences when conditioned on $\mathbf{x}^{0}$, i.e.
$q(\mathbf{x}^{n-1}|\mathbf{x}^{n},\mathbf{x}^{0})$ are tractable given by

$q(\mathbf{x}^{n-1}|\mathbf{x}^{n},\mathbf{x}^{0})={\cal{N}}(\mathbf{x}^{n-1};\tilde{\mu}_{n}(\mathbf{x}^{n},\mathbf{x}^{0}),\tilde{\beta}_{n}\mathbf{I}),$

where

$\tilde{\mu}_{n}(\mathbf{x}^{n},\mathbf{x}^{0}):=\frac{\sqrt{\bar{\alpha}_{n-1}}\beta_{n}}{1-\bar{\alpha}_{n}}\mathbf{x}^{0}+\frac{\sqrt{\alpha_{n}}(1-\bar{\alpha}_{n-1})}{1-\bar{\alpha}_{n}}\mathbf{x}^{n}$

and

$\tilde{\beta}_{n}:=\frac{1-\bar{\alpha}_{n-1}}{1-\bar{\alpha}_{n}}\beta_{n}.$
(5)

Further, (Ho et al., 2020) shows that the KL-divergence between Gaussians can
be written as:

$D_{\mathrm{KL}}(q(\mathbf{x}^{n-1}|\mathbf{x}^{n},\mathbf{x}^{0})||p_{\theta}(\mathbf{x}^{n-1}|\mathbf{x}^{n}))=\\\
\mathbb{E}_{q}\left[\frac{1}{2\Sigma_{\theta}}\|\tilde{\mu}_{n}(\mathbf{x}^{n},\mathbf{x}^{0})-\mu_{\theta}(\mathbf{x}^{n},n)\|^{2}\right]+C,$
(6)

where $C$ is a constant which does not depend on $\theta$. So instead of a
parametrization (1) of $p_{\theta}$ that predicts $\tilde{\mu}$, one can
instead use the property (3) to write
$\mathbf{x}^{n}(\mathbf{x}^{0},\mathbf{\epsilon})=\sqrt{\bar{\alpha}_{n}}\mathbf{x}^{0}+\sqrt{1-\bar{\alpha}_{n}}\mathbf{\epsilon}$
for $\mathbf{\epsilon}\sim{\cal{N}}(\mathbf{0},\mathbf{I})$ and the formula
for $\tilde{\mu}$ to obtain that $\mu_{\theta}$ must predict
$(\mathbf{x}^{n}-\beta_{n}\mathbf{\epsilon}/\sqrt{1-\bar{\alpha}_{n}})/\sqrt{\alpha_{n}}$,
but since $\mathbf{x}^{n}$ is available to the network, we can choose:

$\mu_{\theta}(\mathbf{x}^{n},n)=\frac{1}{\sqrt{\alpha_{n}}}\left(\mathbf{x}^{n}-\frac{\beta_{n}}{\sqrt{1-\bar{\alpha}_{n}}}\mathbf{\epsilon}_{\theta}(\mathbf{x}^{n},n)\right),$

where $\mathbf{\epsilon}_{\theta}$ is a network which predicts
$\mathbf{\epsilon}\sim{\cal{N}}(\mathbf{0},\mathbf{I})$ from $\mathbf{x}^{n}$,
so that the objective simplifies to:

$\mathbb{E}_{\mathbf{x}^{0},\mathbf{\epsilon}}\left[\frac{\beta_{n}^{2}}{2\Sigma_{\theta}\alpha_{n}(1-\bar{\alpha}_{n})}\|\mathbf{\epsilon}-\mathbf{\epsilon}_{\theta}(\sqrt{\bar{\alpha}_{n}}\mathbf{x}^{0}+\sqrt{1-\bar{\alpha}_{n}}\mathbf{\epsilon},n)\|^{2}\right]$
(7)

resembling the loss in Noise Conditional Score Networks (Song & Ermon, 2019,
2020) using score matching. Once trained, to sample from the reverse process
$\mathbf{x}^{n-1}\sim p_{\theta}(\mathbf{x}^{n-1}|\mathbf{x}^{n})$ (1) we can
compute

$\mathbf{x}^{n-1}=\frac{1}{\sqrt{\alpha_{n}}}\left(\mathbf{x}^{n}-\frac{\beta_{n}}{\sqrt{1-\bar{\alpha}_{n}}}\mathbf{\epsilon}_{\theta}(\mathbf{x}^{n},n)\right)+\sqrt{\Sigma_{\theta}}\mathbf{z}$

where $\mathbf{z}\sim{\cal{N}}(\mathbf{0},\mathbf{I})$ for $n=N,\ldots,2$ and
$\mathbf{z}=\mathbf{0}$ when $n=1$. The full sampling procedure for
$\mathbf{x}^{0}$, starting from white noise sample $\mathbf{x}^{N}$, resembles
Langevin dynamics where we sample from the most noise-perturbed distribution
and reduce the magnitude of the noise scale until we reach the smallest one.

##  3 TimeGrad Method

We denote the entities of a multivariate time series by
$x_{i,t}^{0}\in\mathbb{R}$ for $i\in\\{1,\ldots,D\\}$ where $t$ is the time
index. Thus the multivariate vector at time $t$ is given by
$\mathbf{x}_{t}^{0}\in\mathbb{R}^{D}$. We are tasked with predicting the
multivariate distribution some given prediction time steps into the future and
so in what follows consider time series with $t\in[1,T]$, sampled from the
complete time series history of the training data, where we will split this
contiguous sequence into a context window of size $[1,t_{0})$ and prediction
interval $[t_{0},T]$, reminiscent of seq-to-seq models (Sutskever et al.,
2014) in language modeling.

In the univariate probabilistic DeepAR model (Salinas et al., 2019b), the log-
likelihood of each entity $x^{0}_{i,t}$ at a time step $t\in[t_{0},T]$ is
maximized over an individual time series’ prediction window. This is done with
respect to the parameters of some chosen distributional model via the state of
an RNN derived from its previous time step $x^{0}_{i,t-1}$ and its
corresponding covariates $\mathbf{c}_{i,t-1}$. The emission distribution
model, which is typically Gaussian for real-valued data or negative binomial
for count data, is selected to best match the statistics of the time series
and the network incorporates activation functions that satisfy the constraints
of the distribution’s parameters, e.g. a softplus() for the scale parameter of
the Gaussian.

A straightforward time series model for multivariate real-valued data could
use a factorizing output distribution instead. Shared parameters can then
learn patterns across the individual time series entities through the temporal
component — but the model falls short of capturing dependencies in the
emissions of the model. For this, a full joint distribution at each time step
has to be modeled, for example by using a multivariate Gaussian. However,
modeling the full covariance matrix not only increases the number of
parameters of the neural network by $O(D^{2})$, making learning difficult but
computing the loss is $O(D^{3})$ making it impractical. Furthermore,
statistical dependencies for such distributions would be limited to second-
order effects. Approximating Gaussians with low-rank covariance matrices do
work however and these models are referred to as Vec-LSTM in (Salinas et al.,
2019a).

Instead, in this work we propose TimeGrad which aims to learn a model of the
conditional distribution of the future time steps of a multivariate time
series given its past and covariates as:

$q_{\mathcal{X}}(\mathbf{x}_{t_{0}:T}^{0}|\mathbf{x}_{1:t_{0}-1}^{0},\mathbf{c}_{1:T})=\Pi_{t=t_{0}}^{T}q_{\mathcal{X}}(\mathbf{x}_{t}^{0}|\mathbf{x}_{1:t-1}^{0},\mathbf{c}_{1:T}),$
(8)

were we assume that the covariates are known for all the time points and each
factor is learned via a _conditional_ denoising diffusion model introduced
above. To model the temporal dynamics we employ the autoregressive recurrent
neural network (RNN) architecture from (Graves, 2013; Sutskever et al., 2014)
which utilizes the LSTM (Hochreiter & Schmidhuber, 1997) or GRU (Chung et al.,
2014) to encode the time series sequence up to time point $t$, given the
covariates $\mathbf{c}_{t}$, via the updated hidden state $\mathbf{h}_{t}$:

$\mathbf{h}_{t}=\mathrm{RNN}_{\theta}(\mathtt{concat}(\mathbf{x}_{t}^{0},\mathbf{c}_{t}),\mathbf{h}_{t-1}),$
(9)

where $\mathrm{RNN}_{\theta}$ is a multi-layer LSTM or GRU parameterized by
shared weights $\theta$ and $\mathbf{h}_{0}=\mathbf{0}$. Thus we can
approximate (8) by the model

$\Pi_{t=t_{0}}^{T}p_{\theta}(\mathbf{x}_{t}^{0}|\mathbf{h}_{t-1}),$ (10)

where now $\theta$ comprises the weights of the RNN as well as denoising
diffusion model. This model is autoregressive as it consumes the observations
at the time step $t-1$ as input to learn the distribution of, or sample, the
next time step as shown in Figure 1.

###  3.1 Training

Training is performed by randomly sampling context and adjoining prediction
sized windows from the training time series data and optimizing the parameters
$\theta$ that minimize the negative log-likelihood of the model (10):

$\sum_{t=t_{0}}^{T}-\log p_{\theta}(\mathbf{x}_{t}^{0}|\mathbf{h}_{t-1}),$

starting with the hidden state $\mathbf{h}_{t_{0}-1}$ obtained by running the
RNN on the chosen context window. Via a similar derivation as in the previous
section, we have that the conditional variant of the objective (4) for time
step $t$ and noise index $n$ is given by the following simplification of (7)
(Ho et al., 2020):

$\mathbb{E}_{\mathbf{x}_{t}^{0},\epsilon,n}\left[\|\mathbf{\epsilon}-\mathbf{\epsilon}_{\theta}(\sqrt{\bar{\alpha}_{n}}\mathbf{x}^{0}_{t}+\sqrt{1-\bar{\alpha}_{n}}\mathbf{\epsilon},\mathbf{h}_{t-1},n)\|^{2}\right],$

when we choose the variance in (1) to be $\Sigma_{\theta}=\tilde{\beta}_{n}$
(5), where now the $\epsilon_{\theta}$ network is also _conditioned_ on the
hidden state. Algorithm 1 is the training procedure for each time step in the
prediction window using this objective.

Algorithm 1 Training for each time series step $t\in[t_{0},T]$

Input: data $\mathbf{x}_{t}^{0}\sim q_{\cal{X}}(\mathbf{x}_{t}^{0})$ and state
$\mathbf{h}_{t-1}$

repeat

Initialize $n\sim\mathrm{Uniform}({1,\ldots,N})$ and
$\epsilon\sim{\cal{N}}(\mathbf{0},\mathbf{I})$ Take gradient step on
$\nabla_{\theta}\|\mathbf{\epsilon}-\mathbf{\epsilon}_{\theta}(\sqrt{\bar{\alpha}_{n}}\mathbf{x}^{0}_{t}+\sqrt{1-\bar{\alpha}_{n}}\mathbf{\epsilon},\mathbf{h}_{t-1},n)\|^{2}$

until converged

Figure 1: TimeGrad schematic: an RNN _conditioned_ diffusion probabilistic
model at some time $t-1$ depicting the fixed forward process that adds
Gaussian noise and the learned reverse processes.

###  3.2 Inference

After training, we wish to predict for each time series in our data set some
prediction steps into the future and compare with the corresponding test set
time series. As in training, we run the RNN over the last context sized window
of the training set to obtain the hidden state $\mathbf{h}_{T}$ via (9). Then
we follow the sampling procedure in Algorithm 2 to obtain a sample
$\mathbf{x}_{T+1}^{0}$ of the next time step, which we can pass
autoregressively to the RNN together with the covariates $\mathbf{c}_{T+1}$ to
obtain the next hidden state $\mathbf{h}_{T+1}$ and repeat until the desired
forecast horizon has been reached. This process of sampling trajectories from
the “warm-up” state $\mathbf{h}_{T}$ can be repeated many times (e.g. $S=100$)
to obtain empirical quantiles of the uncertainty of our predictions.

Algorithm 2 Sampling $\mathbf{x}_{t}^{0}$ via annealed Langevin dynamics

Input: noise $\mathbf{x}_{t}^{N}\sim{\cal{N}}(\mathbf{0},\mathbf{I})$ and
state $\mathbf{h}_{t-1}$

for $n=N$ to $1$ do

if $n>1$ then

$\mathbf{z}\sim{\cal{N}}(\mathbf{0},\mathbf{I})$

else

$\mathbf{z}=\mathbf{0}$

end if

$\mathbf{x}_{t}^{n-1}=\frac{1}{\sqrt{\alpha_{n}}}(\mathbf{x}^{n}_{t}-\frac{\beta_{n}}{\sqrt{1-\bar{\alpha}_{n}}}\mathbf{\epsilon}_{\theta}(\mathbf{x}^{n}_{t},\mathbf{h}_{t-1},n))+\sqrt{\Sigma_{\theta}}\mathbf{z}$

end for

Return: $\mathbf{x}_{t}^{0}$

###  3.3 Scaling

In real-world data, the magnitudes of different time series entities can vary
drastically. To normalize scales, we divide each time series entity by their
context window mean (or $1$ if it’s zero) before feeding it into the model. At
inference, the samples are then multiplied by the same mean values to match
the original scale. This rescaling technique simplifies the problem for the
model, which is reflected in significantly improved empirical performance as
shown in (Salinas et al., 2019b). The other method of a short-cut connection
from the input to the output of the function approximator, as done in the
multivariate point forecasting method LSTNet (Lai et al., 2018), is not
applicable here.

###  3.4 Covariates

We employ embeddings for categorical features (Charrington, 2018), that allows
for relationships within a category, or its context, to be captured when
training time series models. Combining these embeddings as features for
forecasting yields powerful models like the first place winner of the Kaggle
Taxi Trajectory Prediction111https://www.kaggle.com/c/pkdd-15-predict-taxi-
service-trajectory-i challenge (De Brébisson et al., 2015). The covariates
$\mathbf{c}_{t}$ we use are composed of time-dependent (e.g. day of week, hour
of day) and time-independent embeddings, if applicable, as well as lag
features depending on the time frequency of the data set we are training on.
All covariates are thus known for the periods we wish to forecast.

##  4 Experiments

We benchmark TimeGrad on _six_ real-world data sets and evaluate against
several competitive baselines. The source code of the model will be made
available after the review process.

###  4.1 Evaluation Metric and Data Set

For evaluation, we compute the Continuous Ranked Probability Score (CRPS)
(Matheson & Winkler, 1976) on each time series dimension, as well as on the
sum of all time series dimensions (the latter denoted by
$\mathrm{CRPS}_{\mathrm{sum}}$). CRPS measures the compatibility of a
cumulative distribution function $F$ with an observation $x$ as

$\mathrm{CRPS}(F,x)=\int_{\mathbb{R}}(F(z)-\mathbb{I}\\{x\leq
z\\})^{2}\,\mathrm{d}z,$

where $\mathbb{I}\\{x\leq z\\}$ is the indicator function which is one if
$x\leq z$ and zero otherwise. CRPS is a _proper scoring function_ , hence CRPS
attains its minimum when the predictive distribution $F$ and the data
distribution are equal. Employing the empirical CDF of $F$, i.e.
$\hat{F}(z)=\frac{1}{S}\sum_{s=1}^{S}\mathbb{I}\\{X_{s}\leq z\\}$ with $S$
samples $X_{s}\sim F$ as a natural approximation of the predictive CDF, CRPS
can be directly computed from simulated samples of the conditional
distribution (8) at each time point (Jordan et al., 2019). Finally,
$\mathrm{CRPS}_{\mathrm{sum}}$ is obtained by first summing across the $D$
time-series — both for the ground-truth data, and sampled data (yielding
$\hat{F}_{\mathrm{sum}}(t)$ for each time point). The results are then
averaged over the prediction horizon, i.e. formally
$\mathrm{CRPS}_{\mathrm{sum}}=\mathbb{E}_{t}\left[\mathrm{CRPS}\left(\hat{F}_{\mathrm{sum}}(t),\sum_{i}x_{i,t}^{0}\right)\right]$.
As proved in (de Bézenac et al., 2020) $\mathrm{CRPS}_{\mathrm{sum}}$ is also
a proper scoring function and we use it, instead of likelihood based metrics,
since not all methods we compare against yield analytical forecast
distributions or likelihoods are not meaningfully defined.

For our experiments we use Exchange (Lai et al., 2018), Solar (Lai et al.,
2018),
Electricity222https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014,
Traffic333https://archive.ics.uci.edu/ml/datasets/PEMS-SF,
Taxi444https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page and
Wikipedia555https://github.com/mbohlkeschneider/gluon-
ts/tree/mv_release/datasets open data sets, preprocessed exactly as in
(Salinas et al., 2019a), with their properties listed in Table 1. As can be
noted in the table, we do not need to normalize scales for Traffic.

Table 1: Dimension, domain, frequency, total training time steps and
prediction length properties of the training data sets used in the
experiments.
Data set |  Dim. $D$ | Dom. | Freq. |  Time steps |  Pred. steps
---|---|---|---|---|---
Exchange | $8$ | $\mathbb{R}^{+}$ | day | $6,071$ | $30$
Solar | $137$ | $\mathbb{R}^{+}$ | hour | $7,009$ | $24$
Elec. | $370$ | $\mathbb{R}^{+}$ | hour | $5,833$ | $24$
Traffic | $963$ | $(0,1)$ | hour | $4,001$ | $24$
Taxi | $1,214$ | $\mathbb{N}$ | 30-min | $1,488$ | $24$
Wiki. | $2,000$ | $\mathbb{N}$ | day | $792$ | $30$

###  4.2 Model Architecture

We train TimeGrad via SGD using Adam (Kingma & Ba, 2015) with learning rate of
$1\text{\times}{10}^{-3}$ on the training split of each data set with $N=100$
diffusion steps using a linear variance schedule starting from
$\beta_{1}=$1\text{\times}{10}^{-4}$$ till $\beta_{N}=0.1$. We construct
batches of size $64$ by taking random windows (with possible overlaps), with
the context size set to the number of prediction steps, from the total time
steps of each data set (see Table 1). For testing we use a rolling windows
prediction starting from the last context window history before the start of
the prediction and compare it to the ground-truth in the test set by sampling
$S=100$ trajectories.

The RNN consists of $2$ layers of an LSTM with the hidden state
$\mathbf{h}_{t}\in\mathbb{R}^{40}$ and we encode the noise index
$n\in\\{1,\ldots,N\\}$ using the Transformer’s (Vaswani et al., 2017) Fourier
positional embeddings, with $N_{\max}=500$, into $\mathbb{R}^{32}$ vectors.
The network $\epsilon_{\theta}$ consists of conditional 1-dim dilated ConvNets
with residual connections adapted from the WaveNet (van den Oord et al.,
2016a) and DiffWave (Kong et al., 2021) models. Figure 2 shows the schematics
of a single residual block $i=\\{0,\ldots,7\\}$ together with the final output
from the sum of all the $8$ skip-connections. All, but the last, convolutional
network layers have an output channel size of $8$ and we use a _bidirectional_
dilated convolution in each block $i$ by setting its dilation to $2^{i\%2}$.
We use a validation set from the training data of the same size as the test
set to tune the number of epochs for early stopping.

All experiments run on a single Nvidia V100 GPU with $16$GB of memory.

Figure 2: The network architecture of $\epsilon_{\theta}$ consisting of
$\mathtt{residual\\_layers}=8$ conditional residual blocks with the Gated
Activation Unit $\sigma(\cdot)\odot\tanh(\cdot)$ from (van den Oord et al.,
2016b); whose skip-connection outputs are summed up to compute the final
output. Conv1x1 and Conv1d are 1D convolutional layers with filter size of $1$
and $3$, respectively, circular padding so that the spatial size remains $D$,
and all but the last convolutional layer has output channels
$\mathtt{residual\\_channels}=8$. FC are linear layers used to up/down-sample
the input to the appropriate size for broadcasting.

###  4.3 Results

Table 2: Test set $\mathrm{CRPS}_{\mathrm{sum}}$ comparison (lower is better)
of models on six real world data sets. Mean and standard error metrics for
TimeGrad obtained by re-training and evaluating $10$ times.
Method | Exchange | Solar | Electricity | Traffic | Taxi | Wikipedia
---|---|---|---|---|---|---
VES | $\mathbf{0.005}\scriptstyle{\pm 0.000}$ | $0.9\scriptstyle{\pm 0.003}$ | $0.88\scriptstyle{\pm 0.0035}$ | $0.35\scriptstyle{\pm 0.0023}$ | - | -
VAR | $\mathbf{0.005}\scriptstyle{\pm 0.000}$ | $0.83\scriptstyle{\pm 0.006}$ | $0.039\scriptstyle{\pm 0.0005}$ | $0.29\scriptstyle{\pm 0.005}$ | - | -
VAR-Lasso | $0.012\scriptstyle{\pm 0.0002}$ | $0.51\scriptstyle{\pm 0.006}$ | $0.025\scriptstyle{\pm 0.0002}$ | $0.15\scriptstyle{\pm 0.002}$ | - | $3.1\scriptstyle{\pm 0.004}$
GARCH | $0.023\scriptstyle{\pm 0.000}$ | $0.88\scriptstyle{\pm 0.002}$ | $0.19\scriptstyle{\pm 0.001}$ | $0.37\scriptstyle{\pm 0.0016}$ | - | -
KVAE | $0.014\scriptstyle{\pm 0.002}$ | $0.34\scriptstyle{\pm 0.025}$ | $0.051\scriptstyle{\pm 0.019}$ | $0.1\scriptstyle{\pm 0.005}$ | - | $0.095\scriptstyle{\pm 0.012}$
Vec-LSTM ind-scaling | $0.008\scriptstyle{\pm 0.001}$ | $0.391\scriptstyle{\pm 0.017}$ | $0.025\scriptstyle{\pm 0.001}$ | $0.087\scriptstyle{\pm 0.041}$ | $0.506\scriptstyle{\pm 0.005}$ | $0.133\scriptstyle{\pm 0.002}$
Vec-LSTM lowrank-Copula | $0.007\scriptstyle{\pm 0.000}$ | $0.319\scriptstyle{\pm 0.011}$ | $0.064\scriptstyle{\pm 0.008}$ | $0.103\scriptstyle{\pm 0.006}$ | $0.326\scriptstyle{\pm 0.007}$ | $0.241\scriptstyle{\pm 0.033}$
GP scaling | $0.009\scriptstyle{\pm 0.000}$ | $0.368\scriptstyle{\pm 0.012}$ | $0.022\scriptstyle{\pm 0.000}$ | $0.079\scriptstyle{\pm 0.000}$ | $0.183\scriptstyle{\pm 0.395}$ | $1.483\scriptstyle{\pm 1.034}$
GP Copula | $0.007\scriptstyle{\pm 0.000}$ | $0.337\scriptstyle{\pm 0.024}$ | $0.0245\scriptstyle{\pm 0.002}$ | $0.078\scriptstyle{\pm 0.002}$ | $0.208\scriptstyle{\pm 0.183}$ | $0.086\scriptstyle{\pm 0.004}$
Transformer MAF | $\mathbf{0.005}\scriptstyle{\pm 0.003}$ | $0.301\scriptstyle{\pm 0.014}$ | $0.0207\scriptstyle{\pm 0.000}$ | $0.056\scriptstyle{\pm 0.001}$ | $0.179\scriptstyle{\pm 0.002}$ | $0.063\scriptstyle{\pm 0.003}$
TimeGrad | $0.006\scriptstyle{\pm 0.001}$ | $\mathbf{0.287}\scriptstyle{\pm 0.02}$ | $\mathbf{0.0206}\scriptstyle{\pm 0.001}$ | $\mathbf{0.044}\scriptstyle{\pm 0.006}$ | $\mathbf{0.114}\scriptstyle{\pm 0.02}$ | $\mathbf{0.0485}\scriptstyle{\pm 0.002}$

Using the $\mathrm{CRPS}_{\mathrm{sum}}$ as an evaluation metric, we compare
test time predictions of TimeGrad to a wide range of existing methods
including classical multivariate methods:

  * •

VAR (Lütkepohl, 2007) a mutlivariate linear vector auto-regressive model with
lags corresponding to the periodicity of the data,

  * •

VAR-Lasso a Lasso regularized VAR,

  * •

GARCH (van der Weide, 2002) a multivariate conditional heteroskedastic model
and

  * •

VES a innovation state space model (Hyndman et al., 2008);

as well as deep learning based methods namely:

  * •

KVAE (Fraccaro et al., 2017) a variational autoencoder to represent the data
on top of a linear state space model which describes the dynamics,

  * •

Vec-LSTM-ind-scaling (Salinas et al., 2019a) which models the dynamics via an
RNN and outputs the parameters of an _independent_ Gaussian distribution with
mean-scaling,

  * •

Vec-LSTM-lowrank-Copula (Salinas et al., 2019a) which instead parametrizes a
low-rank plus diagonal covariance via Copula process,

  * •

GP-scaling (Salinas et al., 2019a) which unrolls an LSTM with scaling on each
individual time series before reconstructing the joint distribution via a low-
rank Gaussian,

  * •

GP-Copula (Salinas et al., 2019a) which unrolls an LSTM on each individual
time series and then the joint emission distribution is given by a low-rank
plus diagonal covariance Gaussian copula and

  * •

Transformer-MAF (Rasul et al., 2021) which uses Transformer (Vaswani et al.,
2017) to model the temporal conditioning and Masked Autoregressive Flow
(Papamakarios et al., 2017) for the distribution emission model.

Table 2 lists the corresponding $\mathrm{CRPS}_{\mathrm{sum}}$ values averaged
over $10$ independent runs together with their empirical standard deviations
and shows that the TimeGrad model sets the new state-of-the-art on all but the
smallest of the benchmark data sets. Note that flow based models must apply
continuous transformations onto a continuously connected distribution, making
it difficult to model disconnected modes. Flow models assign spurious density
to connections between these modes leading to potential inaccuracies.
Similarly the generator network in variational autoencoders must learn to map
from some continuous space to a possibly disconnected space which might not be
possible to learn. In contrast EMBs do not suffer from these issues (Du &
Mordatch, 2019).

###  4.4 Ablation

The length $N$ of the forward process is a crucial hyperparameter, as a bigger
$N$ allows the reverse process to be approximately Gaussian (Sohl-Dickstein et
al., 2015) which assists the Gaussian parametrization (1) to approximate it
better. We evaluate to which extent, if any at all, larger $N$ affects
prediction performance, with an ablation study where we record the test set
$\mathrm{CRPS}_{\mathrm{sum}}$ of the Electricity data set for different total
diffusion process lengths $N=2,4,8,\ldots,256$ while keeping all other
hyperparemeters unchanged. The results are then plotted in Figure 3 where we
note that $N$ can be reduced down to $\approx 10$ without significant
performance loss. An optimal value is achieved at $N\approx 100$ and larger
levels are not beneficial if all else is kept fixed.

Figure 3: TimeGrad test set $\mathrm{CRPS}_{\mathrm{sum}}$ for Electricity
data by varying total diffusion length $N$. Good performance is established
already at $N\approx 10$ with optimal value at $N\approx 100$. The mean and
standard errors obtained over $5$ independent runs. We see similar behaviour
with other data sets.

To highlight the predictions of TimeGrad we show in Figure 4 the predicted
median, $50\%$ and $90\%$ distribution intervals of the first $6$ dimensions
of the full $963$ dimensional multivariate forecast of the Traffic benchmark.

Figure 4: TimeGrad prediction intervals and test set ground-truth for Traffic
data of the first $6$ of $963$ dimensions from first rolling-window. Note that
neighboring entities have an order of magnitude difference in scales.

##  5 Related Work

###  5.1 Energy-Based Methods

The EBM of (Ho et al., 2020) that we adapt is based on methods that learn the
gradient of the log-density with respect to the _inputs_ , called Stein Score
function (Hyvärinen, 2005; Vincent, 2011), and at inference time use this
gradient estimate via Langevin dynamics to sample from the model of this
complicated data distribution (Song & Ermon, 2019). These models achieve
impressive results for image generation (Ho et al., 2020; Song & Ermon, 2020)
when trained in an unsupervised fashion without requiring adversarial
optimization. By perturbing the data using multiple noise scales, the learnt
Score network captures both coarse and fine-grained data features.

The closest related work to TimeGrad is in the recent non-autoregressive
conditional methods for high fidelity waveform generation (Chen et al., 2021;
Kong et al., 2021). Although these methods learn the distribution of vector
valued data via denoising diffusion methods, as done here, they do not
consider its temporal development. Also neighboring dimensions of waveform
data are highly correlated and have a uniform scale, which is not necessarily
true for multivariate time series problems where neighboring entities occur
arbitrarily (but in a fixed order) and can have different scales. (Du &
Mordatch, 2019) also use EBMs to model one and multiple steps for a trajectory
modeling task in an non-autoregressive fashion.

###  5.2 Time Series Forecasting

Neural time series methods have recently become popular ways of solving the
prediction problem via univariate point forecasting methods (Oreshkin et al.,
2020; Smyl, 2020) or univariate probabilistic methods (Salinas et al., 2019b).
In the multivariate setting we also have point forecasting methods (Lai et
al., 2018; Li et al., 2019) as well as probabilistic methods, like this
method, which explicitly model the data distribution using Gaussian copulas
(Salinas et al., 2019a), GANs (Yoon et al., 2019), or normalizing flows (de
Bézenac et al., 2020; Rasul et al., 2021). Bayesian neural networks can also
be used to provide _epistemic_ uncertainty in forecasts as well as detect
distributional shifts (Zhu & Laptev, 2018), although these methods often do
not perform as well empirically (Wenzel et al., 2020).

##  6 Conclusion and Future Work

We have presented TimeGrad, a versatile multivariate probabilistic time series
forecasting method that leverages the exceptional performance of EBMs to learn
and sample from the distribution of the next time step, autoregressivly.
Analysis of TimeGrad on six commonly used time series benchmarks establishes
the new state-of-the-art against competitive methods.

We note that while training TimeGrad we do not need to loop over the EBM
function approximator $\epsilon_{\theta}$, unlike in the normalizing flow
setting where we have multiple stacks of bijections. However while sampling we
do loop $N$ times over $\epsilon_{\theta}$. A possible strategy to improve
sampling times introduced in (Chen et al., 2021) uses a combination of
improved variance schedule and an $L_{1}$ loss to allow sampling with fewer
steps at the cost of a small reduction in quality if such a trade-off is
required. A recent paper (Song et al., 2021) generalize the diffusion
processes via a class of non-Markovian processes which also allows for faster
sampling.

The use of normalizing flows for discrete valued data dictates that one
dequantizes it (Theis et al., 2016), by adding uniform noise to the data,
before using the flows to learn. Dequantization is not needed in the EBM
setting and future work could explore methods of explicitly modeling discrete
distributions.

As noted in (Du & Mordatch, 2019) EBMs exhibit better out-of-distribution
(OOD) detection than other likelihood models. Such a task requires models to
have a high likelihood on the data manifold and low at all other locations.
Surprisingly (Nalisnick et al., 2019) showed that likelihood models, including
flows, were assigning higher likelihoods to OOD data whereas EBMs do not
suffer from this issue since they penalize high probability under the model
but low probability under the data distribution explicitly. Future work could
evaluate the usage of TimeGrad for anomaly detection tasks.

For long time sequences, one could replace the RNN with a Transformer
architecture (Rasul et al., 2021) to provide better conditioning for the EBM
emission head. Concurrently, since EBMs are not constrained by the form of
their functional approximators, one natural way to improve the model would be
to incorporate architectural choices that best encode the inductive bias of
the problem being tackled, for example with graph neural networks (Niu et al.,
2020) when the relationships between entities are known.

## References

  * Benidis et al. (2020) Benidis, K., Rangapuram, S. S., Flunkert, V., Wang, B., Maddix, D., Turkmen, C., Gasthaus, J., Bohlke-Schneider, M., Salinas, D., Stella, L., Callot, L., and Januschowski, T.  Neural forecasting: Introduction and literature overview, 2020.
  * Charrington (2018) Charrington, S.  TWiML & AI Podcast: Systems and Software for Machine Learning at Scale with Jeff Dean, 2018.  URL https://bit.ly/2G0LmGg.
  * Chen et al. (2021) Chen, N., Zhang, Y., Zen, H., Weiss, R. J., Norouzi, M., and Chan, W.  WaveGrad: Estimating gradients for waveform generation.  In _International Conference on Learning Representations 2021 (Conference Track)_ , 2021.  URL https://openreview.net/forum?id=NsMLjcFaO8O.
  * Chung et al. (2014) Chung, J., Gulcehre, C., Cho, K., and Bengio, Y.  Empirical evaluation of gated recurrent neural networks on sequence modeling.  In _NIPS 2014 Workshop on Deep Learning, December 2014_ , 2014.
  * de Bézenac et al. (2020) de Bézenac, E., Rangapuram, S. S., Benidis, K., Bohlke-Schneider, M., Kurle, R., Stella, L., Hasson, H., Gallinari, P., and Januschowski, T.  Normalizing Kalman Filters for Multivariate Time series Analysis.  In _Advances in Neural Information Processing Systems_ , volume 33. Curran Associates, Inc., 2020.
  * De Brébisson et al. (2015) De Brébisson, A., Simon, E., Auvolat, A., Vincent, P., and Bengio, Y.  Artificial Neural Networks Applied to Taxi Destination Prediction.  In _Proceedings of the 2015th International Conference on ECML PKDD Discovery Challenge - Volume 1526_ , ECMLPKDDDC’15, pp. 40–51, Aachen, Germany, Germany, 2015. CEUR-WS.org.  URL http://dl.acm.org/citation.cfm?id=3056172.3056178.
  * Dinh et al. (2017) Dinh, L., Sohl-Dickstein, J., and Bengio, S.  Density estimation using Real NVP.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net, 2017.  URL https://openreview.net/forum?id=HkpbnH9lx.
  * Du & Mordatch (2019) Du, Y. and Mordatch, I.  Implicit Generation and Modeling with Energy Based Models.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32, pp. 3608–3618. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/378a063b8fdb1db941e34f4bde584c7d-Paper.pdf.
  * Fraccaro et al. (2017) Fraccaro, M., Kamronn, S., Paquet, U., and Winther, O.  A Disentangled Recognition and Nonlinear Dynamics Model for Unsupervised Learning.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 30, pp. 3601–3610. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper/2017/file/7b7a53e239400a13bd6be6c91c4f6c4e-Paper.pdf.
  * Graves (2013) Graves, A.  Generating Sequences With Recurrent Neural Networks.  _arXiv preprint arXiv:1308.0850_ , 2013.
  * Hinton (2002) Hinton, G. E.  Training Products of Experts by Minimizing Contrastive Divergence.  _Neural Computation_ , 14(8):1771––1800, August 2002.  ISSN 0899-7667.  doi: 10.1162/089976602760128018.  URL https://doi.org/10.1162/089976602760128018.
  * Ho et al. (2020) Ho, J., Jain, A., and Abbeel, P.  Denoising Diffusion Probabilistic Models.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 33. Curran Associates, Inc., 2020.  URL https://papers.nips.cc/paper/2020/file/4c5bcfec8584af0d967f1ab10179ca4b-Paper.pdf.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long Short-Term Memory.  _Neural Computation_ , 9(8):1735–1780, November 1997.  ISSN 0899-7667.  doi: 10.1162/neco.1997.9.8.1735.
  * Hyndman & Athanasopoulos (2018) Hyndman, R. and Athanasopoulos, G.  _Forecasting: Principles and practice_.  OTexts, 2018.  ISBN 9780987507112.
  * Hyndman et al. (2008) Hyndman, R., Koehler, A., Ord, K., and Snyder, R.  _Forecasting with exponential smoothing. The state space approach_ , chapter 17, pp. 287–300.  Springer-Verlag, 2008.  doi: 10.1007/978-3-540-71918-2.
  * Hyvärinen (2005) Hyvärinen, A.  Estimation of Non-Normalized Statistical Models by Score Matching.  _Journal of Machine Learning Research_ , 6(24):695–709, 2005.  URL http://jmlr.org/papers/v6/hyvarinen05a.html.
  * Jordan et al. (2019) Jordan, A., Krüger, F., and Lerch, S.  Evaluating Probabilistic Forecasts with scoringRules.  _Journal of Statistical Software, Articles_ , 90(12):1–37, 2019.  ISSN 1548-7660.  doi: 10.18637/jss.v090.i12.  URL https://www.jstatsoft.org/v090/i12.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In _International Conference on Learning Representations (ICLR)_ , 2015.
  * Kingma & Welling (2019) Kingma, D. P. and Welling, M.  An Introduction to Variational Autoencoders.  _Foundations and Trends in Machine Learning_ , 12(4):307–392, 2019.  doi: 10.1561/2200000056.  URL https://doi.org/10.1561/2200000056.
  * Kong et al. (2021) Kong, Z., Ping, W., Huang, J., Zhao, K., and Catanzaro, B.  DiffWave: A Versatile Diffusion Model for Audio Synthesis.  In _International Conference on Learning Representations 2021 (Conference Track)_ , 2021.  URL https://openreview.net/forum?id=a-xFK8Ymz5J.
  * Lai et al. (2018) Lai, G., Chang, W.-C., Yang, Y., and Liu, H.  Modeling Long- and Short-Term Temporal Patterns with Deep Neural Networks.  In _The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval_, SIGIR ’18, pp. 95–104, New York, NY, USA, 2018. ACM.  ISBN 978-1-4503-5657-2.  doi: 10.1145/3209978.3210006.  URL http://doi.acm.org/10.1145/3209978.3210006.
  * LeCun et al. (2006) LeCun, Y., Chopra, S., Hadsell, R., Ranzato, M., and Huang, F.  A Tutorial on Energy-Based Learning.  In Bakir, G., Hofman, T., Schölkopf, B., Smola, A., and Taskar, B. (eds.), _Predicting Structured Data_. MIT Press, 2006.
  * Li et al. (2019) Li, S., Jin, X., Xuan, Y., Zhou, X., Chen, W., Wang, Y.-X., and Yan, X.  Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting.  In Wallach, H., Larochelle, H., Beygelzimer, A., d’Alché Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 5244–5254. Curran Associates, Inc., 2019.
  * Lütkepohl (2007) Lütkepohl, H.  _New Introduction to Multiple Time Series Analysis_.  Springer Berlin Heidelberg, 2007.  ISBN 9783540262398.  URL https://books.google.de/books?id=muorJ6FHIiEC.
  * Matheson & Winkler (1976) Matheson, J. E. and Winkler, R. L.  Scoring Rules for Continuous Probability Distributions.  _Management Science_ , 22(10):1087–1096, 1976\.
  * Nalisnick et al. (2019) Nalisnick, E., Matsukawa, A., Teh, Y. W., Gorur, D., and Lakshminarayanan, B.  Do Deep Generative Models Know What They Don’t Know?  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=H1xwNhCcYm.
  * Niu et al. (2020) Niu, C., Song, Y., Song, J., Zhao, S., Grover, A., and Ermon, S.  Permutation Invariant Graph Generation via Score-Based Generative Modeling.  In Chiappa, S. and Calandra, R. (eds.), _The 23rd International Conference on Artificial Intelligence and Statistics, AISTATS 2020, 26-28 August 2020, Online [Palermo, Sicily, Italy]_ , volume 108 of _Proceedings of Machine Learning Research_ , pp. 4474–4484. PMLR, 2020\.
  * Oreshkin et al. (2020) Oreshkin, B. N., Carpov, D., Chapados, N., and Bengio, Y.  N-BEATS: Neural basis expansion analysis for interpretable time series forecasting.  In _International Conference on Learning Representations_ , 2020.  URL https://openreview.net/forum?id=r1ecqn4YwB.
  * Papamakarios et al. (2017) Papamakarios, G., Pavlakou, T., and Murray, I.  Masked Autoregressive Flow for Density Estimation.  _Advances in Neural Information Processing Systems 30_ , 2017.
  * Papamakarios et al. (2019) Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan, B.  Normalizing Flows for Probabilistic Modeling and Inference, 2019.
  * Rasul et al. (2021) Rasul, K., Sheikh, A.-S., Schuster, I., Bergmann, U., and Vollgraf, R.  Multivariate Probabilistic Time Series Forecasting via Conditioned Normalizing Flows.  In _International Conference on Learning Representations 2021 (Conference Track)_ , 2021.  URL https://openreview.net/forum?id=WiGQBFuVRv.
  * Salinas et al. (2019a) Salinas, D., Bohlke-Schneider, M., Callot, L., Medico, R., and Gasthaus, J.  High-dimensional multivariate forecasting with low-rank Gaussian Copula Processes.  In Wallach, H., Larochelle, H., Beygelzimer, A., d’Alché Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 6824–6834. Curran Associates, Inc., 2019a.
  * Salinas et al. (2019b) Salinas, D., Flunkert, V., Gasthaus, J., and Januschowski, T.  DeepAR: Probabilistic forecasting with autoregressive recurrent networks.  _International Journal of Forecasting_ , 2019b.  ISSN 0169-2070.  URL http://www.sciencedirect.com/science/article/pii/S0169207019301888.
  * Smyl (2020) Smyl, S.  A hybrid method of exponential smoothing and recurrent neural networks for time series forecasting.  _International Journal of Forecasting_ , 36(1):75–85, 2020.  ISSN 0169-2070.  doi: https://doi.org/10.1016/j.ijforecast.2019.03.017.  URL http://www.sciencedirect.com/science/article/pii/S0169207019301153.  M4 Competition.
  * Sohl-Dickstein et al. (2015) Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., and Ganguli, S.  Deep Unsupervised Learning using Nonequilibrium Thermodynamics.  In Bach, F. and Blei, D. (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 2256–2265, Lille, France, 2015. PMLR.  URL http://proceedings.mlr.press/v37/sohl-dickstein15.html.
  * Song et al. (2021) Song, J., Meng, C., and Ermon, S.  Denoising Diffusion Implicit Models.  In _International Conference on Learning Representations 2021 (Conference Track)_ , 2021.  URL https://openreview.net/pdf?id=St1giarCHLP.
  * Song & Ermon (2019) Song, Y. and Ermon, S.  Generative Modeling by Estimating Gradients of the Data Distribution.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32, pp. 11918–11930. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/3001ef257407d5a371a96dcd947c7d93-Paper.pdf.
  * Song & Ermon (2020) Song, Y. and Ermon, S.  Improved Techniques for Training Score-Based Generative Models.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 33. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/92c3b916311a5517d9290576e3ea37ad-Paper.pdf.
  * Song & Kingma (2021) Song, Y. and Kingma, D. P.  How to Train Your Energy-Based Models.  2021\.  URL https://arxiv.org/abs/2101.03288.
  * Sutskever et al. (2014) Sutskever, I., Vinyals, O., and Le, Q. V.  Sequence to Sequence Learning with Neural Networks.  In Ghahramani, Z., Welling, M., Cortes, C., Lawrence, N., and Weinberger, K. (eds.), _Advances in Neural Information Processing Systems 27_ , pp. 3104–3112. Curran Associates, Inc., 2014.
  * Theis et al. (2016) Theis, L., van den Oord, A., and Bethge, M.  A note on the evaluation of generative models.  In _International Conference on Learning Representations_ , 2016.  URL http://arxiv.org/abs/1511.01844.  arXiv:1511.01844.
  * Tsay (2014) Tsay, R. S.  _Multivariate Time Series Analysis: With R and Financial Applications_.  Wiley Series in Probability and Statistics. Wiley, 2014.  ISBN 9781118617908.
  * van den Oord et al. (2016a) van den Oord, A., Dieleman, S., Zen, H., Simonyan, K., Vinyals, O., Graves, A., Kalchbrenner, N., Senior, A., and Kavukcuoglu, K.  WaveNet: A Generative Model for Raw Audio.  In _The 9th ISCA Speech Synthesis Workshop, Sunnyvale, CA, USA, 13-15 September 2016_ , pp. 125. ISCA, 2016a.  URL http://www.isca-speech.org/archive/SSW_2016/abstracts/ssw9_DS-4_van_den_Oord.html.
  * van den Oord et al. (2016b) van den Oord, A., Kalchbrenner, N., Espeholt, L., kavukcuoglu, k., Vinyals, O., and Graves, A.  Conditional Image Generation with PixelCNN Decoders.  In Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 29, pp. 4790–4798. Curran Associates, Inc., 2016b.  URL https://proceedings.neurips.cc/paper/2016/file/b1301141feffabac455e1f90a7de2054-Paper.pdf.
  * van den Oord et al. (2016c) van den Oord, A., Kalchbrenner, N., and Kavukcuoglu, K.  Pixel Recurrent Neural Networks.  In Balcan, M. F. and Weinberger, K. Q. (eds.), _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48 of _Proceedings of Machine Learning Research_ , pp. 1747–1756, New York, New York, USA, 20–22 Jun 2016c. PMLR.  URL http://proceedings.mlr.press/v48/oord16.html.
  * van der Weide (2002) van der Weide, R.  GO-GARCH: a multivariate generalized orthogonal GARCH model.  _Journal of Applied Econometrics_ , 17(5):549–564, 2002.  doi: 10.1002/jae.688.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L. u., and Polosukhin, I.  Attention is All you Need.  In Guyon, I., Luxburg, U., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 30_ , pp. 5998–6008. Curran Associates, Inc., 2017.  URL http://papers.nips.cc/paper/7181-attention-is-all-you-need.pdf.
  * Vincent (2011) Vincent, P.  A Connection Between Score Matching and Denoising Autoencoders.  _Neural Computation_ , 23(7):1661–1674, 2011\.  URL https://doi.org/10.1162/NECO_a_00142.
  * Wenzel et al. (2020) Wenzel, F., Roth, K., Veeling, B., Swiatkowski, J., Tran, L., Mandt, S., Snoek, J., Salimans, T., Jenatton, R., and Nowozin, S.  How good is the Bayes posterior in deep neural networks really?  In III, H. D. and Singh, A. (eds.), _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pp. 10248–10259. PMLR, 13–18 Jul 2020.  URL http://proceedings.mlr.press/v119/wenzel20a.html.
  * Yoon et al. (2019) Yoon, J., Jarrett, D., and van der Schaar, M.  Time-series Generative Adversarial Networks.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32, pp. 5508–5518. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/c9efe5f26cd17ba6216bbe2a7d26d490-Paper.pdf.
  * Zhu & Laptev (2018) Zhu, L. and Laptev, N.  Deep and Confident Prediction for Time Series at Uber.  In _2017 IEEE International Conference on Data Mining Workshops (ICDMW)_ , volume 00, pp. 103–110, November 2018.  doi: 10.1109/ICDMW.2017.19.  URL doi.ieeecomputersociety.org/10.1109/ICDMW.2017.19.