The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
80,872 characters
From Vector Autoregressions to AI-based Time Series Forecasting: A Review
\maketitle
\begin{abstract}
Forecasting is a central goal of time-series analysis. This review centers on
three major developments in recent AI-based time-series
forecasting: transformers, large
pretrained models for zero-shot forecasting, and diffusion-based generative
forecasters. We connect these methods to the econometric tradition built around
the vector autoregression (VAR) through a common object: the conditional
distribution of the future given the past. The review is organized around three
long-standing challenges: \emph{high dimensionality}, \emph{nonstationarity},
and \emph{nonlinearity}. We argue that modern methods make progress by
expanding the classical forecasting template: they allow more flexible dynamics,
use larger information sets and training corpora, and represent richer
predictive distributions. Yet they often lack the inferential and structural
tools that make classical models useful for testing, explanation, and policy
analysis. We close by outlining open problems where econometric tools remain
important.
\medskip
\noindent\textbf{Keywords:} vector autoregression; transformers; foundation models;
denoising diffusion; probabilistic forecasting; high-dimensional time series.
\end{abstract}
\section{Introduction}
\label{sec:intro}
Forecasting is central to applied econometrics, with uses ranging from monetary
policy analysis and risk management to demand planning and nowcasting. For four
decades, the vector autoregression (VAR) of \citet{sims1980macroeconomics} has
provided a coherent way to study how multiple time series move together. Its
appeal comes not only from forecasting performance, but also from the
interpretive apparatus built around it: impulse-response analysis, Granger
causality \citep{granger1969investigating}, and the cointegration and
error-correction framework of \citet{engle1987cointegration} and
\citet{johansen1991estimation}. Standard graduate treatments
\citep{lutkepohl2005new,hamilton1994time} still organize the field around this
core. For many problems, especially small systems and settings where
interpretation, valid inference, or long-run relationships matter, the VAR
remains the natural starting point.
At the same time, methods imported from machine learning have expanded what is
possible at scale. The transformer architecture \citep{vaswani2017attention} is
a leading example. For time-series forecasting, it can be read as a nonlinear
generalization of the VAR: fixed lag coefficients are replaced by
data-dependent attention weights. In its standard forecasting form, however, it
usually produces point forecasts or marginal quantiles rather than a full
econometric model. Foundation models
\citep{ansari2024chronos,das2024timesfm,woo2024moirai} push this logic further
by pretraining once on a large and diverse collection of time series, then
forecasting new series without task-specific estimation. Denoising diffusion
models \citep{ho2020denoising,rasul2021timegrad} add a generative layer,
sampling future paths from the predictive distribution and thereby representing
forecast uncertainty more directly. Although these developments are often
presented as a break from classical statistics, we argue that they are better
understood as new approaches to the same forecasting problem.
Our aim is not simply to translate between two literatures. The central point is
an asymmetry: modern methods often forecast well, but offer less support for
explanation. They do not naturally produce impulse responses, testable
restrictions, long-run equilibria, or classical sampling uncertainty.
Identification, valid inference, structural interpretation, and
nonstationarity therefore remain open problems. We see these gaps not only as
weaknesses, but also as opportunities. They are precisely the kinds of problems
econometrics is designed to address. Tools such as the bootstrap, conformal
prediction, cointegration, and change-point methods may help connect modern
predictive models with econometric interpretation. The review is written with
this goal in mind: each section places the methods against three recurring
challenges, and the conclusion turns the remaining gaps into a research agenda.
As a running example, consider a small macroeconomic panel containing inflation,
output growth, and a short-term interest rate. A classical VAR uses their joint
history to forecast future values. An identified Structural VAR asks how a monetary policy
shock propagates through inflation and output. A Vector Error Correction Model (VECM) is appropriate when the
variables share a long-run equilibrium. A transformer replaces fixed lag
coefficients with data-dependent weights. A diffusion forecaster describes
uncertainty by sampling possible future macroeconomic paths. We return to this
example throughout to distinguish forecasting performance from the econometric
questions of identification, inference, and nonstationarity.
\paragraph{A unifying view.}
Every model in this review can be read as an answer to a single question: what is the conditional law
of the future $\bm{y}_{t+1:t+H}$ given the past $\bm{y}_{1:t}$ and any covariates? A Gaussian VAR
answers with fixed linear coefficients and a Gaussian predictive density. A transformer answers with
data-dependent nonlinear attention weights, but typically returns only point or quantile summaries.
A diffusion model answers by learning to sample from the predictive law without explicitly
parameterising its density. The differences among these methods are differences in flexibility: they
represent the same object, the conditional distribution, in different ways. Holding this object fixed
clarifies both what is gained as we move from VARs to attention and diffusion, and what is lost:
closed-form inference, orthogonalised shocks, and testable restrictions.
\paragraph{Three challenges.}
The classical linear Gaussian model makes three assumptions that real data often violate. Much of the
history of time series analysis, classical and modern alike, can be read as the progressive relaxation
of these assumptions.
\begin{enumerate}
\item \textbf{High dimensionality.} A VAR$(L)$ on $p$ series has $O(p^2L)$ parameters, so unrestricted estimation
quickly becomes unreliable as $p$ grows.
Classical responses include factor structure, shrinkage, and sparsity through
Lasso-type penalization. Each imposes an explicit restriction before
estimation. Modern architectures address the same problem through channel
independence, cross-variable attention, and pretraining on large collections
of series from many domains. Their distinctive feature is adaptivity: the
restriction is less often specified in advance and more often learned from
data. In this sense, the dimensionality problem is partly relocated from
task-specific estimation to pooled representation learning. A full theoretical
account of when this adaptivity succeeds, and when it fails, is still missing.
\item \textbf{Nonstationarity.} Means, variances, and dependence structures
drift over time, and series may contain unit roots, share long-run
equilibria, or move across regimes. Classical time-series analysis matches
these forms of instability with specific tools: unit-root testing,
differencing, and cointegration for stochastic trends, change-point methods
for abrupt breaks, and time-varying coefficients for gradual drift. The modern
forecasting literature treats nonstationarity more indirectly, relying mainly
on instance normalization, scaling, and conditioning on recent history.
\item \textbf{Nonlinearity.} Linear fixed-coefficient dynamics cannot capture regime-dependent
propagation, asymmetric responses, or non-Gaussian predictive shapes. Machine-learning models can
instead be viewed as nonparametric methods for learning nonlinear dependence. Attention and
diffusion learn nonlinear and non-Gaussian conditional laws without imposing a parametric form in
advance. Attention is especially distinctive because it can compare tokens at arbitrary distances:
classical models typically discount the remote past by construction, while attention weights past
dates according to learned relevance. Distance and dependence are therefore decoupled.
\end{enumerate}
These three challenges organize the rest of the review. Section~\ref{sec:classical}
reviews the classical foundations and locates the three challenges within them.
Section~\ref{sec:transformers} develops the transformer as a data-dependent
generalization of the vector autoregression, and especially
Subsection~\ref{sec:tokenize} focuses on the representational
question inherited from language modelling: what counts as a token of a time
series? Section~\ref{sec:fm} reviews foundation models and zero-shot
forecasting, emphasizing the shift from per-data-set estimation to cross-corpus
pretraining. Section~\ref{sec:diffusion} reviews diffusion-based generative
forecasters, emphasizing the move from marginal summaries to sampled joint
forecast paths.
Section~\ref{sec:synthesis} offers a synthesis and summary. Throughout, we aim
less at exhaustive coverage than at a coherent narrative connecting modern
forecasting methods to econometric concerns. Recent surveys provide broader
coverage of transformer methods \citep{wen2023transformers} and diffusion
models \citep{yang2024survey,meijer2024rise}.
\paragraph{Notation and terminology.}
Throughout, $\bm{y}_t\in\mathbb{R}^{p}$ denotes the vector of $p$ series observed at date $t$, and $T$
denotes the sample length. The symbol $L$ denotes the length of the look-back window. In
Section~\ref{sec:classical}, this is the lag order of the vector autoregression. From
Section~\ref{sec:transformers} onward, it is the context length of the sequence model. In both cases,
the input is a window of the $L$ most recent observations. The symbol $H$ denotes the forecast
horizon, so that $\bm{y}_{t+1:t+H}$ is the object to be predicted, and $h$ indexes individual
forecast steps.
\section{Classical multivariate foundations}
\label{sec:classical}
For a $p\times 1$ vector $\bm{y}_t$, the vector autoregression of order $L$ is
\begin{equation}
\bm{y}_t = \bm{\nu} + \sum_{i=1}^{L} A_i \bm{y}_{t-i} + \bm{\varepsilon}_t,
\qquad \bm{\varepsilon}_t \sim \mathrm{iid}\,(\bm{0}, W),
\label{eq:var}
\end{equation}
where $A_i \in \mathbb{R}^{p\times p}$ are coefficient matrices and $W$ is the innovation covariance matrix,
which is generally non-diagonal. Because every equation shares the same regressors, the system is a
seemingly unrelated regression in which equation-by-equation ordinary least squares is efficient, and
a complete distribution theory is available \citep{lutkepohl2005new,hamilton1994time}. When the roots
of $\det(I_p - A_1 z - \cdots - A_L z^L)$ lie outside the unit circle, the process is stable and
admits the Wold representation
\[
\bm{y}_t = \bm{\mu} + \sum_{h\ge 0} B_h \bm{\varepsilon}_{t-h},
\]
with $B_0 = I_p$. The matrices $B_h$ are the basis of structural analysis. Augmenting \eqref{eq:var} with moving-average terms gives the vector
autoregressive moving-average (VARMA) class. Under the same stability condition,
a VARMA model also admits an infinite moving-average representation. It can
therefore represent current outcomes as weighted accumulations of past shocks,
often more parsimoniously than a finite VAR, although estimation and
identification are harder \citep{lutkepohl2005new}. This shock-accumulation
form provides a useful, though only formal, parallel with the accumulation of
artificial noise in diffusion models (Section~\ref{sec:diffusion}).
In the macroeconomic panel example, $\bm y_t$ might collect inflation, output
growth, and the policy rate. The coefficient matrix $A_i$ records how the
$i$th lag of each variable enters each equation. Lagged inflation may help
forecast the policy rate, lagged policy rates may help forecast output growth,
and so on. This is a reduced-form predictive object. It becomes a structural
object only after additional identifying assumptions are imposed.
Beyond producing forecasts, a fitted VAR is used to answer two further
questions, both directly related to what modern methods do and do not deliver.
The first is \emph{propagation}: what happens after an innovation or shock
enters the system? This leads to impulse-response analysis. The second is
\emph{predictive relevance}: does the history of one series improve forecasts
of another? This leads to Granger-causality testing.
For propagation, the Wold representation gives the basic reduced-form object.
The $(i,j)$ entry of $B_h$ records the response of variable $i$, $h$ periods
after a one-unit innovation to variable $j$. Plotting these entries against the
horizon $h$ gives an impulse-response function. A complication arises because
$W$ is generally not diagonal. The reduced-form errors are correlated across
equations, so an innovation to one variable alone is not yet a well-defined
structural experiment. To interpret an impulse response as the effect of an
economically meaningful shock, such as a monetary policy shock, the correlated
innovations must first be transformed into orthogonal structural shocks. This is
the identification problem. Standard solutions impose additional assumptions,
such as a recursive ordering of contemporaneous effects or restrictions on the
signs or long-run impacts of particular shocks \citep{blanchard1989,kilian2017}.
In the running macroeconomic example, an impulse response might trace the effect
of a monetary policy shock on output growth and inflation over subsequent
quarters. The important point is that the shock is not simply an innovation in
the interest-rate equation. Because reduced-form VAR innovations are correlated,
the monetary policy shock must be identified by restrictions that give it
structural meaning.
For predictive relevance, \citet{granger1969investigating} turned the question
into a testable hypothesis. Let $y_{k,t}$ denote the $k$th component of
$\bm y_t$. In the equation for $y_{k,t}$, the statement that $y_{\ell,t}$ fails
to Granger-cause $y_{k,t}$ means that all coefficients on the lags of
$y_{\ell,t}$ are zero. This is a restriction on the matrices
$A_1,\dots,A_L$, and it can be tested by a standard $F$ or Wald statistic.
For the same panel, the statement that the policy rate Granger-causes inflation
means only that lagged policy rates improve forecasts of inflation, conditional
on the other included lags. It does not mean that changing the policy rate would
have the same effect as a policy intervention. This distinction between
predictive content and structural causation must be preserved when interpreting
modern forecasting models.
Impulse responses and Granger tests share a common requirement: fixed
coefficients. The impulse responses $B_h$ are functions of the coefficient
matrices $A_i$, and the Granger test is a hypothesis about whether particular
entries of those matrices are zero. If the coefficients vary across time or
across input windows, there is no single impulse-response function to compute
and no fixed coefficient on which to impose a zero restriction. This is exactly
the constancy that self-attention discards. Fixed coefficients are replaced by
weights recomputed from each input window. Figure~\ref{fig:var2attn} in
Section~\ref{sec:transformers} illustrates the contrast. Modern methods
therefore gain forecasting flexibility, but lose the fixed objects that make
explanation and testing straightforward.
Nonstationarity is treated directly in classical time-series analysis. The
classical toolkit matches the type of instability to the model. For stochastic trends, many economic series are integrated of order one.
Regressions among unrelated integrated series may therefore be spurious
\citep{GRANGER1974111}. Unit-root tests
\citep{dickey1979,phillips1988} help diagnose this problem. The constructive
response is cointegration: $I(1)$ series may share a stationary long-run
combination, in which case the corresponding representation is the vector
error-correction form
\begin{equation}
\Delta \bm{y}_t = \bm{\alpha}\,\bm{\beta}^{\top}\bm{y}_{t-1}
+ \sum_{j=1}^{L-1} \Gamma_j \,\Delta\bm{y}_{t-j} + \bm{\varepsilon}_t,
\label{eq:vecm}
\end{equation}
where $\Delta$ denotes the first difference, $\bm{\beta}\in\mathbb{R}^{p\times r}$
collects the $r$ cointegrating vectors defining the stationary combinations
$\bm{\beta}^{\top}\bm{y}_t$, $\bm{\alpha}\in\mathbb{R}^{p\times r}$ gives the speeds at
which deviations from those combinations are corrected, and the $\Gamma_j$ are
short-run dynamic coefficients. The cointegrating rank $r$ and the coefficient
matrices can be estimated by the reduced-rank procedure of
\citet{johansen1991estimation,johansen1995}, building on
\citet{engle1987cointegration} and the unit-root distribution theory of
\citet{dickey1979,phillips1988}.
Abrupt instability is handled differently. Structural change becomes an estimation problem: the
number and locations of change points are parameters of the process. Segmentation estimators can be
consistent with near-optimal localization rates even when the number of breaks grows with the sample;
wild binary segmentation is one example \citep{fryzlewicz2014wbs}. More recently, neural networks
have also been trained as change-point detectors with accompanying theory
\citep[for example][]{li2024automatic}. For gradual instability, parameters may instead be allowed
to drift smoothly. This is the idea behind the locally stationary framework of
\citet{dahlhaus1997fitting}, where the process is approximately stationary within short windows and
can be estimated by kernel or local-likelihood methods. Time-varying parameter models in
econometrics follow the same logic.
The practical lesson recurs below. If the data are cointegrated, a model that
uses differences alone is misspecified because it discards the
equilibrium-restoring term in \eqref{eq:vecm}. Neural and generative forecasters
can effectively make this move when long-run structure is not built into the
model. In the macroeconomic example, cointegration would matter if some
variables shared a stable long-run relation, so that deviations from that
relation predict later adjustment. A model estimated only in differences may
forecast short-run movements well while discarding the term that gives the
system its long-run discipline.
\paragraph{The three challenges.}
High dimensionality is already binding in the VAR, since the unrestricted system \eqref{eq:var}
contains $O(p^2L)$ parameters. Classical remedies are explicit and interpretable: shrinkage through
Bayesian priors of the Minnesota family \citep{litterman1986forecasting,banbura2010large};
dimension reduction through dynamic factor models \citep{stock2002forecasting,baing2002}; and,
more recently, sparsity through regularised high-dimensional VARs \citep{basu2015,kock2015}. Each
makes estimation feasible by imposing a stated restriction, in contrast to the implicit,
data-driven regularization of the methods discussed below.
Nonstationarity, as we have just seen, is handled by tools matched to the type of instability:
cointegration for stochastic trends, change-point estimation for abrupt breaks, and smoothly varying
coefficients for gradual drift.
Nonlinearity remains the main limitation. The VAR in \eqref{eq:var} is linear and has constant
coefficients, so regime-dependent propagation, asymmetric responses, and non-Gaussian predictive
shapes lie outside its scope. High dimensionality and nonlinearity are natural strengths of modern machine-learning methods, which
can scale to large panels and learn flexible dependence structures. Nonstationarity remains the axis
on which these methods have the most to learn from econometrics, with its established tools for unit
roots, cointegration, structural breaks, and parameter drift.
\section{The transformer era}
\label{sec:transformers}
Many nonlinear architectures have been brought to time-series prediction,
including recurrent networks, temporal convolutions, and deep basis-expansion
models. These methods relax the linearity of \eqref{eq:var} by replacing fixed
linear dynamics with nonlinear function approximators. Earlier global models had
already moved beyond per-series estimation. DeepAR \citep{salinas2020deepar}
used a recurrent network to produce an autoregressive predictive distribution.
N-BEATS \citep{oreshkin2020nbeats} used stacked residual blocks with trend and
seasonal basis expansions, and N-HiTS \citep{challu2023nhits} extended this idea
with multi-rate pooling and hierarchical interpolation.
This section focuses on the transformer, now the central architecture in neural
forecasting. Its importance does not lie simply in introducing nonlinearity.
That shift had already begun. The transformer era is distinctive because it
makes attention and tokenization the organizing devices for neural forecasting.
We first introduce the basic attention mechanism and show how it replaces the
fixed lag coefficients of a VAR with weights that depend on the input window. We
then discuss the forecasting-specific pressures that shaped early time-series
transformers: long look-back windows, multi-step prediction, seasonality,
distribution shift, and nonstationarity. The final part of the section treats
tokenization as a central design choice, covering point-wise tokens,
patch-based models, and variable-token models. The probabilistic dimension
returns in Section~\ref{sec:diffusion}, where we discuss generative models.
\subsection{The basic transformer mechanism}
The transformer of \citet{vaswani2017attention} replaces recurrence with
self-attention. For a time-series reader, the central object is the attention
weight matrix. It plays a role analogous to lag weights in an autoregression,
but these weights are recomputed from the current input window.
Figure~\ref{fig:var2attn} previews this comparison. The summation structure of
the VAR is retained, while fixed coefficient matrices are replaced by
input-dependent weights.
\begin{figure}[t]
\centering
\includegraphics[width=0.98\textwidth]{fig_var_to_attention.png}
\caption{The central relaxation of this review. Left: the VAR forecasts with
fixed coefficient matrices $A_i$, estimated once and applied to every window.
This constancy is what defines impulse responses and Granger tests. Right:
attention forecasts with the same summation structure but with weights
$a_i(\cdot)$ that are functions of the input window and are recomputed at every
date. A seasonal window and a trending window can therefore induce different lag
profiles from the same fitted model. The gain is adaptivity and nonlinearity.
The loss is the fixed object on which the classical inferential apparatus rests.
Equation \eqref{eq:attn-var} states the correspondence precisely.}
\label{fig:var2attn}
\end{figure}
Let the input be a sequence of $L$ token vectors collected in the rows of
$X\in\mathbb{R}^{L\times d}$. Rows index token positions, while columns index embedding
coordinates. The meaning of a ``token'' in time series is deferred to
Subsection~\ref{sec:tokenize}. Attention forms queries, keys, and values through
learned projections,
\[
Q=XW_Q,\qquad K=XW_K,\qquad V=XW_V,
\]
with $W_Q,W_K\in\mathbb{R}^{d\times d_k}$ and $W_V\in\mathbb{R}^{d\times d_v}$. It then returns
\begin{equation}
\mathrm{Attention}(Q,K,V)
=
\mathrm{softmax}\!\left(\frac{QK^\top}{\sqrt{d_k}}\right)V.
\label{eq:attn}
\end{equation}
The softmax is applied to each row of its matrix argument. For a vector
$z=(z_1,\dots,z_L)^\top$, it is defined by
\begin{equation}
\mathrm{softmax}(z)_i
=
\frac{\exp(z_i)}{\sum_{j=1}^{L}\exp(z_j)},
\qquad i=1,\dots,L.
\label{eq:softmax}
\end{equation}
Each row of the resulting weight matrix is therefore nonnegative and sums to
one. Each output row of \eqref{eq:attn} is consequently a convex combination of
the rows of $V$.
Three properties are worth isolating. First, the row weights in \eqref{eq:attn}
are data dependent. Unlike fixed regression coefficients, they are
recomputed for each input from that input's own content. Second, the map is
global. Since rows of $X$ index token positions, any
two rows can interact directly through the $L\times L$ score matrix $QK^\top$.
Dependence between distant token positions therefore need not pass through
intermediate states. Third, the
normalization by $\sqrt{d_k}$ is substantive. When the entries of $Q$ and $K$
are approximately independent with unit variance, the inner products
$q^\top k$ have variance of order $d_k$. Without this scaling, the
softmax can put almost all weight on a single key, making the attention weights
hard to adjust during training.
A complete transformer block applies \eqref{eq:attn} in $M$ parallel heads,
typically with head dimension $d/M$. The head outputs are concatenated and mixed
by an output projection. A position-wise feed-forward network, residual
connections, and layer normalization complete the standard block. Because attention itself is permutation equivariant, ordering is supplied
separately through positional encodings. For forecasting, a causal mask sets the scores for future positions to
$-\infty$. After the softmax, those entries receive zero weight, so each position
can attend only to itself and earlier positions. Figure~\ref{fig:attnarch} summarizes the standard transformer block.
\begin{figure}[t]
\centering
\includegraphics[width=0.95\textwidth]{Attention.png}
\caption{Self-attention as a learned weighted estimator. The input $X$ is projected to queries, keys,
and values; the scaled scores $QK^{\top}/\sqrt{d_{k}}$ pass through a row-wise softmax to form the
weight matrix $A$, and the output $AV$ is a weighted average of the values, followed by a residual
connection, normalization, and a position-wise feed-forward map, over $M$ heads. Two econometric
readings are available: $A$ collects Nadaraya--Watson kernel weights under a learned kernel,
and equivalently the columns of $V$ act as learned common factors with the rows of $A$ as
data-dependent loadings.}
\label{fig:attnarch}
\end{figure}
\subsection{Adapting transformers to forecasting}
The connection to \eqref{eq:var} is exact in form. To make the comparison
explicit, let $E$ denote the embedding map that turns an observed vector into a
token vector. For a one-step forecast at date $t$, write
\[
x_s=E(\bm y_s),\qquad
X_t=[x_{t-L},x_{t-L+1},\dots,x_{t-1}]^\top\in\mathbb{R}^{L\times d}.
\]
Thus $X_t$ is the embedded look-back window. Raw observations enter the
attention computation only through this embedding. The construction of the
embedding, whether from single observations, patches, or quantized values, is
the subject of Subsection~\ref{sec:tokenize}.
For a one-step forecast, we use the attention output at the last observed
position, since the input window ends at $t-1$. The forecast can be written as
\begin{equation}\label{eq:attn-var}
\hat{\bm{y}}_{t}=\sum_{i=1}^{L}
\underbrace{a_i(\bm y_{t-L:t-1})}_{\text{data dependent}}
\,v(\bm y_{t-i}),
\qquad
v(\bm y_{t-i})=W_V^\top E(\bm y_{t-i}).
\end{equation}
Then the weight on lag $i$ is
\begin{equation}
a_i(\bm y_{t-L:t-1})
=
\frac{
\exp\!\left(q_t^\top k_{t-i}/\sqrt{d_k}\right)
}{
\sum_{j=1}^{L}
\exp\!\left(q_t^\top k_{t-j}/\sqrt{d_k}\right)
},
\qquad q_t=W_Q^\top E(\bm y_{t-1}),
\quad
k_{t-j}=W_K^\top E(\bm y_{t-j}).
\label{eq:aweight}
\end{equation}
Each $a_i(\bm y_{t-L:t-1})$ is nonnegative, and the weights sum to one across
$i=1,\dots,L$ by the softmax construction.
The query $q_t$ is computed from the most recent available token,
$E(\bm y_{t-1})$, and plays the role of the forecasting position at the end of
the look-back window. The keys $k_{t-j}$ are computed from the lagged tokens and
summarize the candidate past dates. The score $q_t^\top k_{t-j}$ therefore
measures how strongly the forecasting position attends to lag $j$. In this single-layer representation, $q_t$ depends only on
the last observed token $E(\bm y_{t-1})$. In a stacked transformer,
the representation at that position has already mixed information from earlier
tokens, so the query can depend on the wider history. Each attention weight also
depends on the full window through the softmax normalization across all lags.
Equation \eqref{eq:attn-var} is a weighted average of transformed lags. The
weights $a_i(\bm y_{t-L:t-1})$ play the same formal role as the coefficient
matrices $A_i$ in \eqref{eq:var}, but the analogy has two limits. In a VAR, the
coefficient matrices are fixed after estimation and are applied to every
window. In attention, the weights are functions of the current look-back window
and are recomputed at each date $t$. The attention representation is also
nonlinear. Nonlinearity enters through the value map $v(\cdot)$ and through the
softmax that produces the weights.
In the running macro example, the attention weights in \eqref{eq:attn-var} may give
different importance to past inflation or interest-rate episodes depending on the
current state of the economy. A recent high-inflation window and a tranquil window
can therefore induce different effective lag profiles from the same fitted model.
This adaptivity is useful for forecasting, but the resulting weights are not fixed
coefficients and should not be read as impulse responses.
The same operation also admits two useful econometric readings, as previewed in
Figure~\ref{fig:attnarch}: attention as a learned weighted estimator and
attention as adaptive pooling.
Firstly, for a single head, \eqref{eq:aweight} has the form of a
Nadaraya--Watson kernel smoother. The query plays the role of the evaluation
point, the keys play the role of design points, and the values are the
quantities being averaged. The similarity kernel is
\[
K(q,k)=\exp(q^\top k/\sqrt{d_k}).
\]
Attention differs from a classical smoother because similarity is not measured
in the raw input space. The matrices $W_Q$ and $W_K$ learn the representation in
which tokens are compared. The quantities being
averaged are also learned through $W_V$.
The second interpretation is adaptive pooling. Stacking \eqref{eq:attn} over
the sequence gives
\[
\mathrm{Attention}=AV,
\qquad
A=\mathrm{softmax}(QK^\top/\sqrt{d_k}), \qquad
V=XW_V .
\]
Here $A$ is a row-stochastic $L\times L$ weight matrix. Each row of $A$ gives
the weights used at one token position, and the rows of $V$ are the learned
value vectors being pooled. The operation therefore combines information across
tokens through weights that depend on the current input. This is related to the
pooling logic of factor models, but the analogy is limited. Dynamic factor
models \citep{stock2002forecasting,baing2002} summarize a high-dimensional
panel through a small number of common components with loadings fixed after
estimation. Attention instead forms input-dependent weighted averages of value
vectors. The resulting representation need not be low-dimensional. Multi-head attention repeats this adaptive pooling in parallel, and stacked
layers apply it repeatedly.
Because \eqref{eq:attn} constructs an $L\times L$ matrix of scores, its
computational cost scales quadratically with the context length. This becomes
costly in forecasting applications, where long look-back windows are often used
to capture seasonality, cycles, and other persistent dependence. An early wave
of time-series transformers therefore focused on reducing this cost.
\citet{zhou2021informer} retain only the most informative queries, reducing the
cost to order $L\log L$, and decode the forecast horizon in a single pass.
\citet{wu2021autoformer} embed a seasonal--trend decomposition and replace
inner-product attention with an auto-correlation operator.
\citet{zhou2022fedformer} move the computation to a frequency basis and obtain
linear complexity. These models reported improvements over the unmodified
transformer on standard long-horizon benchmarks, but the gains proved less
settled than they first appeared.
Distribution shift is a separate problem, and it is addressed by a separate
device. \citet{kim2022revin} propose reversible instance normalization (RevIN),
which removes the location and scale of each window before modelling and
restores them afterwards. It is a neural counterpart of the detrending,
demeaning, and rescaling that classical practice applies before estimation.
Later work models the shift rather than only removing it, feeding the discarded
information back into the attention computation \citep{liu2022nonstationary}.
The limits of these early transformer forecasters were made explicit by
\citet{zeng2023transformers}. They showed that DLinear, essentially a single
linear layer applied after seasonal--trend decomposition, matched or
outperformed these architectures on the same benchmarks. Their diagnosis was
that point-wise tokens carry limited information, autoregressive or iterated
decoding can accumulate errors, and some benchmarks reward smooth extrapolation,
where simple linear maps are already strong.
The lesson is that making attention cheaper or more elaborate is not enough if
the object being attended to is poorly chosen. The transformer approach was
therefore not abandoned. Instead, the literature reconsidered the unit over
which attention is computed. That unit is the token, and it is the subject of
the next subsection.
\subsection{Tokenizing time series}
\label{sec:tokenize}
Every architecture in this review operates on tokens, a representation inherited
from natural language processing. The borrowing is worth making explicit,
because it is both a source of these methods' power and a source of their
characteristic difficulties. In language modelling, a token is the elementary
unit consumed by the model. Since neural networks compute on vectors rather than
symbols, raw text must first be segmented into discrete units. Modern systems
typically use subword units, so a frequent word may be one token while a rare
word is assembled from several pieces.
Two steps should be kept distinct. Tokenization assigns a discrete index, which
carries no meaning by itself. Embedding maps that index to a learned vector, and
the usable representation resides in this vector space. A sentence therefore
becomes a sequence of embeddings, stacked as $X\in\mathbb{R}^{L\times d}$. This array is
the input to the attention operation in \eqref{eq:attn}
(Figure~\ref{fig:nlptoken}).
For a real-valued series there is no canonical token. The earliest
time-series transformers treated a single observation as a token. The DLinear
critique exposed the weakness of this choice: point-wise tokens carry little
local structure, and simple linear extrapolation can be hard to beat on smooth
long-horizon benchmarks. Later work therefore reconsidered the unit over which
attention is computed.
\citet{nie2023patchtst} propose PatchTST, which groups $P$ contiguous
observations into a single token. Each token then carries local shape rather
than a single value, and the sequence shortens from $L$ to $L/P$. This reduces
the cost of attention by a factor of roughly $P^2$. PatchTST also adopts
channel independence, using a shared univariate model across series. This avoids
fitting spurious cross-series correlations, but it also declines to model
cross-series dependence directly.
\begin{figure}[t]
\centering
\includegraphics[width=0.95\textwidth]{fig_nlp_tokenization.png}
\caption{Tokenization in language modelling. A sentence is segmented into subword tokens, here
``Forecasting'' splitting into \texttt{Fore}, \texttt{cast}, and \texttt{ing}; each token is assigned
an index in a fixed vocabulary, which carries no meaning, and each index is mapped by a learned
embedding to a vector, in which meaning resides. The embeddings, stacked as $X\in\mathbb{R}^{L\times d}$, are
the input to attention. Foundation models for time series replace the text tokens of this pipeline
with quantized values or real-valued patches (Figure~\ref{fig:tokenize}).}
\label{fig:nlptoken}
\end{figure}
\citet{liu2024itransformer} take the opposite view. In iTransformer, each
variable's whole history is treated as a token, so attention runs across the
cross-section rather than across time. Temporal structure is then handled mainly
inside the feed-forward component. The contrast with PatchTST is useful: one
design protects against spurious cross-series dependence by separating channels,
while the other makes cross-series dependence the object of attention.
A third approach follows the language pipeline more literally.
\citet{ansari2024chronos} quantize each scaled observation into a fixed
vocabulary of bins. The series becomes a sequence of discrete symbols, and an
unmodified language-model architecture can be trained on it. The price is a
resolution ceiling imposed by the binning. Quantized vocabularies and
real-valued patches are the two constructions that organize much of the
foundation-model literature, and they are contrasted in Figure~\ref{fig:tokenize}.
The choice is not merely computational. The token fixes the objects over which
attention is computed. For the macro panel, the model may compare individual
quarterly observations, short patches of local dynamics, or entire variable
histories. A patch token may encode a local inflation cycle. A variable token may
ask how the inflation series as a whole relates to the policy-rate series. The
token therefore fixes the economic object being compared. This is why
tokenization provides the natural bridge from time-series transformers to
foundation models.
\begin{figure}[t]
\centering
\includegraphics[width=0.95\textwidth]{fig_tokenization.png}
\caption{Two tokenization schemes for pretrained forecasters. Left: each scaled value is quantized into
one of a fixed set of bins, so that the series becomes a sequence of vocabulary symbols that a language
model continues (Chronos). Right: values are retained and contiguous timesteps are embedded as
real-valued patches (TimesFM, Moirai). Quantization imposes a resolution ceiling; real-valued patches
are compact and probabilistic by construction.}
\label{fig:tokenize}
\end{figure}
Read against the three challenges, the transformer era is uneven. High
dimensionality is handled mainly through implicit design choices. Channel
independence sidesteps the cross-section by not modelling it directly.
Cross-variable attention estimates dependence between series, but it provides no
test of which series matter. Tokenization also shapes the effective dimension of
the problem, since point-wise, patch-wise, and variable-wise tokens define
different objects over which attention is computed.
Nonstationarity is handled mostly through normalization. This can work well in
practice, but it is a device rather than a model. It says nothing about unit
roots or cointegration, so a model fitted to normalized data has no notion of a
long-run equilibrium. Later work models distribution shift rather than simply
removing it, feeding the discarded information back into the attention
computation \citep{liu2022nonstationary}. The change-point and locally
stationary methods of Section~\ref{sec:classical} have not yet been fully
integrated into these architectures.
Nonlinearity is where the main gain lies. Equation~\eqref{eq:attn} gives a
flexible operator whose weights adapt to the input, and it performs well when
data are plentiful. The same pattern persists in the models below: clear gains
against nonlinearity, ad hoc treatment of nonstationarity, and only implicit
treatment of dimensionality.
\section{Foundation models and zero-shot forecasting}
\label{sec:fm}
Until recently, forecasting was usually organized around a per-data-set
estimation workflow. For each application, the researcher selected a model,
estimated its parameters on the available series, and then re-estimated the
model as new data arrived. This workflow is difficult when the target series is
short, because there is little data from which to estimate the model. Foundation
models change the order of operations. The model is first pretrained on a large
and heterogeneous collection of series. It is then applied to a new series by
conditioning on that series' observed history, without re-estimating the model's
parameters. In this sense, forecasting becomes less an exercise in estimating a
new model for each data set and more an exercise in using a pretrained model
with the right context.
The analogy with language modelling is direct. A language model is pretrained on
many texts and then uses the prompt to generate a continuation. A time-series
foundation model is pretrained on many series and then uses the observed history
to generate a forecast. The time-series case, however, is harder in several
ways. Series differ in scale, sampling frequency, seasonal pattern, and degree
of nonstationarity. A model trained across many domains must therefore learn how
to compare and use histories that may live on very different numerical scales.
In the macro example, zero-shot forecasting would mean applying a pretrained
model to a new inflation-output-interest-rate panel without re-estimating its
parameters. The econometric concern is whether the model is truly learning from
the supplied history, or whether very similar macro series were already present
during pretraining. This is the contamination problem. It is why post-cutoff
evaluation is central: the evaluation series should come from a period that the
pretraining corpus could not have included.
The operative mechanism is in-context learning. The observed history is
given to the model as input, and the model adapts its forecast from that input
alone. Its weights are not updated. This distinction is useful. Zero-shot
forecasting refers to the setting in which the model forecasts a series that was
not used for training. In-context learning refers to how the model adapts in
that setting, by using the supplied history during a single forward pass.
Context length is therefore central. A longer context lets the model condition
on more of the target series. In multivariate settings, it can also include
related series supplied alongside the target. The forecast then depends not only
on what the model learned during pretraining, but also on how much relevant
history is available at prediction time.
The tokenizations introduced in Section~\ref{sec:tokenize} also organize the
pretrained models in Table~\ref{tab:fm}. Chronos \citep{ansari2024chronos}
represents the quantized-vocabulary approach. It reuses the language-model
apparatus and treats forecasting as classification over value bins. This makes
the connection to language modelling especially direct, but it also imposes a
resolution ceiling on continuous forecasts.
TimesFM \citep{das2024timesfm} and Moirai \citep{woo2024moirai} represent the
real-valued-patch approach. TimesFM is a decoder-only, long-context model.
Moirai is an any-variate, mixed-frequency model designed to handle different
numbers of variables and different sampling frequencies. Real-valued patches
have become a common design because they are compact, fast, and naturally suited
to continuous-valued prediction.
Other designs modify this template. Lag-Llama \citep{rasul2023lagllama} encodes
each date through lagged values. TimeGPT \citep{garza2023timegpt} is an early
closed-weight commercial system. Time-MoE \citep{shi2024timemoe} pursues
mixture-of-experts scaling. A recent development is the attachment of
generative output heads, based on flow matching or diffusion, to pretrained
backbones. This anticipates the convergence discussed in
Section~\ref{sec:diffusion}.
\begin{table}[t]
\centering
\caption{Time-series foundation models discussed in the text, classified by tokenization,
architecture, predictive output, and availability of weights. The two predominant tokenizations,
quantized vocabularies and real-valued patches, are developed in Section~\ref{sec:tokenize}.}
\label{tab:fm}
\small
\begin{tabular}{@{}lllll@{}}
\toprule
Model & Tokenization & Architecture & Output & Weights \\
\midrule
Chronos & quantized vocabulary & encoder--decoder & categorical & open \\
TimesFM & real-valued patches & decoder-only & quantiles & open \\
Moirai & real-valued patches & any-variate encoder & distributional & open \\
Lag-Llama & lagged covariates & decoder-only & distributional & open \\
TimeGPT & undisclosed & undisclosed & quantiles & closed \\
Time-MoE & point-wise real values & decoder-only MoE & point & open \\
\bottomrule
\end{tabular}
\end{table}
The empirical evidence is encouraging, but it should be read with care. On
contamination-aware benchmarks such as GIFT-Eval \citep{aksu2024gifteval},
leading foundation models often outperform tuned statistical methods and many
task-specific deep models in zero-shot mode (Table~\ref{tab:gifteval}).
Fine-tuning can improve performance, but the gains are often modest. This is
itself informative, since it suggests that zero-shot forecasting should be tried
before building a task-specific model.
The main caveat is contamination. A pretrained model's zero-shot score is
inflated if the evaluation series, or close relatives of it, appeared in the
pretraining corpus. This leakage is hard to verify because pretraining data are
often collected upstream and are not always fully disclosed. In econometric
terms, it is a leakage problem with no close analogue in the usual per-data-set
estimation setting. Benchmarks such as GIFT-Eval are therefore useful because
they make contamination a central evaluation issue rather than an afterthought.
Other limitations remain. Covariate support is uneven, irregular sampling is
still difficult, and very long horizons remain challenging. In adversarial or
strongly nonstationary domains, such as finance, the advantage of pretraining
may disappear entirely. The operational rule is therefore unchanged: evaluate on
genuinely held-out data, ideally post-cutoff and proprietary, and compare
against strong classical baselines.
\begin{table}[t]
\centering
\caption{Representative results from GIFT-Eval \citep{aksu2024gifteval}, computed from the public
result files in the GIFT-Eval repository as accessed on 14 July 2026. Values are normalized weighted
quantile losses (geometric mean across 97 benchmark tasks, relative to the seasonal-naive baseline),
so seasonal naive equals one and lower is better. The upper block reports task-specific supervised
baselines; the lower block reports zero-shot foundation models. TimesFM-2.5 denotes a later checkpoint
of the model introduced by \citet{das2024timesfm}, released without a separate publication. The
leaderboard is updated over time, so these values should be read as a dated snapshot.}
\label{tab:gifteval}
\small
\begin{tabular}{@{}lc@{}}
\toprule
Model & Normalized WQL ($\downarrow$) \\
\midrule
\multicolumn{2}{@{}l}{\emph{Task-specific baselines}}\\
\midrule
Seasonal naive & 1.000 \\
Auto-ARIMA \citep{hyndman2008automatic} & 0.912 \\
DeepAR \citep{salinas2020deepar} & 0.853 \\
DLinear \citep{zeng2023transformers} & 0.846 \\
N-BEATS \citep{oreshkin2020nbeats} & 0.816 \\
\midrule
\multicolumn{2}{@{}l}{\emph{Zero-shot foundation models}}\\
\midrule
Moirai-2.0 \citep{liu2025moirai2} & 0.516 \\
TimesFM-2.5 \citep{das2024timesfm} & 0.490 \\
Chronos-2 \citep{ansari2024chronos2} & 0.485 \\
\bottomrule
\end{tabular}
\end{table}
Against the three challenges, foundation models address high dimensionality
mainly through pretraining. Broad cross-corpus training lets the model pool
information across many domains, and any-variate designs \citep{woo2024moirai}
allow the number of input series to vary. This differs from regularized or
factor-based autoregressions, where parsimony is imposed through explicit and
interpretable restrictions.
Their treatment of nonstationarity is weaker. Scaling and normalization help
models operate across heterogeneous series, but they do not replace the
econometric apparatus of unit roots, cointegration, and error correction.
Nonlinearity is inherited from the transformer backbone. The central innovation
is therefore the replacement of in-sample estimation by cross-corpus
pretraining, whose statistical properties remain only partly understood.
Recent theory outside time series interprets in-context learning as implicit
Bayesian inference. In that view, the prompt updates a prior learned during
pretraining, and risk falls as context length grows
\citep{xie2022icl,muller2022pfn,wies2023learnability}. These results, however,
usually assume independent examples drawn from stationary task mixtures. A
forecasting context is different: it is a single dependent path, often
nonstationary. Extending the theory to that setting remains an open problem.
\section{Diffusion and generative forecasting}\label{sec:diffusion}
A useful forecast describes a range of possible futures, not just a single path. Decisions about inventory, capacity, and risk rely on quantiles and scenarios rather than a single point. Foundation models that output separate quantiles describe marginal uncertainty
at each horizon, but these quantiles do not by themselves characterize the joint
distribution across horizons and series. As a result, they cannot describe a coherent future scenario over the forecast horizon. Generative models address this limitation by sampling whole future paths (Figure~\ref{fig:why-gen}). Among generative methods, diffusion has emerged as a leading framework for time-series forecasting, building on earlier successes in image and audio synthesis.
In the running example, a diffusion forecaster would not merely produce separate
quantiles for future inflation, output growth, and the policy rate. It would sample
joint future paths, so that a high-inflation scenario comes with a corresponding
path for output and interest rates. The attraction is scenario coherence; the
econometric question is whether those scenarios respect the structural and long-run
relations thought to govern the macroeconomic system.
\begin{figure}[t]\centering
\includegraphics[width=0.95\textwidth]{fig_diff_scenarios.png}
\caption{Why generative forecasting. A point forecast is a single trajectory and conveys little information about uncertainty (left). A diffusion model samples many coherent future paths, which together represent the predictive distribution (right).}
\label{fig:why-gen}
\end{figure}
The contrast with Section~\ref{sec:classical} is direct.
In a Gaussian VARMA, the $h$-step predictive distribution is available in closed form,
\begin{equation}\label{eq:varma-pred}
y_{t+h}\mid y_{1:t}\ \sim\ \mathcal{N}\!\big(\hat y_{t+h},\,\Sigma_h\big),
\end{equation}
so future uncertainty is summarized by a forecast mean and a covariance matrix. Diffusion takes a different route:
\begin{equation}\label{eq:diff-pred}
y_{t+1:t+H}\mid y_{1:t}\ \sim\ p_\theta\!\big(\,\cdot \mid y_{1:t}\big),
\end{equation}
and learns to sample directly from the predictive distribution. It therefore does not impose a parametric constraint on its shape. The resulting distribution may be multimodal, heavy-tailed, or asymmetric in ways that a Gaussian VARMA \eqref{eq:varma-pred} cannot capture. In this sense, diffusion targets the same forecasting object as the classical model, but represents it more flexibly.
Density forecasting has been approached with several classes of generative models, including variational autoencoders, normalizing flows, and generative adversarial networks. Each learns to represent the predictive distribution by a different mechanism. We focus here on denoising diffusion, which has become a leading approach for time series forecasting and is the subject of the remainder of this section. Broader reviews of generative forecasting methods are available elsewhere \citep{meijer2024rise,yang2024survey}.
A key limitation should be noted at the start. Neither diffusion nor the foundation models of Section~\ref{sec:fm} enforce cointegration or long-run equilibrium by construction. Where such structure must hold, the vector error-correction model in \eqref{eq:vecm} remains the appropriate benchmark.
\subsection{Denoising diffusion: mechanism and foundations}
\label{sec:diff-mech}
To fix ideas, we begin with the vanilla diffusion mechanism and introduce forecasting-specific conditioning only in Section~\ref{sec:diff-models}.
Denoising diffusion combines a fixed forward process that gradually corrupts a target sample with a learned reverse process that reconstructs it (Figure~\ref{fig:diff-mech}). The basic idea goes back to \citet{sohldickstein2015}, who formulated generation as the finite-time reversal of a Markov chain that diffuses data into a tractable prior. It became practically useful with the denoising diffusion probabilistic models of \citet{ho2020denoising}. Throughout this subsection, let $x^{(0)}$ denote the target sample to be generated, with the superscript indexing diffusion step rather than calendar time.
In the forward process, each diffusion step shrinks the current state and adds Gaussian noise according to a prescribed variance schedule:
\begin{equation*}
x^{(n)}=\sqrt{\alpha_n}\,x^{(n-1)}+\sqrt{\beta_n}\,\epsilon_n,
\qquad \epsilon_n\sim\mathcal N(0,I).
\end{equation*}
Here $\{\beta_n\}_{n=1}^N$ is a prescribed sequence that controls the amount of noise added at each diffusion step, and $\alpha_n=1-\beta_n$. Writing $\bar{\alpha}_n=\prod_{s=1}^n \alpha_s$, iterating this recursion yields the closed form
\begin{equation}\label{eq:forward}
x^{(n)} = \sqrt{\bar\alpha_n}\,x^{(0)} + \sqrt{1-\bar\alpha_n}\,\epsilon,
\qquad \epsilon\sim\mathcal N(0,I).
\end{equation}
This representation allows training at any diffusion step without simulating the full chain. As $n$ grows, the distribution approaches the standard normal prior from which generation begins.
\begin{figure}[t]\centering
\includegraphics[width=0.95\textwidth]{fig_diff_forward.png}
\caption{The diffusion mechanism. A fixed forward process gradually adds noise to the series from left to right. A learned reverse process removes that noise from right to left. Generation starts from pure noise and ends in a sample. Conditioning on the observed history turns this generative process into forecasting.}\label{fig:diff-mech}
\end{figure}
For sufficiently small per-step variance, the reverse transition can be approximated by a Gaussian, with a mean predicted by the network. After the noise-prediction reparameterization of \citet{ho2020denoising}, the variational bound yields a timestep-weighted squared-error objective. In practice, one often optimizes the simplified objective
\begin{equation}\label{eq:diff-loss}
\mathbb{E}_{n,\,x^{(0)},\,\epsilon}
\big\|\,\epsilon-\epsilon_\theta\!\big(x^{(n)},n\big)\big\|^2,
\end{equation}
where the expectation is taken over a diffusion step $n$ drawn uniformly from $\{1,\dots,N\}$, a clean sample $x^{(0)}$, and Gaussian noise $\epsilon\sim\mathcal N(0,I)$. This objective asks the model to recover the noise injected in the forward process. Under this parameterization, the reverse update draws $x^{(n-1)}$ from a Gaussian with mean
\begin{equation}\label{eq:reverse}
\mu_\theta\!\big(x^{(n)},n\big)
= \frac{1}{\sqrt{\alpha_n}}\!\left(x^{(n)}
- \frac{\beta_n}{\sqrt{1-\bar\alpha_n}}\,\epsilon_\theta\!\big(x^{(n)},n\big)\right),
\end{equation}
so that generation proceeds by repeatedly applying the learned denoiser from level $N$ back to the data.
Up to a known scale factor, predicting the noise is equivalent to estimating the score, that is, the gradient of the log-density of the noised data. This is the perspective developed by \citet{song2019score}, who learn that gradient directly at a range of noise levels and generate samples by annealed Langevin dynamics. In that formulation, generation repeatedly follows the learned score field back toward high-density regions of the data distribution.
These two views, the variational one of \citet{ho2020denoising} and the score-based one of \citet{song2019score}, were unified by \citet{song2021scorebased} within a single stochastic differential equation (SDE) framework. In this formulation, the forward SDE gradually injects noise, while the reverse-time SDE, driven by the learned score, defines the generative process. The same framework also admits a deterministic probability-flow ordinary differential equation (ODE) with the same marginals. In continuous time, this ODE plays the same role as the deterministic denoising diffusion implicit model (DDIM) sampler of \citet{song2021ddim}.
A main practical drawback of classical diffusion sampling is that it requires tens to hundreds of denoising steps. The probability-flow ODE and DDIM already reduce this cost by replacing stochastic reverse sampling with more efficient deterministic trajectories. Subsequent developments push this idea further, including flow matching \citep{lipman2023flow}, which replaces the noising chain with a learned velocity field along a prescribed probability path, and, more recently, consistency models \citep{song2023consistency} and rectified flow \citep{liu2023rectified}. These methods can reduce generation to a handful of function evaluations, making deployment far cheaper.
The forward process also gives formal content to the analogy drawn in
Section~\ref{sec:classical}, though the comparison should be interpreted with
care. A single forward step follows the recursion
$x^{(n)} = \sqrt{\alpha_n}\,x^{(n-1)} + \sqrt{\beta_n}\,\epsilon_n$, whose
unrolling is exactly the closed form \eqref{eq:forward}. The noised state is
therefore a weighted accumulation of Gaussian terms $\{\epsilon_s\}_{s\le n}$,
with the same moving-average algebra as the Wold representation, in which the
VARMA outcome accumulates past innovations. The resemblance is formal
rather than substantive. VARMA innovations are disturbances of the data generating process, indexed by calendar time and carrying economic meaning. Diffusion terms, by contrast, are artificial noise injected by the modeller. They are indexed by noise level rather than calendar time and are used only to define the generative mechanism. What the two constructions share is a linear Gaussian recursion and the tractability it affords. In the VARMA, that recursion yields a closed-form predictive density. In diffusion, it yields the closed form \eqref{eq:forward}, which makes training feasible at arbitrary noise levels. The key departure is that diffusion learns to invert this accumulation, so the reverse map is nonlinear.
\subsection{Conditional diffusion for time-series forecasting}
\label{sec:diff-models}
To adapt vanilla diffusion to forecasting, let the target sample be the future segment
\[
x^{(0)} := \bm y_{t+1:t+H}.
\]
The forward noising process is unchanged. It perturbs only this future segment, while the observed history $\bm y_{1:t}$ enters separately as conditioning information.
In the unconditional setting of Section~\ref{sec:diff-mech}, the reverse chain learns transitions of the form
$p_\theta(x^{(n-1)} \mid x^{(n)}).$
For forecasting, these become conditional transitions
\[
p_\theta\!\big(x^{(n-1)} \mid x^{(n)}, \bm y_{1:t}\big).
\]
Equivalently, one may write the model in terms of a conditional denoiser
$\epsilon_\theta(x^{(n)}, n \mid \bm y_{1:t}).$
Running the reverse chain from the Gaussian prior therefore produces a sample from the conditional predictive distribution
\[
p_\theta\!\big(\bm y_{t+1:t+H} \mid \bm y_{1:t}\big),
\]
which is the forecasting object introduced in \eqref{eq:diff-pred}. In this sense, conditional diffusion forecasting preserves the vanilla generative mechanism and changes only the information supplied to the reverse map.
With this conditional formulation in place, the rest of the subsection organizes
diffusion forecasters into four groups.
Because time series have distinctive
structure, including temporal dependence, seasonality, trend, multiscale
patterns, and cross-variable interactions, plain diffusion algorithms are often
modified or re-engineered to fit these features.
The first group is organized by how the
forecast horizon is generated. Autoregressive models such as TimeGrad and
ScoreGrad produce one multivariate observation at a time, whereas whole-segment
models such as CSDI and SSSD generate the entire future segment jointly. The
second group focuses on how conditioning enters. TimeDiff strengthens the
conditioning signal during training, while TSDiff trains an unconditional model
and introduces the observed history only at sampling through self-guidance. The
third group adds time-series structure to the diffusion process itself.
Multiscale and decomposition models connect denoising to temporal aggregation,
seasonal--trend structure, or frequency-based components. The final group
changes the role of diffusion more substantially, using it either to support
denoising and latent representation learning in short noisy series or to model
dispersion around a transformer point forecast. Table~\ref{tab:diff-models}
summarizes the models discussed below.
\begin{figure}[t]\centering
\includegraphics[width=0.95\textwidth]{fig_diff_ar_vs_block.png}
\caption{Autoregressive versus whole-segment generation. TimeGrad generates
one date at a time, so the predictive band widens and error compounds along the
horizon (left); CSDI and SSSD generate the whole horizon jointly, avoiding recursive feedback of previously sampled values (right).}\label{fig:ar-vs-joint}
\end{figure}
\paragraph{Autoregressive versus whole-segment generation.}
The future horizon can be produced in two ways. A model may generate one date at a time, feeding each
sampled value back as an input for the next. Or it may generate the whole segment
$\bm y_{t+1:t+H}$ at once, in a single pass of the reverse chain.
The distinction is not new to diffusion. It is the generative version of an older choice in
forecasting, between iterated and direct multi-step prediction
\citep{marcellino2006comparison}. Iteration exploits the recursive structure of the process. It
is efficient when the one-step model is correct, but error compounds along the horizon, because
each step conditions on values the model has itself generated. Whole-segment generation avoids
that accumulation. It is more robust to a misspecified one-step model, but it must capture
dependence across the entire horizon in a single map.
\citet{rasul2021timegrad} propose TimeGrad, the first diffusion model for multivariate
forecasting. It summarizes the observed history in the hidden state of a recurrent network. That
state is passed to the denoiser as conditioning information, and the model generates one future
date at a time.
The samples are coherent across series. But the design inherits both costs of iteration. A full
denoising loop is run at every forecast date and error accumulates along the
horizon (Figure~\ref{fig:ar-vs-joint}, left). \citet{yan2021scoregrad} introduce ScoreGrad, the continuous-time counterpart
of this autoregressive design. The Recurrent Neural Network (RNN) still summarizes the history and the
forecast is still generated one date at a time, but the conditional distribution
at each date is modeled by a score-based stochastic differential equation. This replaces the fixed denoising loop with numerical
solvers for the reverse-time SDE, while leaving the autoregressive horizon
factorization unchanged. ScoreGrad therefore inherits the same possibility of
horizon-wise error accumulation.
A different line of work begins with CSDI \citep{tashiro2021csdi}, a conditional score-based diffusion model originally developed for imputation. The key observation is that forecasting can be viewed as imputation in which the missing segment lies at the end of the window. CSDI therefore masks the future and generates the entire horizon jointly, conditioned on the observed past. This removes autoregressive drift and allows the same model to handle ordinary missing-data imputation (Figure~\ref{fig:ar-vs-joint}, right). \citet{alcaraz2023sssd} propose SSSD, which keeps this whole-segment, mask-based design and
changes only the backbone. Attention is replaced by the structured state-space (S4) layers of
\citet{gu2022s4}. These capture long-range dependence at near-linear cost and scale more
naturally to long sequences.
\paragraph{Conditioning refinements.}
A second set of models changes not the horizon factorization, but the way conditioning is
supplied.
\citet{shen2023timediff} propose TimeDiff, a non-autoregressive forecaster that strengthens the conditioning signal used by the denoiser. Two devices do the work. Future mixup blends ground-truth future values into the conditioning signal during training, so that the denoiser sees future-compatible structure beyond what the history alone provides. When forecasting, however, these ground-truth values are unavailable and the condition is replaced by the history-based representation. Autoregressive initialization then starts the reverse chain from a cheap linear forecast rather than from pure noise. The model therefore refines an informed trajectory instead of building the future segment from scratch. Both devices are forecasting-specific design choices that improve empirical performance.
\citet{kollovieh2023tsdiff} take the opposite route and introduce TSDiff, an
unconditional diffusion model for time-series windows. Unlike conditional
forecasters, TSDiff is trained with no forecasting condition supplied during
training. The observed history enters only at sampling time through
self-guidance. The motivation is the decomposition
\begin{equation}
\nabla_x \log p\big(x \mid y_{1:t}\big)
= \nabla_x \log p\big(x\big) + \nabla_x \log p\big(y_{1:t} \mid x\big),
\label{eq:guidance}
\end{equation}
which shows that an unconditional score can be converted into a conditional one
by adding a guidance correction. TSDiff approximates this correction at each
denoising step by first using the denoiser to estimate the clean window implied
by the current noisy sample. It then compares the history part of that estimate
with the observed history and nudges the reverse update toward better agreement.
The resulting sampler produces realistic windows whose past matches the observed
history, and the future part of each window is read as the forecast.
Guidance buys generality: one trained model can be steered toward forecasting, forecast
refinement, or synthetic generation without retraining. It costs exactness, since the correction
in \eqref{eq:guidance} is approximated rather than computed, so the guided sampler is not
guaranteed to draw from the true conditional predictive law. Conditioning supplied during
training targets that law more directly, but ties the model to the task it was trained for.
\paragraph{Multiscale and decomposition.}
A further group builds temporal structure into diffusion by linking the noise
ladder to coarse-to-fine representations of the series
(Figure~\ref{fig:multiscale}). The intuition is that high noise levels preserve
only broad low-frequency structure, while lower noise levels recover local
fluctuations and high-frequency detail. This makes the reverse process resemble
a coarse-to-fine reconstruction of the time series.
\citet{fan2024mgtsd} introduce MG-TSD, which derives coarse-grained versions of the series by temporal aggregation and uses them as intermediate guidance targets during diffusion training. The motivation is that both diffusion noising and temporal aggregation remove fine-scale variation, leaving smoother low-frequency structure. The reverse process is therefore encouraged to recover broad temporal patterns before refining high-frequency detail.
\citet{shen2024mrdiff} propose mr-Diff, which makes the hierarchy explicit by
generating a seasonal--trend decomposition from coarse to fine, with each finer
stage conditioned on the coarser one. \citet{yuan2024diffts} introduce
Diffusion-TS and place the decomposition inside the denoiser itself: each
denoising estimate is split into interpretable trend and Fourier-based seasonal
components. Together, these methods make diffusion forecasting less purely generic by tying
denoising to temporal scale and decomposition.
\begin{figure}[t]\centering
\includegraphics[width=0.9\textwidth]{fig_diff_coarse_to_fine.png}
\caption{Multiscale structure in diffusion. MG-TSD uses coarse-grained series
as guidance targets at intermediate diffusion levels, whereas mr-Diff explicitly
generates the forecast from coarse trends to finer details.}\label{fig:multiscale}
\end{figure}
\paragraph{Distinct designs.}
The preceding models differ mainly in how they condition or guide a diffusion
sampler whose target is the future segment. The following two models change the
role of diffusion itself.
\citet{li2022d3vae} propose D$^3$VAE for short and noisy time series, where whole-segment diffusion models can be difficult to fit with limited data. In this model, diffusion is
used less as a standalone generator and more as a denoising and augmentation
device. A coupled forward process perturbs the history and the future together,
creating noisy versions of a small training sample. A bidirectional variational
autoencoder then learns to invert this process. The objective combines a
denoising term, which removes both injected and intrinsic noise, with a
disentanglement penalty that separates latent factors.
\citet{li2024tmdm} propose TMDM, which combines a transformer point forecaster
with a diffusion model for residual uncertainty. The transformer first produces a
conditional mean forecast. Rather than diffusing toward a standard normal prior,
the forward process diffuses toward a Gaussian centered at that forecast,
$\mathcal{N}(\hat y,I)$, and the same transformer output also modulates the
reverse chain. Diffusion therefore models plausible deviations around a strong
conditional mean, rather than generating the future segment from scratch. This
design connects the point-forecasting models of
Sections~\ref{sec:transformers}--\ref{sec:fm} with the generative
forecasting methods of this section.
\begin{table}[t]\centering
\caption{Diffusion forecasters, classified by generation scheme, core design,
and conditioning mechanism.}
\label{tab:diff-models}
\begin{tabular}{llll}
\toprule
Model & Generation & Core design & Conditioning \\
\midrule
TimeGrad & autoregressive & RNN & recurrent state \\
ScoreGrad & autoregressive & RNN + SDE & recurrent state \\
CSDI & whole-segment & attention & masked observations \\
SSSD & whole-segment & state-space (S4) & masked observations \\
TimeDiff & whole-segment & convolutional & future mixup + AR init \\
TSDiff & unconditional & state-space (S4) & sampling-time self-guidance \\
MG-TSD & autoregressive & RNN (multi-granularity) & coarse-grained guidance \\
mr-Diff & whole-segment & seasonal--trend & coarse-to-fine stages \\
Diffusion-TS & whole-segment & decomposition attn. & guidance \\
D$^3$VAE & whole-segment & diffusion + BVAE & latent \\
TMDM & whole-segment & transformer + diff. & conditional mean \\
\bottomrule
\end{tabular}
\end{table}
The diffusion literature therefore answers the three challenges unevenly. It is strongest where the problem is joint uncertainty, less settled where the
issue is nonstationary structure, and broadly flexible for nonlinear and
non-Gaussian predictive shapes.
Diffusion is most distinctively valuable for
high dimensionality. By learning a joint predictive law over horizons and series,
diffusion produces coherent scenario paths rather than separate marginal
quantiles. This is precisely what quantile heads alone cannot provide.
For nonstationarity, the picture is less settled. For nonstationarity, the picture is less settled. Recent work has begun to
incorporate time-series characteristics into diffusion design. For example,
\citet{lee2024ant} select the noise schedule using autocorrelation-based measures
of nonstationarity. This remains an adaptive scheduling device rather than a
formal treatment of unit roots or cointegration. The vector
error-correction form \eqref{eq:vecm} therefore remains the appropriate
guardrail when a long-run equilibrium must hold.
Nonlinearity and non-Gaussianity are handled more naturally. Diffusion imposes no
parametric form on the predictive distribution, so it can represent multimodal,
heavy-tailed, and asymmetric forecast laws beyond the Gaussian VARMA. The broader
trajectory is to place such generative heads on top of pretrained backbones,
combining transformer representations, pretraining-based generalization, and
diffusion-based scenario uncertainty.
\section{Synthesis and outlook}
\label{sec:synthesis}
Table~\ref{tab:synth} summarizes the argument. Read across the challenge columns, no family yet
meets high dimensionality, nonstationarity, and nonlinearity with the completeness applied work
requires, and the pattern is uneven. Progress on nonlinearity and non-Gaussian predictive shape
has been rapid. High dimensionality is handled implicitly, through pooling, channel choices,
pretraining, or joint sampling, rather than through the interpretable restrictions of factor and
shrinkage methods. Nonstationarity is met with normalization or adaptation, not with the unit-root
and cointegration apparatus that economic series often require. The structural and inferential
machinery of the classical column, orthogonalized shocks, testable restrictions, and
equilibrium-restoring terms, has no clear counterpart further right.
\begin{table}[t]
\centering
\small
\begin{tabular}{@{}p{2.4cm}p{3.0cm}p{3.0cm}p{3.0cm}p{2.6cm}@{}}
\toprule
\textbf{Family} & \textbf{High dimensionality} & \textbf{Nonstationarity} & \textbf{Nonlinearity}
& \textbf{Inference} \\
\midrule
VAR / VECM & factor, shrinkage, sparse priors (explicit) & unit roots, cointegration (explicit) &
linear, fixed coefficients & full distribution theory \\
\addlinespace
Transformers & channel independence / cross-variable attention & instance normalization (RevIN) &
data-dependent attention & none standard \\
\addlinespace
Foundation models & cross-corpus pretraining; any-variate & scaling across regimes (implicit) &
inherited from backbone & none standard \\
\addlinespace
Diffusion & joint predictive law; no explicit dimension reduction & limited explicit treatment & non-Gaussian,
multi-modal & sampling-based only \\
\bottomrule
\end{tabular}
\caption{How each method family treats the three challenges, and what inferential apparatus it
provides. ``Explicit'' denotes an interpretable, testable restriction; ``implicit'' a data-driven
mechanism without such guarantees.}
\label{tab:synth}
\end{table}
These advances make the modern architectures flexible representations of the conditional law of the
future given the past. Flexibility is not by itself an econometric standard. For a method to inform
policy, risk management, or scientific interpretation, its use must also be governed by credible
evaluation, calibrated uncertainty, explicit identification, and a clear account of the
nonstationarity it is meant to handle.
Evaluation should mimic the information set available at the forecast origin. This rules out random
train--test splits for dependent data, in favour of rolling-origin or expanding-window designs. For
pretrained models a further concern is contamination. A zero-shot score is uninformative if the
target series, or close variants of it, entered the pretraining corpus, so the strongest evidence
comes from post-cutoff or proprietary data, or from benchmarks that separate the pretraining and
evaluation corpora. Modern models should be compared against strong classical baselines, including
seasonal naive forecasts, ARIMA, VAR or Bayesian VAR, and a VECM where cointegration matters.
A forecast is often useful only with a statement of uncertainty. Classical Gaussian VARs supply
this through an explicit sampling theory. Modern neural forecasters usually do not. Quantile heads,
diffusion samplers, and ensembles produce empirical predictive distributions, but these should not
be read as calibrated without a check. Calibration can be assessed through interval coverage,
probability integral transform diagnostics, and horizon-by-horizon reliability. Where
distributional assumptions are weak, conformal prediction offers model-agnostic finite-sample
coverage under exchangeability or approximate-stability conditions. The operative question is not
whether a model gives sharp intervals, but whether those intervals have their advertised coverage
under the instability faced in the application.
Forecast accuracy does not identify structural effects. The reduced-form VAR forecasts well, but its innovations are correlated, so an impulse
response becomes economically meaningful only after identifying restrictions, such as a recursive
ordering, sign restrictions, external instruments, or long-run restrictions. The flexible methods inherit this problem and add nothing
to its solution. Attention weights, learned kernels, and diffusion trajectories are functions of the
fitted predictive system, not structural shocks. A data-dependent weight is not an identified parameter, and its sampling distribution is
unknown. Structural use therefore requires two steps that accuracy cannot supply. First, define the
estimand, for example a nonlinear analogue of an impulse response or a counterfactual policy path.
Second, supply assumptions that connect it to the data, by imposing econometric restrictions on the
flexible model or by combining it with external identification. In the macro panel, if the question
is the effect of a monetary tightening, the model must be linked to an identified policy shock. A
diffusion path in which the interest rate rises does not by itself identify one. Absent such
structure, modern forecasters are predictive devices rather than structural models.
Nonstationarity is not a single problem. Series may carry stochastic or deterministic trends,
structural breaks, evolving volatility, or gradual parameter drift, and classical econometrics
distinguishes these cases because each calls for a different response. Modern systems tend to
address them together, through scaling, normalization, long context windows, or pretraining across
regimes. These devices improve robustness but do not replace an explicit account of the
instability, and differencing or normalizing away low-frequency movement can destroy long-run
equilibrium information when series are cointegrated.
Taken together, the move from vector autoregressions to transformers, foundation models, and
diffusion is best understood as an expansion of the forecasting class, not a replacement of
econometric structure. Self-attention keeps the VAR's weighted sum of past values, but lets the
weights depend on the input and adds nonlinearity. Foundation models replace per-data-set estimation
with pretraining across many series and zero-shot adaptation to new ones. Diffusion replaces a single
parametric predictive density with samples of coherent future paths. Each step adds flexibility,
scale, and richer predictive distributions. None of them, by itself, supplies the explicit
restrictions, identification, and inferential discipline that structural interpretation and policy
analysis require. In particular, the burden of dimensionality is partly shifted from estimating many
parameters on one data set to learning shared representations from large and heterogeneous corpora.
It is striking that modern machine-learning methods address the dimensionality problem in such a
natural way: rather than specifying low-dimensional restrictions by hand, they learn reusable
representations through pooling, attention, and pretraining across many related series.
The task ahead is to combine the predictive strength of modern architectures with the structural
and inferential tools of econometrics. Developing reliable inference and interpretability for these
models is a promising direction, because it would clarify when their predictive distributions are
calibrated, when learned restrictions can be interpreted econometrically, and how far their
black-box representations can support policy or scientific claims rather than forecasts alone.
Determining when implicit restrictions, learned through pretraining and pooling, outperform explicit
econometric structure under dependence, heterogeneity, and instability remains, in our view, the
central open problem this literature poses.
\bibliographystyle{plainnat}
\bibliography{references}