# Forecasting Electricity Market Signals via Generative AI††thanks: This work
was supported in part by the National Science Foundation under Award EECS
2218110.

Xinyi Wang , Qing Zhao , Lang Tong

###### Abstract

This paper presents a generative artificial intelligence approach to
probabilistic forecasting of electricity market signals, such as real-time
locational marginal prices and area control error signals. Inspired by the
Wiener-Kallianpur innovation representation of nonparametric time series, we
propose a weak innovation autoencoder architecture and a novel deep learning
algorithm that extracts the canonical independent and identically distributed
innovation sequence of the time series, from which future time series samples
are generated. The validity of the proposed approach is established by proving
that, under ideal training conditions, the generated samples have the same
conditional probability distribution as that of the ground truth. Three
applications involving highly dynamic and volatile time series in real-time
market operations are considered: (i) locational marginal price forecasting
for self-scheduled resources such as battery storage participants, (ii)
interregional price spread forecasting for virtual bidders in interchange
markets, and (iii) area control error forecasting for frequency regulations.
Numerical studies based on market data from multiple independent system
operators demonstrate the superior performance of the proposed generative
forecaster over leading classical and modern machine learning techniques under
both probabilistic and point forecasting metrics.

###### Index Terms:

Probabilistic forecasting, electricity price forecasting, representation
learning, scenario generation, generative artificial intelligence, bidding
strategies for energy and ancillary service markets, and risk-sensitive market
operations.

##  I Introduction

A defining feature of generative artificial intelligence (AI) is its ability
to produce artificial samples that resemble reality. In particular, generative
AI learns the underlying structure of a phenomenon from examples, from which
an arbitrarily large number of artificial samples having the same properties
exhibited in the examples are produced. The success of generative AI is well
documented [1].

The classical subject of forecasting has a natural connection with generative
AI through probabilistic forecasting. Given past observations, probabilistic
forecasting aims to obtain the conditional probability distribution of future
time series given the past observations. Once such a distribution is obtained,
Monte-Carlo samples of the future time series can be generated. The pre-deep-
learning era approaches to probabilistic forecasting are mostly model-based.
By imposing a specific parametric structure on the conditional distribution,
e.g., Gaussian distribution with mean and variance as parameters,
probabilistic forecasting becomes estimating parameters of the probability
distribution from past observations using point estimation or Bayesian
techniques. However, such parametric approaches cannot adequately capture
temporal dependencies for highly complex time series such as natural language
time series. To this end, generative AI techniques involving sophisticated
neural networks and novel machine learning methods have achieved remarkable
performance in real-world applications unmatched by conventional methods.

This paper develops a generative probabilistic forecasting (GPF) solution
using a generative AI approach to forecasting electricity market signals such
as locational marginal prices (LMPs), area control errors (ACEs), and more. In
particular, we develop a multivariate nonparametric forecasting technique that
directly generates future time series samples based on current and past
observations, bypassing computational and sample complexity challenges of
estimating the conditional probability distribution. Currently, such
nonparametric GPF techniques with provable performance are lacking.

###  I-A Literature Review

The literature on parametric and nonparametric probabilistic forecasting is
vast. See earlier reviews on the subject in [2, 3, 4]. Our review will focus
on some of the classic and modern machine learning-based techniques over the
last decade. To this end, we pay particular attention to methods suitable for
forecasting market operational quantities such as LMP and power flows.
Techniques for load and renewable generation forecasting are not included.

Mainstream probabilistic forecasting techniques are from parametric or
nonparametric families. Parametric forecasting predicts a parameterized
conditional probability distribution of future time series variables, reducing
an infinite-dimensional inference problem to one of finite-dimensional
estimation. Popular approaches include autoregressive and moving average
models [5, 6, 7, 8, 9], Gaussian models [10], Student’s $t$-distribution [11],
and others [12]. A classic benchmark with strong performance is SNARX [4],
which we include in our comparison. While restricting probability
distributions to be parametric results in tractable computations, it comes
with the loss of accuracy from model mismatch.

Nonparametric forecasting has a long history. See [2] for a review up to 1997.
These techniques estimate the underlying probability distribution or its
derived properties (e.g., quantile) without assuming a parametric form.
Classical techniques [13] face formidable sample and computational
complexities, further magnified when the underlying time series has arbitrary
temporal dependencies. Quantile regression is one of the most popular
techniques for electricity prices. By estimating multiple quantiles, one
obtains a probability mass function as an approximation. Over the last decade,
modern machine learning techniques are widely adopted to compute the
conditional quantiles. Examples include Extreme Learning Machine (ELM) [14],
Recurrent Neural Network (RNN) [15], and attention mechanism [16].

Most machine learning-based GPF methods adopt an autoencoder architecture from
which samples of future time series are generated. Examples include
Variational AutoEncoder (VAE) based methods [17, 18, 19, 20], normalizing
flows [21], and denoising diffusion probabilistic models [22]. These
techniques employ likelihood-based learning, implicitly assuming independent
data segments to compute a likelihood upper bound. For time series with
temporal dependencies, the independent data assumption is tenuous.

GAN-based GPF methods [23, 24] have also been proposed. They involve a
generator to map past observations and exogenous IID samples to future time
series samples. The autoencoder is trained by matching generated samples with
historical data. Theoretical guarantees of such techniques to produce the
correct samples are missing.

Since GPF can produce almost any type of point forecasts, we review a few
point forecasters with established reputations in some applications. Large
language model (LLM)-based techniques have been developed for point
forecasting of time series [25, 26]. We include transformer-based point
techniques as benchmarks for comparison in Sec. IV, given their demonstrated
superior performance for natural languages and other time series. Informer
[27] and Pyraformer [28] are point forecasting techniques that adopted the
popular transformer model to capture the temporal dependencies. For their
award-winning contributions, we include [27] and [28] in our comparative study
in this paper.

###  I-B Summary of Contributions

The main contribution of this work is a weak innovation autoencoder (WIAE)
based GPF, abbreviated as WIAE-GPF. Inspired by the classical Wiener-
Kallianpur (strong) innovation representation of time series [29], WIAE-GPF
distinguishes itself from existing autoencoder-based forecasters by
constraining the latent process to be an independent and identical uniformly
distributed (IID-uniform) innovation sequence, which is shown to be Bayesian
sufficient for forecasting future time series. From this structure, we
establish the validity of WIAE-GPF in Theorem 1, showing that samples
generated by WIAE-GPF have the same conditional distribution as the ground
truth under ideal training conditions. Such a theoretical guarantee is rare
for nonparametric GPF methods.

Because practical implementations of WIAE are necessarily finite-dimensional,
we establish a structural convergence property that the conditional
distribution of the WIAE-GPF output converges to that of the conditional
distribution of the time series. These theoretical contributions make the
proposed WIAE forecasting the first approach that achieves asymptotic
consistency of GPF under idealized learning conditions***To establish
asymptotic optimality, we assume that the Wasserstein-GAN learning achieves
global optimality with infinite training data set.. Sec. III for details.

There have been some but limited applications of generative AI techniques in
power system operations despite their accelerated advances in representing and
learning time series models. Missing in particular are the validation and
comparative studies using real-world market data. We fill this gap by
comparing the WIAE-GPF with leading traditional and machine-learning
algorithms for three applications: real-time LMP forecasting for energy
markets, interregional LMP spread forecasting for interchange markets, and ACE
forecasting for regulation markets. Such comparisons are essential because
these real-time market signals have characteristics not present in media
signals such as video and natural language time series, where machine learning
techniques have demonstrated success. Both LMP and ACE are highly dynamic time
series with frequent spikes. Our comparison study offers a compelling case for
WIAE-GPF across multiple performance measures for point and probabilistic
forecasters.

The idea of WIAE-GPF was first presented in a preliminary version of this work
[30], from which the current paper makes substantial new contributions†††Based
on Turnitin comparison, this paper exhibits less than 15% percent overall
similarity and less than 4% similarity to the preliminary version.. In
particular, the Bayesian sufficiency theorem in Sec. III-A is significantly
stronger than that in [30]. Also new are a theorem (Theorem 1 in Sec. III-A)
that establishes formally the validity of WIAE-GPF and a theorem on the
structural convergence (Theorem 2 in Sec. III-B). We considered three specific
real-time market applications, two of the three were not considered in [30].
The numerical results for all three applications as well as the analysis and
discussions are all new.

###  I-C Organization and Notations

This paper is organized as follows. Sec. II presents a nonparametric
multivariate model and innovation representations. Sec. III develops WIAE-GPF,
the proposed GPF algorithm. Sec. IV presents the application of WIAE-GPF in
three market operations and the comparison studies of major forecasting
benchmarks.

The notations used in this paper are standard. Random variables are in capital
letters and their realizations in lowercases. We use $({\bm{X}}_{t})$ to
denote a multivariate random time series, where column vector
${\bm{X}}_{t}=(X_{1t},\cdots,X_{dt})$ is the time series at time $t$, and
$(X_{it})$ the $i$th sub-time series of $({\bm{X}}_{t})$. In this paper,
${\bm{X}}_{t_{1}:t_{2}}$ denotes the segment of $({\bm{X}}_{t})$ from $t_{1}$
to $t_{2}$, i.e.,
${\bm{X}}_{t_{1}:t_{2}}=({\bm{X}}_{t_{1}},\cdots,{\bm{X}}_{t_{2}})$. For two
random vectors ${\bm{X}}$ and ${\bf Y}$,
${\bm{X}}\stackrel{{\scriptstyle\mbox{\tiny a.s.}}}{{=}}{\bf Y}$ means the two
random variables equal almost surely, and
${\bm{X}}\stackrel{{\scriptstyle\mbox{\tiny d}}}{{=}}{\bf Y}$ means the two
equal in distribution. An IID random sequence with marginal distribution
cumulative distribution $F$ is denoted by
${\bm{X}}_{t}\stackrel{{\scriptstyle\mbox{\tiny\sf IID}}}{{\sim}}F$. Table I
shows the major designated symbols used in the paper.

TABLE I: Abbreviations and mathematical notations used in this paper.

GPF | Generative Probabilitic Forecasting.
---|---
WIAE | Weak Innovation AutoEncoder.
ACE | Area Control Error.
LLM | Large Language Model.
NMSE | Normalized Mean Square Error.
NMAE | Normalized Mean Absolut Error.
MASE | Mean Absolute Scaled Error.
sMAPE | Symmetric Mean Absolute Percentage Error.
CRPS | Continuous Ranked Probability Score.
CP | Coverage Probability.
CPE | Coverage Probability Error.
NCW | Normalized Coverage Width.
$({\bm{X}}_{t})$ | The random process of predictive interests.
$({\bf V}_{t})$ | The innovation sequence.
$({\bf U}_{t})$ | An IID sequence of uniform distribution.
$(\hat{{\bm{X}}}_{t})$ | The rescontruction sequence output by WIAE decoder.
$\left(\hat{{\bf V}}_{t}^{(m)}\right)$ | The weak innovation sequence estimated by a $m$-dimensional WIAE.
$\left(\hat{{\bm{X}}}_{t}^{(m)}\right)$ | The reconstruction sequence estimated by a $m$-dimensional WIAE.
$({\bm{x}}_{t})$ | A sequence of real numbers indicating the past realizations of $({\bm{X}}_{t})$.
$(\hbox{\boldmath$\nu$\unboldmath}_{t})$ | A sequence of real numbers indicating the past realizations of $({\bf V}_{t})$.
$G$ | WIAE encoder function.
$H$ | WIAE decoder function.
$G_{\theta}$ | A neural network approximation of $G$ parameterized by $\theta$.
$H_{\eta}$ | A neural network approximation of $H$ parameterized by $\eta$.
$D_{\gamma}$ | Innovation discriminator that measures the distance between $({\bf V}_{t})$ and $({\bf U}_{t})$.
$D_{\omega}$ | Reconstruction discriminator that measures the distance between $({\bm{X}}_{0:t+T})$ and $({\bm{X}}_{0:t},\hat{{\bm{X}}}_{t+1:t+T})$.
${\cal U}[0,1]^{d}$ | The continuous $d$-dimensional uniform distribution on $[0,1]$.

##  II Innovation Representation Learning

###  II-A Strong and Weak Innovation Representations

In 1958, Wiener and Kallianpur proposed an innovation representation of scalar
time series [29]. In the parlance of modern machine learning, an innovation
representation is a causal autoencoder shown in Fig. 1 with the latent process
$({\bf V}_{t})$ being an IID-uniform innovation sequence. In particular, ${\bf
V}_{t}$ represents the new information (innovation) contained in
${\bm{X}}_{t}$ independent of the past
${\bm{X}}_{0:t-1}=({\bm{X}}_{t-1},{\bm{X}}_{t-2},\cdots)$. Mathematically, the
innovation representation of the time series is defined by causal mappings
$(G,H)$ and $({\bf V}_{t})$:

$\displaystyle{\bf V}_{t}$
$\displaystyle=G({\bm{X}}_{t},{\bm{X}}_{t-1},\cdots),$ (.1)
$\displaystyle({\bf V}_{t})\stackrel{{\scriptstyle\mbox{\sf\tiny
IID}}}{{\sim}}{\cal U}[0,1]^{d},$ (.2) $\displaystyle\hat{{\bm{X}}}_{t}$
$\displaystyle=H({\bf V}_{t},{\bf V}_{t-1},\cdots),$ (.3)

Figure 1: An autoencoder interpretation of innovation representations.

The Wiener-Kallianpur’s innovation autoencoder requires further that the
decoder output $(\hat{{\bm{X}}}_{t})$ reconstructs the input $({\bm{X}}_{t})$
(almost surely), i.e., $(\hat{{\bm{X}}_{t}})\stackrel{{\scriptstyle\mbox{\tiny
a.s.}}}{{=}}({\bm{X}}_{t})$, which makes Wiener-Kallianpur’s autoencoder a
strong innovation Autoencoder. The perfect causal reconstruction implies that
the innovation sequence $({\bf V}_{t})$ is a sufficient statistic for all
decision-making based on $({\bm{X}}_{t})$. Therefore, using the IID-uniform
$({\bf V}_{t})$ for decision-making incurs no performance loss.

However, Rosenblatt showed that the Wiener-Kallianpur (strong) innovation
representation does not exist for broad classes of random processes, including
some of the widely used finite-state Markov chains [31]. Rosenblatt suggested
a weak innovation representation, requiring that the autoencoder output
$(\hat{{\bm{X}}}_{t})$ matches its input $({\bm{X}}_{t})$ only in
distribution:

$\displaystyle({\bm{X}}_{0:t},\hat{{\bm{X}}}_{t+1:t+T})\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{=}}({\bm{X}}_{0:t+T}),\forall t.$ (2)

Herein, we call the autoencoder $(G,H)$ for the weak innovation representation
the Weak Innovation Auto Encoder (WIAE).

###  II-B Innovation Representation Learning

Beyond the Gaussian and additive Gaussian models, there is no known algorithm
to obtain WIAE, especially when the underlying time series is nonparametric
with an unknown probability structure. In [32], the authors proposed a GAN-
based learning of the strong innovation representation by jointly minimizing
the Wasserstein distance of the latent process from the uniform IID process
and the mean squared error ($l_{2}$ distance) of the autoencoder output as the
estimate of the input. However, strong innovation representation applies only
to a restricted class of time series, and the joint optimization of
autoencoder with mixed Wasserstein and $l_{2}$ distance measures can be
challenging. Finally, learning a scalar innovation representation limits the
ability to incorporate multiple time series observations. The WIAE learning
proposed below overcomes these shortcomings.

###  II-C WIAE Learning

We present a deep learning approach to learn a WIAE for the weak innovation
representation defined in (2). Shown in Fig. 2 is the schematic that
highlights key components of the WIAE learning.

Figure 2: Training schematics of WIAE. Dash lines indicate the flow of
training information.

The encoder $G_{\theta}$ and decoder $H_{\eta}$ are causal convolutional
neural networks parameterized by coefficients $\theta$ and $\eta$,
respectively. The weak innovation representation, at its core, matches the
input-output distributions and constrains the latent process $({\bf V}_{t})$
to be IID-uniform. To this end, we introduce two neural network
discriminators, the innovation discriminator $D_{\gamma}$ and the
reconstruction discriminator $D_{\omega}$ with parameters $\gamma$ and
$\omega$ respectively, to enforce (.2) and (2). In particular, the innovation
discriminator $D_{\gamma}$ compares the distributions of $(\hat{{\bf V}}_{t})$
and $({\bf V}_{t})$, and the reconstruction discriminator $D_{\omega}$ the
compares joint distributions of ${\bm{X}}_{0:t+T}$ and
$({\bm{X}}_{0:t},\hat{{\bm{X}}}_{t+1:t+T})$ with sufficiently large $T$. These
discriminators produce error signals to update neural network parameters
$(\theta,\eta,\gamma,\omega)$. In this work, we adopt the Wasserstein
discriminator proposed in [33] to compute the Wasserstein distance between
distributions.

Because the two discriminators are both based on the Wasserstein-distance
measure, their parameters $(\theta,\eta,\omega,\gamma)$ can be obtained via a
single optimization:

$L:=\min_{\theta,\eta}\max_{\gamma,\omega}\big{(}\mathbb{E}[D_{\gamma}\left(({\bf
U}_{t})\right)]-\mathbb{E}[D_{\gamma}((\hat{{\bf V}}_{t}))]\\\
+\lambda(\mathbb{E}[D_{\omega}({\bm{X}}_{0:t+T})]-\mathbb{E}[D_{\omega}(({\bm{X}}_{0:t},\hat{{\bm{X}}}_{t+1:t+T}))])\big{)},$
(3)

where $\lambda$ is a real number that scales the two Wasserstein distances.
The two parts of the inner maximization loop of loss function (3) regularize
$(G_{\theta},H_{\eta})$ according to (.2) and (2). It’s evident that
minimizing the inner loop with respect to $\theta$ and $\eta$ is equivalent to
enforcing $({\bf V}_{t})$ being IID uniform, and
$({\bm{X}}_{0:t},\hat{{\bm{X}}}_{t+1:t+T})$ having the same distribution as
${\bm{X}}_{0:t+T}$. The training of the four neural networks is standard. Here
we used the off-the-shelf Adam optimizer.

In a practical implementation of WIAE, finite (input) dimensional neural
networks are used. The training of a finite-dimensional WIAE via (3) must also
be implemented by finite segments of the random processes involved. In Sec.
III-B, we consider the implications of such practical restrictions.

##  III WIAE-GPF and its Properties

In this section, we introduce WIAE-GPF. Specifically, given past observations
${\bm{x}}_{0:t}$, WIAE-GPF produces (arbitrarily many) samples of
$\hat{{\bm{X}}}_{t^{*}}$ that has the conditional distribution of
${\bm{X}}_{t^{*}}$. We present next the structure of WIAE-GPF, the Bayesian
sufficiency of WIAE, and a structure convergence when WIAE is implemented with
finite-dimensional implementations.

###  III-A Structure of WIAE-GPF

The structure of the proposed WIAE-GPF forecaster is shown in Fig. 3. At time
$t$, given the realization of ${\bm{X}}_{0:t}={\bm{x}}_{0:t}$ and autoencoder
$(G,H)$ trained by (3), ${\bm{x}}_{0:t}$ up to time $t$, the WIAE encoder $G$
generates the innovation sequence $\hbox{\boldmath$\nu$\unboldmath}_{0:t}$.
The WIAE decoder $H$ maps $\hbox{\boldmath$\nu$\unboldmath}_{0:t}$ and
independently generated IID-uniform pseudo innovations $\tilde{{\bf
V}}_{t+1:t^{*}}\stackrel{{\scriptstyle\mbox{\sf\tiny IID}}}{{\sim}}{\cal
U}[0,1]^{d}$ to produce a sample
$\tilde{{\bm{X}}}_{t^{*}}=\tilde{{\bm{x}}}_{t^{*}}$.

Note that in forecasting ${\bm{X}}_{t^{*}}$, we do not have realizations for
random samples of ${\bm{X}}_{t+1:t^{*}}$. The salient feature of WIAE-GPF is
to replace samples from the unknown and arbitrarily distributed
${\bm{X}}_{t+1:t^{*}}$ by realizations of pseudo innovations $\tilde{{\bf
V}}_{t+1:t^{*}}$ known to be IID-uniform. Note that once the autoencoder is
trained, generating random samples with the conditional distribution of
${\bm{X}}_{t^{*}}$ is trivial.

Figure 3: Forecasting pipeline for WIAE-GPF.

We now establish the validity of WIAE-GPF by showing that the WIAE-GPF output
$\tilde{{\bm{X}}}_{t^{*}}$ has the same conditional distribution as
${\bm{X}}_{t^{*}}$. This is not obvious because the input of $H$ includes an
exogenous random sequence $(\tilde{{\bf V}}_{t})$ and weak innovation $({\bf
V}_{0:t})$ that may not be a sufficient statistic.

We first show that the weak innovation sequence is Bayesian
sufficient‡‡‡$T(X)$ is Bayesian sufficient for the estimation of a random
variable $Y$ if the posterior distribution of $Y$ given $X$ is the same as the
one given $T(X)$ [34]., which implies that any stochastic decision involving
future time series ${\bm{X}}_{t^{*}}$ can be made without loss based on the
innovations ${\bf V}_{0:t}$. The same result was first presented in [30] under
the more restrictive setting of $H$ being injective.

######  Lemma 1 (Bayesian Sufficiency of Multivariate Weak Innovations)

Let $({\bm{X}}_{t})$ be a stationary time series for which the weak innovation
representation exists. Let $({\bf V}_{t})$ be the weak innovation
representation of $({\bm{X}}_{t})$. Then, for all ${\bm{x}}$ and
${\bm{X}}_{0:t}={\bm{x}}_{0:t}$,

$\displaystyle\Pr[{\bm{X}}_{t^{*}}\leq{\bm{x}}|{\bm{X}}_{0:t}={\bm{x}}_{0:t}]=\Pr\left[\hat{{\bm{X}}}_{t^{*}}\leq{\bm{x}}|{\bf
V}_{0:t}=\hbox{\boldmath$\nu$\unboldmath}_{0:t}\right].$ (4)

Proof: By the definition of weak innovation representation,

$\displaystyle\Pr[{\bm{X}}_{t^{*}}\leq{\bm{x}}|{\bm{X}}_{0:t}$
$\displaystyle={\bm{x}}_{0:t}]=\Pr[\hat{{\bm{X}}}_{t^{*}}\leq{\bm{x}}|{\bm{X}}_{0:t}={\bm{x}}_{0:t}]$
(5)
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\Pr[\hat{{\bm{X}}}_{t^{*}}\leq{\bm{x}}|G({\bm{X}}_{0:t})=G({\bm{x}}_{0:t}))]$
$\displaystyle=\Pr[\hat{{\bm{X}}}_{t^{*}}\leq{\bm{x}}|{\bf
V}_{0:t}=\hbox{\boldmath$\nu$\unboldmath}_{0:t}],$

where $(a)$ is from the Markovian structure of the autoencoder,i.e.,
${\bm{X}}_{0:t}\stackrel{{\scriptstyle G}}{{\rightarrow}}\hat{{\bf
V}}_{0:t}\stackrel{{\scriptstyle H}}{{\rightarrow}}\hat{{\bm{X}}}_{0:t}$. By
definition of Bayesian statistics [34], ${\bf V}_{0:t}=G({{\bm{X}}_{0:t}})$ is
a sufficient statistics for ${\bm{X}}_{t^{*}}$ for all $t^{*}>t$. $\square$

The validity of WIAE-GPF is shown next.

######  Theorem 1 (Validity of WIAE-GPF)

For all $t^{*}>t$, the conditional distribution of the WIAE-GPF output
$\tilde{{\bm{X}}}_{t^{*}}$ given ${\bm{X}}_{0:t}={\bm{x}}_{0:t}$ is identical
to that of ${\bm{X}}_{t+T}$ (given ${\bm{X}}_{0:t}={\bm{x}}_{0:t}$), i.e., ,

$\displaystyle\Pr[{{\bm{X}}}_{t^{*}}\leq{\bm{x}}|{\bm{X}}_{0:t}={\bm{x}}_{0:t}]=\Pr[\tilde{{\bm{X}}}_{t^{*}}\leq{\bm{x}}|{\bm{X}}_{0:t}={\bm{x}}_{0:t}].$
(6)

Proof: By Lemma 1,

$\displaystyle\Pr[{\bm{X}}_{t^{*}}\leq{\bm{x}}|{\bm{X}}_{0:t}={\bm{x}}_{0:t}]=\Pr[\hat{{\bm{X}}}_{t^{*}}\leq{\bm{x}}|{\bf
V}_{0:t}=\hbox{\boldmath$\nu$\unboldmath}_{0:t}],$ (7)

where ${\bf V}_{0:t}=G({\bm{X}}_{0:t})$ and
$\hbox{\boldmath$\nu$\unboldmath}_{0:t}=G({\bm{x}}_{0:t})$. Now consider

$\begin{cases}\hat{{\bm{X}}}_{t^{*}}=G({\bf V}_{0:t},{\bf V}_{t+1:t^{*}})\\\
\tilde{{\bm{X}}}_{t^{*}}=G({\bf V}_{0:t},\tilde{{\bf
V}}_{t+1:t^{*}}),\end{cases}$

where, by definition, $({\bf V}_{0:t},{\bf V}_{t+1:t^{*}},\tilde{{\bf
V}}_{t+1:t^{*}})$ are jointly independent IID uniform sequences, and
$\tilde{{\bf V}}_{t+1:t^{*}}\stackrel{{\scriptstyle\mbox{\tiny d}}}{{=}}{\bf
V}_{t+1:t^{*}}$. Therefore,

$\displaystyle\Pr[\hat{{\bm{X}}}_{t+T}\leq{\bm{x}}|{\bf
V}_{0:t}=\hbox{\boldmath$\nu$\unboldmath}_{0:t}]=\Pr[\tilde{{\bm{X}}}_{t+T}\leq{\bm{x}}|{\bf
V}_{0:t}=\hbox{\boldmath$\nu$\unboldmath}_{0:t}].$ (8)

Combining (7) and (8), we have (6). $\square$

###  III-B Structural Convergence of WIAE-GPF

This section focuses on the practical issue of finite-dimensional
implementations of WIAE and discriminators in Fig. 2. It is evident that no
machine learning technique guarantees that a finite-dimensional implementation
can extract weak innovations, even if the amount of historical samples
available is unbounded. Here we present a structural convergence result to
show that, under the ideal training conditions with unbounded training samples
and global convergence of training, the innovations generated from a finite-
dimensional WIAE converge in distribution the the true weak innovations.

The structural convergence analysis assumes that the training samples are
unbounded, and the training algorithm converges to the global minimum. Let
$G_{\theta^{*}}^{(m)}$ be the optimally trained finite (input) dimensional CNN
encoder that takes time-shifted $m$ consecutive observations
${\bm{X}}_{t-m+1:t}$ and produces the latent process $\left(\hat{{\bf
V}}^{(m)}_{t}\right)$. Likewise, let $H^{(m)}_{\eta^{*}}$ be the optimally
trained $m$-dimensional CNN decoder that produces the WIAE output sequence
$\left(\hat{{\bm{X}}}^{(m)}\right)$. Similarly defined are the finite
dimensional discriminators that take $n$ consecutive inputs, denoted by
$(D_{\omega}^{(n)},D_{\gamma}^{(n)})$. In this paper, we choose $n=m$.

To analyze the asymptotic property of finite (input) dimensional WIAE-GPF, we
make the following assumptions:

  1. A1

Existence: The random process $({\bm{X}}_{t})$ has a weak innovation
representation defined in (.1 \- .3) & (2), and there exists a causal WIAE
with continuous $G$ and $H$.

  2. A2

Feasibility: There exists a sequence of finite input dimension auto-encoder
functions $(G_{\bar{\theta}}^{(m)},H_{\bar{\eta}}^{(m)})$ that converges
uniformly to $(G,H)$ under the mean-squared distance metric.

  3. A3

Training: The training sample sizes are infinite. The training algorithm for
all finite-dimensional WIAE using finite-dimensional training samples
converges almost surely to the global optimum.

###### Theorem 2

Under (A1-A3),

$\displaystyle({\bm{X}}_{0:t},\hat{{\bm{X}}}_{t^{*}}^{(m)})\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{\rightarrow}}({\bm{X}}_{0:t},{\bm{X}}_{t^{*}}),~{}\forall t$ (9)

as $m$ goes to infinity.

Proof: See Sec. -A in the appendix.

###  III-C From GPF to Point and Quantile Forecasting

GPF produces samples of the conditional probability distribution, from which
point and quantile forecasts can be easily computed. Here we outline
techniques to compute forecasts of several popular point and quantile
forecasters. To this end, let
$\left\\{\hat{{\bm{x}}}_{t^{*}}^{(k)},k=1,\cdots,K\right\\}$ be the set of GPF
generated samples from the probability distribution of the time series at time
$t^{*}$ conditioned on past observations up to time $t$. For the simplicity of
mathematical expressions, we assume that
$\left\\{\hat{{\bm{x}}}_{t^{*}}^{(k)}\right\\}$ is sorted in the ascending
order.

  * •

Minimum Mean Squared Error (MMSE) Forecast: The MMSE forecast is the mean of
the conditional distribution. The MMSE forecast $\hat{{\bm{x}}}^{\mbox{\tiny
MMSE}}_{t^{*}}$ by a GPF is given by the conditional sample mean

$\hat{{\bm{x}}}^{\mbox{\tiny
MMSE}}_{t^{*}}=\frac{1}{K}\sum_{k=1}^{K}\hat{{\bm{x}}}_{t^{*}}^{(k)}.$

  * •

Minimum Mean Absolute Error (MMAE) Forecast: The MMAE forecast is the median
of the conditional distribution. The MMAE forecast
$\hat{{\bm{x}}}^{\mbox{\tiny MMAE}}_{t^{*}}$ by a GPF is given by the
conditional sample median

$\hat{{\bm{x}}}^{\mbox{\tiny
MMAE}}_{t^{*}}=\begin{cases}\hat{{\bm{x}}}_{t^{*}}^{\left((K+1)/2\right)},&\mbox{if
$K$ is odd}\\\
0.5\left(\hat{{\bm{x}}}_{t^{*}}^{\left(K/2\right)}+\hat{{\bm{x}}}_{t^{*}}^{\left(K/2+1\right)}\right),&\mbox{if
$K$ is even}.\\\ \end{cases}$

  * •

Quantile Forecast: The GPF forecast of $q$-quantile
$\hat{{\bm{x}}}^{\mbox{\tiny$q$-quantile}}_{t^{*}}$ is given by:

$\hat{{\bm{x}}}^{\mbox{\tiny$q$-quantile}}_{t^{*}}=\begin{cases}\hat{{\bm{x}}}_{t^{*}}^{\left(qK\right)},&\mbox{if
$qK$ is an integer}\\\
0.5\left(\hat{{\bm{x}}}_{t^{*}}^{\left([qK]\right)}+\hat{{\bm{x}}}_{t^{*}}^{\left([qK]+1\right)}\right),&\mbox{otherwise},\\\
\end{cases}$

where $[a]$ indicates the greatest integer not exceeding $a$.

##  IV WIAE-GPF for Market Operations

We conducted extensive experiments to compare leading GPF and point
forecasting techniques based on a suite of performance metrics. This section
summarizes our findings for three market applications where GPF is
particularly valuable to system operators and market participants: (a) LMP
forecasting for the optimal bidding in NYISO’s 5-minute real-time energy
market, (b) GPF for the interregional LMP spread for the Coordinated
Transaction Scheduling (CTS) [35] market between NYISO and PJM, and (c) ACE
forecasting for regulation services using PJM’s 15-second ACE data. Common to
these applications is that the forecasted variables are endogenously
determined by the market operations. In contrast to exogenous variables such
as wind/solar generations and inelastic demands, the LMP and ACE values are
the results of dispatch and commitment optimization, where binding constraints
introduce spikes in dual variables from which LMPs are computed. They are
highly dynamic as shown in Fig. 4.

Figure 4: Real-time, day-ahead LMPs, and load at Long Island, and real-time
LMP at NYC July 2023.

###  IV-A Baseline Methods in Comparison

TABLE II: Comparison of the baselines. Algorithm | Forecasting Type | Time Series Model | Forecastor Output | ML Models
---|---|---|---|---
SNARX [4] | Probabilistic | Semiparametric AR | AR Model Parameters | Kernel Estimation
WIAE-GPF | Probabilistic | Nonparametric | Generative | CNN + WIAE
TLAE [17] | Probabilistic | Parametric | Generative | RNN + VAE
DeepVAR [36] | Probabilistic | Parametric (AR Model) | Model Parameters | LSTM
BWGVT [16] | Probabilistic | Nonparametric | Forecasted Quantiles | LLM + Quantile Regression
Pyraformer [28] | Point | Nonparametric | Point Estimate | LLM
Informer [27] | Point | Nonparametric | Point Estimate | LLM

We compared WIAE-GPF with six leading forecasters based on their relevance to
power system applications and their established reputations. See Table II for
attributes of these techniques and references. WIAE-GPF is the only
nonparametric GPF forecaster. Because there are limited nonparametric GPF
techniques, we also included in our comparison popular machine-learning-based
parameterized GPF and point-forecasting techniques.

We started with SNARX [4], a classical parametric forecasting technique based
on an auto-regressive moving-average (ARMA) model. An early study showed that
SNARX performed the best among a list of parametric techniques [4]. For deep-
learning techniques, we compared with DeepVAR [36], a multivariate
generalization of the popular DeepAR that has become a key baseline for time
series forecasting in multiple applications. Temporal Latent AutoEncoder
(TLAE) [17] is an autoencoder-based parametric GPF where the conditional
distribution of the future time series variables is obtained through a
transformed Gaussian random process with mean and variance as parameters. Once
these parameters are estimated from observed realizations, Gaussian Monte
Carlo samples are fed into the decoder to generate samples of the forecasted
random variable.

We also included three popular forecasting techniques based on LLMs. One is
the award-winning technique Informer [27]; the other is Pyraformer [28] that
captures temporal dependency at multiple granularities§§§Pyraformer is the
keynote presentation at The International Conference on Learning
Representations (ICLR) 2022.. Pyraformer showed superior performance over a
wide range of LLM-based point forecasting techniques. Specifically developed
for LMP forecasting, BWGVT¶¶¶We have named the algorithm by the first letters
of authors’ last names [16] combines quantile regression with a transformer
architecture derived for LLMs.

###  IV-B Evaluation Metrics

Comparing probabilistic forecasting methods is difficult due to the lack of
ground truth for the underlying conditional distribution. However, because a
GPF can produce arbitrarily many Monte Carlo samples∥∥∥We used $1000$ Monte-
Carlo samples and sample average to obtain point estimates., it can be
evaluated by all point forecasting metrics. More importantly, an ideal
probabilistic forecaster that produces the correct conditional distribution
will perform well under any point estimator metric under regularity
conditions. Therefore, evaluating a GPF method based on a set of point-
forecasting techniques is a credible way to assess its performance.

We used four popular point forecasting and two widely used probabilistic
forecasting metrics. See Appendix for their definitions. Normalized mean
squared error (NMSE) measures the error associated with the mean of the
estimated conditional probability distribution. Normalized Absolute Error
(NMAE) measures the error associated with the median. Mean absolute scaled
error (MASE) is the ratio of the NMAE of a method over that of the (naive)
persistent predictor that uses the latest observation available as the
forecast. Symmetric mean absolute percentage error (sMAPE) averages and
symmetrizes the percentage error computed at each time stamp and is less
sensible to outliers. The probabilistic forecasting metrics were the
continuous ranked probability score (CRPS) [37], Coverage Probability Error
(CPE), and Normalized Coverage Width (NCW). CRPS evaluates the quadratic
difference between the predicted empirical cumulative density function
(c.d.f.) with an indicator c.d.f. based on the ground truth. CPE and NCW are
often used to evaluate prediction intervals. CPE is the deviation of the
coverage probability (CP) from the nominal confidence level $\beta\%$, whereas
NCW represents the width of the prediction intervals. At similar level of CP,
the method with smaller NCW shows better accuracy in prediction interval
estimation. In this paper, we computed the CPE and NCW of 10%, 50%, and 90%
intervals predicted by each probabilistic method. The mathematical definition
of those metrics can be found in the appendix.

For probabilistic forecasting techniques, we used their conditional means as
the point forecasts when evaluated by NMSE, whereas the conditional median is
used for NMAE, MASE, and sMAPE. For the quantile regression technique, BWGVT,
we use the estimated 0.5-quantile as its point forecast for all metrics since
it’s unclear how to compute empirical mean from quantiles. When producing
interval forecasts for GPF methods, we took the empirical quantiles from the
Monte-Carlo forecasts of future values. In particular, the beta-coverage
interval was defined as the $\beta$-width interval symmetric around the sample
median.

###  IV-C LMP forecasting for Energy Market Participation

TABLE III: Point evaluation of forecasting results for real-time price forecasting at Long Island. The numbers in the parentheses were the ranking of the algorithm. The columns under the label of LONGIL are the GPF performance of the 12-step foresting of the LMP at LONGIL. The columns under LONGIL & NYC are the performance of LMP forecasting at LONGIL using both LONGIL and NYC observations.  Methods | NMSE | NMAE | MASE | sMAPE
---|---|---|---|---
LONGIL | LONGIL & NYC | LONGIL | LONGIL & NYC | LONGIL | LONGIL & NYC | LONGIL | LONGIL & NYC
SNARX [4] | (7) $0.9852$ | (7) $0.5029$ | (7) $0.9733$ | (7) $0.7318$ | (7) $0.7053$ | (7) $0.4548$ | (7) $0.4319$ | (7) $0.4162$
WIAE-GPF | (1) $\mathbf{0.0585}$ | (2) $0.0487$ | (1) $\mathbf{0.2074}$ | (1) $\mathbf{0.1186}$ | (1) $\mathbf{0.1503}$ | (1) $\mathbf{0.0737}$ | (1) $\mathbf{0.0839}$ | (1) $\mathbf{0.0316}$
TLAE [17] | (4) $0.2956$ | (1) $\mathbf{0.0232}$ | (6) $0.4186$ | (2) $0.1366$ | (6) $0.3034$ | (2) $0.0849$ | (3) $0.2720$ | (2) $0.3190$
DeepVAR [36] | (6)$0.3919$ | (6) $0.4060$ | (5) $0.4088$ | (6) $0.4097$ | (5) $0.2963$ | (6) $0.2546$ | (5) $0.3629$ | (3) $0.3703$
BWGVT [16] | (3) $0.2670$ | (5) $0.2528$ | (4) $0.3158$ | (5) $0.3280$ | (4) $0.2289$ | (5) $0.2038$ | (4) $0.2817$ | (6) $0.3900$
Pyraformer [28] | (5) $0.3128$ | (3) $0.1382$ | (2) $0.3074$ | (3) $0.2556$ | (2) $0.2228$ | (3) $0.1588$ | (2) $0.1059$ | (4) $0.3765$
Informer [27] | (2) $0.1912$ | (4) $0.1829$ | (3) $0.3147$ | (4) $0.2830$ | (3) $0.2281$ | (4) $0.1759$ | (6) $0.3761$ | (5) $0.3827$
TABLE IV: Evaluation of probabilistic forecasting results of real-time prices
at Long Island. The numbers inside the square parentheses are the
corresponding NCW.

Methods | CRPS | CPE (90%) [NCW] | CPE (50%) [NCW] | CPE (10%) [NCW]
---|---|---|---|---
LONGIL | LONGIL & NYC | LONGIL | LONGIL & NYC | LONGLIL | LONGIL & NYC | LONGLIL | LONGIL & NYC
SNARX [4] | (3) $18.2864$ | (3) $10.5286$ | (5) $-0.2263$ $[0.6620]$ | (5) $-0.3779$ $[0.7990]$ | (4) $-0.2451$ $[0.5574]$ | (2) $-0.0746$ $[0.9161]$ | (3) $-0.0641$ $[0.8145]$ | (5) $0.1870$ $[0.8082]$
WIAE-GPF | (1) $\mathbf{4.6029}$ | (1) $\mathbf{1.1519}$ | (1) $~{}\mathbf{0.0245}$ $[0.9936]$ | (1) $~{}\mathbf{0.0157}$ $[0.9958]$ | (1) $\mathbf{0.0216}$ $[0.9681]$ | (1) $\mathbf{0.0184}$ $[0.8742]$ | (1) $\mathbf{0.0242}$ $[0.9836]$ | (1) $\mathbf{-0.0190}$ $[0.8207]$
TLAE [17] | (2) $6.1874$ | (2) ${2.8449}$ | (2) $~{}0.0325$ $[0.8179]$ | (4) $~{}0.0592$ $[0.8176]$ | (2) $0.0423$ $[0.6357]$ | (3) $0.1286$ $[0.7637]$ | (4) $-0.0848$ $[0.6704]$ | (4) $-0.0834$ $[0.6516]$
DeepVAR [36] | (4) $23.1460$ | (4) $22.9343$ | (4) $-0.0575$ $[0.9738]$ | (2) $-0.0592$ $[0.8467]$ | (3) $-0.1637$ $[0.5709]$ | (4) $-0.1488$ $[0.6614]$ | (2) $-0.0616$ $[0.5889]$ | (3) $-0.0431$ $[0.4260]$
BWGVT [16] | (5) $24.2595$ | (5) $24.3065$ | (3) $~{}0.0999$ $[7.3500]$ | (3) $~{}0.0829$ $[3.4864]$ | (5) $0.4815$ $[2.7822]$ | (5) $0.4322$ $[3.0100]$ | (5) $0.1034$ $[7.2325]$ | (2) $0.0311$ $[1.1181]$

(a) WIAE-GPF

(b) Pyraformer

(c) SNARX

(d) DeepVAR

Figure 5: Trajectories of the real-time price at LONGIL, and its prediction
generated by selective methods.

For a self-scheduled resource submitting a quantity bid to the energy market,
the ability to forecast future prices is essential in constructing its bids
and offers. With GPF generating future LMP realizations, the problem of
optimal offer/bid strategies can be formulated as scenario-based stochastic
optimization [38]. Our experiment was based on a use case of a merchant
storage owner submitting quantity offers and bids to a deregulated wholesale
market, using LMP from NYISO as the hypothetical price realizations.

The real-time market of NYISO closes sixty minutes ahead of actual delivery,
which means that the forecasting horizon needs to be longer than 60 minutes.
Two experiments were conducted to produce probabilistic forecasts of 60-minute
ahead LMPs at the Long Island (LONGIL) using (a) the day-ahead prices and the
current and past real-time LMPs at LONGIL, along with the system load up to
the time of submitting the bid; (b) the neighboring NYC real-time LMP in
addition to the data in (a).

Fig. 4 shows the real-time LMP trajectories at both LONGIL and NYC along with
the demand and the day-ahead LMP at LONGIL. The real-time LMPs at LONGIL and
NYC showed apparent spatial dependencies, while the dependency between day-
ahead and real-time LMPs at LONGIL. The dependencies between load and real-
time LMP were less obvious. The real-time LMPs and load were collected every 5
minutes and day-ahead LMPs every hour. We used the first 25 days for training
and validation, and last 6 days of July for evaluation.

Test results are shown in Table III and IV with boldfaced numbers being the
best performance. We observed that WIAE-GPF performed the best for all cases
except for the forecasting of LONGIL LMP using both LONGIL and NYC data, for
which WIAE-GPF is close second under NMSE. The strong performance of WIAE-GPF
comes from that it aims to match conditional distribution, with a validity
guarantee that ensures the Monte-Carlo samples generated have the same
conditional distribution as that of the actual time series variable. Overall,
the second-best is TLAE; the VAE-trained RNN autoencoder achieved noticeable
gain by combining data from two locations, especially for the joint
forecasting at LONGIL & NYC. However, under the probabilistic forecasting
measure of CRPS, TLAE was 50% to 100% higher than that of WIAE-GPF. Because
sampling the correlated Gaussian latent process is nontrivial, TLAE used a re-
parameterization heuristic that could cause an accumulation of biases.

The three LLM-based forecasters (BWGVT, Pyraformer, and Informer) performed
roughly the same, worse than WIAE-GPF and TLAE but better than SNARX and
DeepVAR. Note that Pyraformer and Informer are both point forecasters trained
to minimize the mean squared forecasting error. They did not perform well
under NMAE, MASE, and sMAPE. BWGVT is adapted to be a probabilistic
forecaster, although it did not perform well under the CRPS score. Its
quantile prediction was outlier-sensitive, especially when compared with the
GPF methods that depend on the stochastic latent process. We also observed
that BWGVT predicted larger intervals as shown by its high NCW for all cases.
Consequently, their CPE’s were always positive, indicating that the prediction
intervals it generated covered more than the nominal percentages. Hence, both
their point estimation results and CRPS were worse than TLAE and WIAE-GPF.

For the two parametric techniques DeepVAR and SNARX, being significantly
simpler than LLM forecasters, performed slightly worse. Both techniques did
not perform well under point estimation metrics. Shown by the CPE’s in Table
IV, their prediction interval estimation is also worse than other
probabilistic forecasting methods. Their NCW showed that they predicted narrow
intervals that failed to cover the majority of the ground truths. The
inaccurate semiparametric-AR model assumption appears to be culpable. Their
performance was influenced by model mismatch, where ARMA models typically
expect smoother trajectories. For the real-time LMP that exhibits high
volatility, these methods were slow catch the rapid changes, leading to
predicting shifted peaks that resulted in huge forecasting errors. This
phenomena were corroborated by their high sMAPE, indicating higher percentage
error at each time stamp.

It is interesting to observe that using LMPs in neighboring locations improved
the forecasting performance of all methods except DeepVAR and BWGVT, as shown
by the LONGIL and LONGIL & NYC columns under the DeepVAR and BWGVT rows in
Table III. This implies that these two methods didn’t account for spatial
correlations optimally. We conclude this to the difficulty of training of
transformer and LSTM models.

To gain insights into the performance of WIAE-GPF and other benchmark
techniques, we plotted the ground truth trajectories (black) and trajectory
forecasts generated by WIAE-GPF (red) and a competing algorithm (blue) in Fig.
5. Note that the spikes were not predicted by any methods. This was not
surprising given the nature of how these spikes were produced. Aside from
these spikes, these figures show clearly that WIAE-GPF (red) tracked the
ground truth (black) the closest, which was supported by the fact that WIAE-
GPF has the smallest NMAE. We also observed that WIAE-GPF had the smallest
variation, which is supported by the fact that WIAE-GPF had the smallest NMSE.
Furthermore, WIAE-GPF was the least affected by the price spikes. This was
because, as a GPF method, the Monte-Carlo samples used to produce the MMSE
point estimate were less likely to include extreme samples. However, for
SNARX, as shown in Fig. 5(c), the spikes caused significant deviations from
the forecasted trajectory.

###  IV-D Interregional LMP spread for Interchange Markets

The interchange market aims to improve overall economic efficiency across ISOs
by allowing virtual bidders to arbitrage price differences at proxy buses of
two neighboring ISOs. This experiment was based on the use case of a virtual
bidder bidding into the CTS market between NYISO and PJM. The proxy buses of
this market were Sandy Point of NYISO and Neptune of PJM.

The CTS market closes 75 minutes ahead of delivery and is cleared every 15
minutes. A virtual bidder submits a price-quantity bid along with the
direction of the virtual trade from the source of the proxy with low LMP to
the destination proxy with high LMP. Once the market is cleared, the
settlement is based on the actual LMP spread between the two proxies and the
cleared quantity. The bidder profits if the virtual trade direction matches
the direction of the real-time LMP spread. Otherwise, the bidder incurs a
loss. Therefore, the ability to predict the LMP spread direction is especially
important.

We performed a 75-minute ahead LMP spread forecasting using the interface
power flow and LMP spread data between NYISO and PJM at the Neptune proxy,
collected in February 2024. The interface power flow samples were collected
every 5 minutes, and LMP spread every 15 minutes. We used the first 24 days
for training and validation, and the last 5 days of February for testing.

We added Prediction Error Rate (PER) as a measure for the accuracy of the
virtual trading direction prediction, given that the sign of spread is of
great importance to profitability. PER indicates the percentage of forecasts
that don’t have the same direction as the ground truth. For point forecasts,
we compared the signs of the forecasts with the signs of the ground truth. For
probabilistic forecasting, we compared the direction of the ground truth with
that of the minimum error-probability prediction of the LMP spread, which is
the sign of the conditional median. For GPF, we compare the sample median with
the sign of the ground truth.

TABLE V: Evaluation of forecasting results for spread forecasting between
NYISO and PJM.

Methods | NMSE | NMAE | MASE | sMAPE | PER
---|---|---|---|---|---
SNARX [4] | (7) $2.4531$ | (7) $1.3415$ | (7) $1.1847$ | (7) $0.4958$ | (7) $0.6781$
WIAE-GPF | (1) $\mathbf{0.0098}$ | (1) $\mathbf{0.2738}$ | (1) $\mathbf{0.2418}$ | (1) $\mathbf{0.4493}$ | (1) $\mathbf{0.0606}$
TLAE [17] | (5) $0.9592$ | (5) $0.9785$ | (5) $0.8641$ | (4) $0.4785$ | (4) $0.3692$
DeepVAR [36] | (6)$1.8986$ | (3) $0.7224$ | (3) $0.6380$ | (5)$0.4806$ | (3) $0.3505$
BWGVT [16] | (3) $0.9053$ | (4) $0.8525$ | (4) $0.7529$ | (3) $0.4674$ | (2) $0.2313$
Pyraformer [28] | (4) $0.9478$ | (6) $1.2674$ | (6) $1.1193$ | (6) $0.4909$ | (6) $0.6738$
Informer [27] |  $0.8045$ (2) |  $0.4185$ (2) |  $0.4252$ (2) |  $0.4580$ (2) |  $0.5487$(5)

TABLE VI: Probabilistic forecasting results of spread forecasting between
NYISO and PJM.

Methods | CRPS | CPE (90%) [NCW] | CPE (50%) [NCW] | CPE (10%) [NCW]
---|---|---|---|---
SNARX [4] | (5) $120.0403$ | (5) $-0.5443$ $[5.0419]$ | (5) $-0.2958$ $[1.4555]$ | (2) $-0.0195$ $[0.6015]$
WIAE-GPF | (1) $\mathbf{4.0329}$ | (1) $\mathbf{0.0215}$ $[0.4427]$ | (2) ${-0.0274}$ $[0.5255]$ | (3) $0.0212$ $[0.1692]$
TLAE [17] | (2) $15.5195$ | (3) $-0.0443$ $[0.7745]$ | (1) $\mathbf{0.0052}$ $[0.8784]$ | (5) $-0.0933$ $[0.3133]$
DeepAR [36] | (4) $32.8296$ | (2) $-0.0355$ $[2.0279]$ | (3) ${-0.1480}$ $[1.4739]$ | (1) $\mathbf{-0.0021}$ $[0.4198]$
BWGVT [16] | (3) $31.5660$ | (4) $0.0989$ $[5.1788]$ | (4) $0.1835$ $[6.0939]$ | (4) $0.0656$ $[4.0473]$

(a) WIAE-GPF

(b) Pyraformer

(c) SNARX

(d) DeepVAR

Figure 6: Trajectories of the interregional LMP spread between NYISO and PJM,
and its prediction generated by selective methods.

Seen from Table V and VI, WIAE-GPF performed better than all other techniques
in all metrics. TLAE performed the second-best in CRPS ($15.5195$) but
slightly worse than BWGVT when evaluated under point estimation metrics. Its
sequential sampling of the latent Gaussian process added to its numerical
instability. BWGVT was the overall second-best performing probabilistic
technique. Its transformer architecture with enhanced capability of capturing
long-term temporal dependency didn’t offer much gain for the training
difficulty imposed by the increasing number of deep-learning parameters, see
Sec. IV-F. BWGVT also exhibited the tendency to predict a wide interval
covering more than the nominal percentage. Point estimation techniques, namely
Pyraformer and Informer, were not competitive when evaluated under point
estimation metrics other than NMSE. Among probabilistic methods, DeepVAR
performed similarly to the LLM methods, and SNARX had the most difficulties.
These (semi) parametric methods suffered from model mismatch, and were
sensitive to sudden changes, Thus, shifted peaks and valleys were often
witnessed in their predictions.

Same observation can also be made through Fig. 6. WIAE-GPF has the most stable
prediction of interregional LMP spreads, which is corroborated by its smallest
NMSE and NMAE. Pyraformer also follows the trend of LMP spreads accurately but
with higher variance. The two AR-based parametric models, SNARX and DeepVAR
exhibited the tendency to predict shifted spikes and failures to catch the
rapid and dramatic change of LMP spread.

###  IV-E Area Control Error Forecasting for Reserve Market Participants

ACE is defined as the difference between actual and scheduled load-generation
imbalance, adjusted by the area frequency deviation [39]. It is the control
signal for frequency regulation, and its probabilistic forecasting is
especially important for the operator to procure resources and market
participants to bid in the regulation ancillary service market.

In this subsection, we present the simulation results of a 5-minute ahead
forecasting of ACE. We utilized the ACE data from Jan 24th to 26th, collected
by PJM. The ACE signal is measured every 15 seconds and can be quite volatile,
as shown by the trajectory in Fig. 7.

Figure 7: Trajectory of ACE at PJM, Jan. 24th - 26th, 2023. TABLE VII: Estimation Results of ACE forecasting for PJM. The prediction step is 5-minute ahead. Methods | NMSE | NMAE | MASE | sMAPE
---|---|---|---|---
SNARX [4] | (6) $1.1922$ | (6) $1.1303$ | (6) $0.7029$ | (6) $0.4605$
WIAE-GPF | (1) $\mathbf{0.5957}$ | (1) $\mathbf{0.7555}$ | (1) $\mathbf{0.4698}$ | (1) $\mathbf{0.1059}$
TLAE [17] | (5) $1.1727$ | (5) $1.0605$ | (5) $0.6595$ | (3) $0.2782$
DeepVAR [36] | (7) $1.4431$ | (7) $1.1750$ | (7) $0.7307$ | (5) $0.3952$
BWGVT [16] | (3) $0.9562$ | (2) $0.9793$ | (2) $0.6090$ | (4) $0.3168$
Pyraformer [28] | (4) $0.9783$ | (4) $0.9948$ | (4) $0.6186$ | (7) $0.4986$
Informer [27] | (2) $0.6006$ | (3) $0.9819$ | (3) $0.6106$ | (2) $0.2247$
TABLE VIII: Probabilistic Estimation Results of ACE forecasting for PJM. The
prediction step is 5-minute ahead.

Methods | CRPS | CPE (90%) [NCW] | CPE (50%) [NCW] | CPE (10%) [NCW]
---|---|---|---|---
SNARX | (5) $2.1007$ | (5) $-0.8260$ $[1.6759]$ | (5) $-0.4797$ $[1.6803]$ | (2) $0.0343$ $[5.5768]$
WIAE-GPF | (1) $\mathbf{0.0081}$ | (1) $\mathbf{-0.0016}$ $[0.9199]$ | (1) $\mathbf{0.0321}$ $[0.9336]$ | (1) $\mathbf{-0.0132}$ $[0.8885]$
TLAE | (4) $1.5541$ | (4) $-0.7857$ $[0.0004]$ | (4) $-0.4489$ $[0.0005]$ | (4) $-0.0957$ $[0.0027]$
DeepAR | (3) $1.2947$ | (3) $-0.3526$ $[0.5665]$ | (3) $-0.2560$ $[0.5296]$ | (3) $-0.0521$ $[0.5434]$
BWGVT | (2) $1.2488$ | (2) $0.0065$ $[1.8309]$ | (2) $0.0754$ $[2.0385]$ | (5) $0.0996$ $[2.4261]$

WIAE-GPF achieved better performance than other methods, with CRPS less than
$0.01$ and sMAPE less than $11\%$, as shown by the WIAE-GPF row. We credited
the strong performance of WIAE-GPF to the simplicity of its latent process,
and its Bayesian sufficiency. BWGVT ranked second among all methods since the
ACE data has few outliers. But its CRPS at $1.2488$ is dramatically larger
than that of WIAE-GPF. Its CPE and NCW for 10% confidence interval prediction
also showed that it cannot accurately predict a narrow interval. Pyraformer
and Informer, trained with NMSE, had better performance under NMSE but worse
under NMAE. With NMSE over $110\%$, TLAE had the worst performance among GPF
methods. DeepVAR and SNARX performed worse than the other forecasting methods,
with NMSE and NMAE larger than $110\%$, possibly due to model mismatch.

###  IV-F Discussion: On using LGM

TABLE IX: Statistics that models long-range dependency of time series. Metrics | Real Time | Interchange Spread | ACE
---|---|---|---
Hurst Exponent | $0.5257$ | $0.5301$ | $0.5351$
DFA | $0.6053$ | $0.6614$ | $0.8609$

The success of LLM-based prediction in natural processing ignited broad
interest in adopting LLM models in various applications, including electricity
price forecasting with BWGVT, Pyraformer, and Informer. Our experiments showed
that the innovation-based method (WIAE) performed uniformly better than the
three LLM techniques, except for the real-time prediction at LONGIL-NYC under
NMSE, where Pyraformer was the best among all forecasters. Note that the
innovation representation used in WIAE can model but not explicitly long-range
dependencies of the random process. WIAE does not include attention modeling.

As LLM forecasters techniques, Pyraformer, Informer and BWGVT performed better
than the more conventional SNARX (see Table III, V, VII). Compared with the
more straightforward deep learning method of DeepAR, LLM-based forecasting did
not show clear advantages. The same can be said when they were compared with
TLAE.

To understand if long-range dependencies matter in the probabilistic
forecasting of electricity market signals, we examined the characteristics of
LMP signals using the Hurst exponent and Detrended Fluctuation Analysis (DFA)
as indicators for the long-range dependencies of LMP; both parameters had the
range [0,1], with deviation from 0.5 indicating symptoms of long-range
dependencies.

Table IX shows the estimated Hurst exponent and DFA. The Hurst Exponent and
DFA slope displayed a slight deviation from 0.5. The English and Korean
literature [40, 41] are known to have long-range dependencies with the Hurst
exponents ranging from 0.64 to 0.73. In comparison, the long-term effect of
real-time electricity market signals is minimal.

While further studies are necessary, the use of LLM may not be suitable for
electricity market signals where long-range dependencies are not evident.
Indeed, real-time LMPs are computed either on an interval-by-interval basis or
as part of short sliding window economic dispatch. Any temporal coupling is a
result of temporal dependencies of demand and supplies, neither shown to have
long-range dependencies. An unproven hypothesis is that the model complexity
of LLM may offset any benefit it may bring to price forecasting.

##  V Conclusion

This paper presents WIAE-GPF, a generative AI approach to probabilistic
forecasting of nonparametric time series based on the innovation
representation pioneered by Wiener, Kallianpur, and Rosenblatt six decades
earlier. Two contributions stand out. One is that the innovation
representation ensures that WIAE-GPF produces the correct conditional
probability distributions under perfect learning condition assumptions. To our
knowledge, WIAE-GPF is the first nonparametric GPF technique with such a
theoretical guarantee. Second, WIAE-GPF demonstrated superior performance
against traditional and major machine learning-based probabilistic forecasting
techniques in comprehensive numerical studies using real-world market data in
three market applications.

Somewhat surprisingly that direct implementations of LLM-based generative AI
techniques did not work well. One hypothesis is that the physics of generating
LMP and similar market signals do not have the long memory characteristics
that language models have because of the ways LMPs are computed, as we
demonstrated in [42]. Consequently, any benefits from LSTM and attention
mechanisms may be offset by complications of the training of significantly
more complex networks. Instead of using blindly LLM-based transformers,
considering how LMPs are generated would be more appropriate. In other words,
physics-based generative AI can potentially achieve significantly better
outcomes.

Although WIAE-GPF performed well for electricity market signals, the
underlying innovation representation as a canonical representation of time
series is a powerful model for other time series models. The application of
innovation representation to anomaly detection in power systems can be found
in [43, 44, 45].

###  -A Proof of Theorem 2

Let $\left(\bar{{\bf V}}_{t}^{(m)}\right)$ and
$\left(\bar{{\bm{X}}}_{t}^{(m)}\right)$ denote the latent process and the
reconstruction sequence, under weights $\bar{\theta}_{m}$ and $\bar{\eta}_{m}$

$\displaystyle\bar{{\bf
V}}_{t}^{(m)}=G_{\bar{\theta}}^{(m)}({\bm{X}}_{t},{\bm{X}}_{t-1},\cdots,{\bm{X}}_{t-m+1}),$
$\displaystyle\bar{{\bm{X}}}_{t}^{(m)}=H_{\bar{\eta}}^{(m)}(\bar{{\bf
V}}_{t},\bar{{\bf V}}_{t-1},\cdots,\bar{{\bf V}}_{t-m+1}).$

We define the loss of a WIAE pair $(G_{\theta},H_{\eta})$ achieved under a
$m$-dimensional discriminator pairs as

$L^{(m)}(\theta,\eta):=\max_{\gamma,\eta}\big{(}\mathbb{E}[D_{\gamma}^{(m)}\left({\bf
U}_{t:t-m+1}\right)]-\mathbb{E}[D_{\gamma}^{(m)}(\hat{{\bf V}}_{t:t-n+1})]\\\
+\lambda(\mathbb{E}[D_{\omega}^{(m)}({\bm{X}}_{t-n+2:t^{*}})]-\mathbb{E}[D_{\omega}^{(m)}(({\bm{X}}_{t-n+2:t},\hat{{\bm{X}}}_{t+1:t^{*}}))])\big{)}.$

We first show that $L^{(m)}(\theta_{m}^{*},\eta_{m}^{*})\rightarrow 0$ as
$m\rightarrow\infty$, where $(\theta_{m}^{*},\eta_{m}^{*})$ denotes the
optimal weights of $(G_{\theta}^{(m)},H_{\eta}^{(m)})$ obtained by minimizing
(3).

Following the line of [46], we defined the distance between two random
processes $({\bm{X}}_{t})$ and $({\bf Y}_{t})$ by the expected $\ell_{\infty}$
norm:

$d\left(({\bm{X}}_{t}),({\bf
Y}_{t})\right):=\mathbb{E}\left[\sup_{t}\lvert{\bm{X}}_{t}-{\bf
Y}_{t}\rvert\right].$

The uniform convergence assumed in assumption A2 is also defined on metric
spaces with distance measure $d(\cdot,\cdot)$. Hence, by assumption A2,
$G_{\bar{\theta}}^{(m)}\rightarrow G$ uniformly, which implies that,
$\forall\epsilon$, there exists a $M_{1}$ such that $\forall m>M_{1}$,
$d\left((\bar{{\bf V}}_{t}^{(m)}),({\bf V}_{t})\right)<\epsilon.$ Thus, for
$\forall F:\ell^{\infty}(T)\to\mathbb{R}$, $F$ bounded and continuous,

$d\left(F(({\bf V}_{t})),F((\bar{{\bf
V}}_{t}^{(m)}))\right)<\delta(\epsilon).$

In other words,

$\lim_{m\rightarrow\infty}\mathbb{E}[F((\bar{{\bf
V}}^{(m)}_{t}))]=\mathbb{E}[F(({\bf V}_{t}))],$

which fulfills the definition of weak convergence. Therefore,

$\displaystyle\bar{{\bf V}}_{t:t-m+1}^{(m)}\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{\rightarrow}}{\bf V}_{t:t-m+1},$ (10)

due to the fact that convergence in expectation implies convergence in
distribution.

Similarly, by the uniform convergence of $H_{\bar{\eta}}^{(m)}$ to $H$, we
have that that $\forall m>M_{2}$,

$\displaystyle d\left((\hat{{\bm{X}}}_{t}),(H_{\bar{\eta}}^{(m)}(({\bf
V}_{t})))\right)<\epsilon,$

where $(H_{\bar{\eta}}^{(m)}(({\bf V}_{t}))))$ represent the random sequence
generated by passing $({\bf V}_{t})$ through $H_{\bar{\eta}}$. Thus, for
$\forall F:\ell^{\infty}(T)\to\mathbb{R}$, $F$ bounded and continuous,

$d\left(F((\hat{{\bm{X}}}_{t})),(H_{\bar{\eta}}^{(m)}(({\bf
V}_{t}))))\right)<\delta(\epsilon).$

Hence we have $(H_{\bar{\eta}}^{(m)}(({\bf V}_{t}))))$ converges in
distribution to $(\hat{{\bm{X}}}_{t})$. Since $H$ is continuous and
$H_{\bar{\eta}}^{(m)}$ converges uniformly to $H$, $H_{\bar{\eta}}^{(m)}$ is
also continuous. Thus, by continuous mapping theorem,

$\bar{{\bf V}}_{t-m+1:t}^{(m)}\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{\rightarrow}}{\bf
V}_{t-m+1:t}\stackrel{{\scriptstyle}}{{\Rightarrow}}H_{\bar{\eta}}^{(m)}(\bar{{\bf
V}}_{t-m+1:t}^{(m)})\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{\rightarrow}}H_{\bar{\eta}}^{(m)}\left({\bf V}_{t-m+1:t}\right),$

that is,
$\bar{{\bm{X}}}_{t}^{(m)}\stackrel{{\scriptstyle\mbox{\tiny$d$}}}{{\rightarrow}}H_{\bar{\eta}}^{(m)}({\bf
V}_{t},\cdots,{\bf V}_{t-m+1})$. Therefore,

$\displaystyle({\bm{X}}_{t-m+2:t},\bar{{\bm{X}}}_{t^{*}}^{(m)})\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{\rightarrow}}({\bm{X}}_{t-m+2:t},\hat{{\bm{X}}}_{t^{*}})\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{=}}({\bm{X}}_{t-m+2:t},{\bm{X}}_{t^{*}}).$ (11)

By (10)&(11), $L^{(m)}(\bar{\theta}_{m},\bar{\eta}_{m})\rightarrow 0$. Since
$\theta^{*}_{m}$ and $\eta^{*}_{m}$ are the optimal parameters obtained by
minimizing (3) evaluated by $m$-dimensional discriminators
$(D_{\omega}^{(m)},D_{\gamma}^{(m)})$,

$\displaystyle
L^{(m)}(\theta_{m}^{*},\eta_{m}^{*}):=\min_{\theta,\eta}L^{(m)}(\theta,\eta)\leq
L^{(m)}(\bar{\theta}_{m},\bar{\eta}_{m})\rightarrow 0.$

Because $L^{(m)}(\theta_{m}^{*},\eta_{m}^{*})\rightarrow 0$ as
$m\rightarrow\infty$,
$\bm{V}_{t:t-m+1}^{(m)}\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{\rightarrow}}\bm{V}_{t:t-m+1}^{(m)}$ and
$({\bm{X}}_{t-m+2:t},\hat{\bm{X}}_{t^{*}}^{(m)})\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{\rightarrow}}\bm{(}{\bm{X}}_{t-m+2:t},\bm{X}_{t^{*}}^{(m)})$ follow
directly from the equivalence of convergence in Wasserstein distance and
convergence in distribution [47]. Since the discriminator dimensionality also
goes to $\infty$, we have
$({\bm{X}}_{0:t},\hat{{\bm{X}}}_{t^{*}}^{(m)})\stackrel{{\scriptstyle\mbox{\tiny
d}}}{{\rightarrow}}({\bm{X}}_{0:t},{\bm{X}}_{t^{*}})$. Further, the
conditional distribution of
$\hat{{\bm{X}}}_{t^{*}}^{(m)}|{\bm{X}}_{0:t}={\bm{x}}_{0:t}$ converges in
distribution to ${\bm{X}}_{t^{*}}|{\bm{X}}_{0:t}={\bm{x}}_{0:t}$ follows from
a simple application of the Bayes rule. $\square$

###  -B Definition of Metrics for Time Series Forecasting

Given the original time series $({\bm{x}}_{t})$, the forecasts
$(\tilde{{\bm{x}}}_{t})$, $N$ the size of datasets, and $T$ the prediction
step, the point estimation metrics can be calculated through:

$\displaystyle\mbox{NMSE}=\frac{\frac{1}{N}\sum_{t=1}^{N}({\bm{x}}_{t}-\tilde{{\bm{x}}}_{t})^{2}}{\frac{1}{N}\sum_{t=1}^{N}{\bm{x}}_{t}^{2}},$
$\displaystyle\mbox{NMAE}=\frac{\frac{1}{N}\sum_{t=1}^{N}|{\bm{x}}_{t}-\tilde{{\bm{x}}}_{t}|}{\frac{1}{N}\sum_{t=1}^{N}|{\bm{x}}_{t}|},$
$\displaystyle\mbox{MASE}=\frac{\frac{1}{N}\sum_{t=1}^{N}|{\bm{x}}_{t}-\tilde{{\bm{x}}}_{t}|}{\frac{1}{N-T}\sum_{t=T+1}^{N}|{\bm{x}}_{t}-{\bm{x}}_{t-T}|},$
$\displaystyle\mbox{sMAPE}=\frac{1}{N}\sum_{t=1}^{N}\frac{|{\bm{x}}_{t}-\tilde{{\bm{x}}}_{t}|}{(|{\bm{x}}_{t}|+|\tilde{{\bm{x}}}_{t}|)/2}.$

The purpose of adopting multiple metrics is to comprehensively evaluate the
forecasting performance. NMSE and NMAE evaluate the overall performance, and
MASE reflects the relative performance to the naive forecaster. Methods with
MASE smaller than $1$ outperform the naive forecaster. sMAPE is the symmetric
counterpart of mean absolute percentage error (MAPE) that can be both upper
bounded and lower bounded. Since for electricity datasets, the actual values
can be very close to $0$, thus nullifies the effectiveness of MAPE, we regard
sMAPE as the better metric.

For probabilistic methods, we further evaluates their CRPS. CRPS can be
computed from

$\displaystyle\mbox{CRPS}=\frac{1}{N-T}\sum_{t=1}^{N-T}\left(\int_{\mathbb{R}}\left(\tilde{F}_{t+T|t}({\bm{x}}|{\bm{x}}_{0:t})-\mathbb{I}_{{\bm{x}}_{t+T}\leq{\bm{x}}}\right)^{2}d{\bm{x}}\right),$

where $\mathbb{I}$ is the indicator function and
$\tilde{F}_{t+T|t}({\bm{x}}|{\bm{x}}_{0:t})$ the empirical cumulative density
function (c.d.f.) of $\tilde{{\bm{X}}}_{t+T}$ conditioned on
${\bm{X}}_{0:t}-{\bm{x}}_{0:t}$ predicted by probabilistic forecasting
methods. CRPS is equivalent to comparing the empirical conditional c.d.f.
forecasted by probabilistic methods with the indicator c.d.f.
$\mathbb{I}_{\tilde{{\bm{x}}}_{t}>{\bm{x}}_{t}}$ of the true value
${\bm{x}}_{t}$. It can be viewed as a generalization of MAE to probabilistic
methods.

The coverage probability (CP) of an confidence interval predictor is the
(estimated) probability that the ground truth falls within the predicted
interval. For a $T$-step prediction of $\beta\%$-intervals, we denote the
upper and lower bound by $\hat{U}_{t+T|t,\beta}$ and $\hat{L}_{t+T|t,\beta}$.
CP can be computed through

$\displaystyle\mbox{CP}(\beta\%)=\frac{1}{N-T}\sum_{t=1}^{N-T}\mathbb{I}_{{\bm{x}}_{t+T}\in[\hat{L}_{t+T|t,\beta},\hat{U}_{t+T|t,\beta}]}.$

The closer the CP to its nominal value $\beta\%$, the more accurate the
prediction is. Thus, the coverage probability error (CPE) is often adopted for
evaluation. CPE measures the deviation of CP from its nominal value $\beta\%$

$\mbox{CPE}(\beta\%)=\mbox{CP}(\beta\%)-\beta\%.$

The value of CPE closer to zero means the prediction interval estimation is
more accurate.

Although CP and CPE are widely adopted for its simplicity, since they only
estimate the unconditional coverage, they do not measure the accuracy the
coverage based on the forecasted conditional probability distribution. Its
limitation was discuss in [48].

In particular, while a good forecaster produces small CPE and a forecaster
with high CPE must be a poor forecaster, a forecaster producing small CPE may
not be a good forecaster. To this end, the normalized coverage width (NCW) can
be used as a secondary measure. NCW is defined as

$\displaystyle\mbox{NCW}(\beta\%)=\frac{1}{N-T}\sum_{t=1}^{N-T}\frac{\hat{U}_{t+T|t,\beta}-\hat{L}_{t+T|t,\beta}}{\hat{U}_{\beta}-\hat{L}_{\beta}},$

where $\hat{U}_{\beta}$ and $\hat{L}_{\beta}$ are the prediction interval
estimated from the empirical quantile of the testing data. For instance, when
predicting a $90\%$ interval, $\hat{U}_{90}$ is the empirical $0.95$-quantile
of the testing set, whereas $\hat{L}_{90}$ is the empirical $0.05$-quantile of
the testing set. As a result, NCW is the average width of intervals predicted
normalized by the width of the interval estimated through the empirical
marginal distribution of the testing set. One would expect that, conditional
on observations, one would get a more concentrated prediction interval than
the interval estimated based on unconditional distribution. Hence, a method
with NCW smaller than $1$ estimates prediction interval more accurately than
the unconditional estimation. At similar level of CP, the method with smaller
NCW shows better accuracy in prediction interval estimation.

## References

  * [1] H. GM, M. K. Gourisaria, M. Pandey, and S. S. Rautaray, “A comprehensive survey and analysis of generative models in machine learning,” _Computer Science Review_ , vol. 38, p. 100285, 2020.
  * [2] W. Härdle, H. Lütkepohl, and R. Chen, “A review of nonparametric time series analysis,” _International Statistical Review / Revue Internationale de Statistique_ , vol. 65, no. 1, pp. 49–72, 1997, publisher: [Wiley, International Statistical Institute (ISI)].
  * [3] J. Nowotarski and R. Weron, “Computing electricity spot price prediction intervals using quantile regression and forecast averaging,” _Computational Statistics_ , vol. 30, no. 3, pp. 791–803, Sep. 2015.
  * [4] R. Weron and A. Misiorek, “Forecasting spot electricity prices: A comparison of parametric and semiparametric time series models,” _International Journal of Forecasting_ , vol. 24, no. 4, pp. 744–763, 2008.
  * [5] M. Zhou, Z. Yan, Y. X. Ni, G. Li, and Y. Nie, “Electricity price forecasting with confidence-interval estimation through an extended arima approach,” _IEE Proc.-Gener.Transmiss.Distrib_ , vol. 153, no. 2, pp. 187–195, 2006.
  * [6] J. P. González, A. M. S. Muñoz San Roque, and E. A. Pérez, “Forecasting functional time series with a new hilbertian armax model: Application to electricity price forecasting,” _IEEE Transactions on Power Systems_ , vol. 33, no. 1, pp. 545–556, 2018.
  * [7] B. Uniejewski and R. Weron, “Regularized quantile regression averaging for probabilistic electricity price forecasting,” _Energy Economics_ , vol. 95, p. 105121, 2021.
  * [8] L. M. Lima, P. Damien, and D. W. Bunn, “Bayesian predictive distributions for imbalance prices with time-varying factor impacts,” _IEEE Transactions on Power Systems_ , vol. 38, no. 1, pp. 349–357, 2023.
  * [9] S. Chai, Z. Xu, and Y. Jia, “Conditional density forecast of electricity price based on ensemble elm and logistic emos,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 3, pp. 3031–3043, 2019.
  * [10] G. Dudek, “Multilayer perceptron for GEFCom2014 probabilistic electricity price forecasting,” _International Journal of Forecasting_ , vol. 32, no. 3, pp. 1057–1060, Jul. 2016.
  * [11] D. Lee, H. Shin, and R. Baldick, “Bivariate probabilistic wind power and real-time price forecasting and their applications to wind power bidding strategy development,” _IEEE Transactions on Power Systems_ , vol. 33, no. 6, pp. 6087–6097, 2018.
  * [12] J. Nowotarski and R. Weron, “Recent advances in electricity price forecasting: A review of probabilistic forecasting,” _Renewable and Sustainable Energy Reviews_ , vol. 81, pp. 1548–1568, Jan. 2018.
  * [13] D. J. Sheskin, _Handbook of Parametric and Nonparametric Statistical Procedures, Fifth Edition_. Chapman and Hall/CRC., 2011.
  * [14] C. Zhang and Y. Fu, “Probabilistic electricity price forecast with optimal prediction interval,” _IEEE Transactions on Power Systems_ , pp. 1–10, 2023.
  * [15] J.-F. Toubeau, T. Morstyn, J. Bottieau, K. Zheng, D. Apostolopoulou, Z. De Grève, Y. Wang, and F. Vallée, “Capturing spatio-temporal dependencies in the probabilistic forecasting of distribution locational marginal prices,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 3, pp. 2663–2674, 2021.
  * [16] J. Bottieau, Y. Wang, Z. De Grève, F. Vallée, and J.-F. Toubeau, “Interpretable transformer model for capturing regime switching effects of real-time electricity prices,” _IEEE Transactions on Power Systems_ , vol. 38, no. 3, pp. 2162–2176, 2023.
  * [17] N. Nguyen and B. Quanz, “Temporal Latent Auto-Encoder: A Method for Probabilistic Multivariate Time Series Forecasting,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 10, pp. 9117–9125, May 2021, number: 10.
  * [18] Z. Zheng, L. Wang, L. Yang, and Z. Zhang, “Generative Probabilistic Wind Speed Forecasting: A Variational Recurrent Autoencoder Based Method,” _IEEE Transactions on Power Systems_ , vol. 37, no. 2, pp. 1386–1398, Mar. 2022, conference Name: IEEE Transactions on Power Systems.
  * [19] L. Li, J. Zhang, J. Yan, Y. Jin, Y. Zhang, Y. Duan, and G. Tian, “Synergetic Learning of Heterogeneous Temporal Sequences for Multi-Horizon Probabilistic Forecasting,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 10, pp. 8420–8428, May 2021, number: 10.
  * [20] M. Khodayar, S. Mohammadi, M. E. Khodayar, J. Wang, and G. Liu, “Convolutional Graph Autoencoder: A Generative Deep Neural Network for Probabilistic Spatio-Temporal Solar Irradiance Forecasting,” _IEEE Transactions on Sustainable Energy_ , vol. 11, no. 2, pp. 571–583, Apr. 2020, conference Name: IEEE Transactions on Sustainable Energy.
  * [21] K. Rasul, A.-S. Sheikh, I. Schuster, U. M. Bergmann, and R. Vollgraf, “Multivariate Probabilistic Time Series Forecasting via Conditioned Normalizing Flows,” Feb. 2022.
  * [22] Y. Li, X. Lu, Y. Wang, and D. Dou, “Generative time series forecasting with diffusion, denoise, and disentanglement,” in _Advances in Neural Information Processing Systems_ , S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, Eds., vol. 35. Curran Associates, Inc., 2022, pp. 23 009–23 022.
  * [23] A. Koochali, P. Schichtel, A. Dengel, and S. Ahmed, “Probabilistic Forecasting of Sensory Data With Generative Adversarial Networks – ForGAN,” _IEEE Access_ , vol. 7, pp. 63 868–63 880, 2019, conference Name: IEEE Access.
  * [24] K. Yeo, Z. Li, and W. Gifford, “Generative Adversarial Network for Probabilistic Forecast of Random Dynamical Systems,” _SIAM Journal on Scientific Computing_ , vol. 44, no. 4, pp. A2150–A2175, Aug. 2022, publisher: Society for Industrial and Applied Mathematics.
  * [25] Y. Li, Y. Ding, Y. Liu, T. Yang, P. Wang, J. Wang, and W. Yao, “Dense skip attention based deep learning for day-ahead electricity price forecasting,” _IEEE Transactions on Power Systems_ , vol. 38, no. 5, pp. 4308–4327, 2023.
  * [26] H. Xu, F. Hu, X. Liang, and M. A. Gunmi, “Attention mechanism multi-size depthwise convolutional long short-term memory neural network for forecasting real-time electricity prices,” _IEEE Transactions on Power Systems_ , pp. 1–12, 2024.
  * [27] H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang, “Informer: Beyond efficient transformer for long sequence time-series forecasting,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 12, pp. 11 106–11 115, May 2021.
  * [28] S. Liu, H. Yu, C. Liao, J. Li, W. Lin, A. X. Liu, and S. Dustdar, “Pyraformer: Low-Complexity Pyramidal Attention for Long-Range Time Series Modeling and Forecasting,” Mar. 2022.
  * [29] N. Wiener, _Nonlinear Problems in Random Theory_. Cambridge, MA: Technology Press of Massachusetts Institute of Technology, 1958.
  * [30] X. Wang, L. Tong, and Q. Zhao, “Generative probabilistic time series forecasting and applications in grid operations,” _to appear in the Proceedings of Conference on Information Sciences and Systems_ , 2024. [Online]. Available: https://arxiv.org/abs/2402.13870
  * [31] M. Rosenblatt, “Stationary Processes as Shifts of Functions of Independent Random Variables,” _Journal of Mathematics and Mechanics_ , vol. 8, no. 5, pp. 665–681, 1959.
  * [32] X. Wang and L. Tong, “Innovations Autoencoder and its Application in One-class Anomalous Sequence Detection,” _Journal of Machine Learning Research_ , vol. 23, no. 49, pp. 1–27, 2022.
  * [33] M. Arjovsky, S. Chintala, and L.Bottou, “Wasserstein GAN,” Jan. 2017, arXiv:1701.07875.
  * [34] P. J. Bickel and K. A. Doksum, _Mathematical statistics: basic ideas and selected topics. (2nd ed.)_. Upper Saddle River, N.J.: Pearson Prentice Hall, 2007, vol. 1.
  * [35] M. White, R. Pike, C. Brown, R. Coutu, B. Ewing, S. Johnson, and C. Mendrala, “White paper: Inter-regional interchange scheduling analysis and options,” ISO New England and New York ISO, Tech. Rep., January 2011.
  * [36] D. Salinas, M. Bohlke-Schneider, L. Callot, R. Medico, and J. Gasthaus, “High-dimensional multivariate forecasting with low-rank gaussian copula processes,” in _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds., vol. 32. Curran Associates, Inc., 2019.
  * [37] T. Gneiting and M. Katzfuss, “Probabilistic forecasting,” _Annual Review of Statistics and Its Application_ , vol. 1, no. 1, pp. 125–151, 2014.
  * [38] E. Tómasson, M. R. Hesamzadeh, and F. A. Wolak, “Optimal offer-bid strategy of an energy storage portfolio: A linear quasi-relaxation approach,” _Applied Energy_ , vol. 260, p. 114251, 2020.
  * [39] NERC, “Balancing and frequency control,” NERC Resource Subcommittee, Priceton,NJ, Tech. Rep., January 2011. [Online]. Available: https://www.nerc.com/comm/OC/BAL0031_Supporting_Documents_2017_DL/NERC%20Balancing%20and%20Frequency%20Control%20040520111.pdf
  * [40] M. A. Montemurro and P. A. Pury, “Long-range fractal correlations in literary corpora,” _Fractals_ , vol. 10, no. 4, pp. 451–461, 2002.
  * [41] J. Bhan, S. Kim, J. Kim, Y. Kwon, S. il Yang, and K. Lee, “Long-range correlations in korean literary corpora,” _Chaos, Solitons & Fractals_, vol. 29, no. 1, pp. 69–81, 2006.
  * [42] X. Wang and L. Tong, “Generative probabilistic forecasting with applications in market operations,” 2024, arXiv:2403.05743. [Online]. Available: https://arxiv.org/abs/2403.05743
  * [43] ——, “Innovations autoencoder and its application in one-class anomalous sequence detection,” _J. Mach. Learn. Res._ , vol. 23, no. 1, Jan 2022.
  * [44] K. R. Mestav, X. Wang, and L. Tong, “A deep learning approach to anomaly sequence detection for high-resolution monitoring of power systems,” _IEEE Transactions on Power Systems_ , vol. 38, no. 1, pp. 4–13, 2023.
  * [45] L. Tong, X. Wang, and Q. Zhao, “Grid monitoring and protection with continuous point-on-wave measurements and generative ai,” 2024, arXiv:2403.06942. [Online]. Available: https://arxiv.org/abs/2403.06942
  * [46] J. Hoffmann-Jørgensen, _Stochastic Processes on Polish Spaces_. Aarhus, Denmark: Aarhus Universitet, Matematisk Institut., 1991.
  * [47] C. Villani, _The Wasserstein distances_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2009, pp. 57–75.
  * [48] P. F. Christoffersen, “Evaluating Interval Forecasts,” _International Economic Review_ , vol. 39, no. 4, pp. 841–862, Nov. 1998.