11institutetext: Porter School of the Environment and Earth Sciences, Raymond
and Beverly Sackler Faculty of Exact Sciences, Tel Aviv University, Tel Aviv,
6997801, Israel
11email<EMAIL_ADDRESS>22institutetext: Department of Particle Physics
and Astrophysics, Weizmann Institute of Science, Rehovot 7610001, Israel

# Adaptation of the phase distance correlation periodogram to account for
measurement uncertainties

A. Binnenfeld Adaptation of the phase distance correlation periodogram to
account for measurement uncertaintiesAdaptation of the phase distance
correlation periodogram to account for measurement uncertainties S. Shahaf
Adaptation of the phase distance correlation periodogram to account for
measurement uncertaintiesAdaptation of the phase distance correlation
periodogram to account for measurement uncertainties S. Zucker Adaptation of
the phase distance correlation periodogram to account for measurement
uncertaintiesAdaptation of the phase distance correlation periodogram to
account for measurement uncertainties

(Accepted XXX. Received YYY)

We present an improvement of the phase distance correlation (PDC) periodogram
to account for uncertainties in the time-series data. The PDC periodogram
introduced in our previous papers is based on the statistical concept of
distance correlation. By viewing each measurement and its accompanying error
estimate as a probability distribution, we are able to use the concept of
energy distance to design a distance function (metric) between measurement-
uncertainty pairs. We used this metric as the basis for the PDC periodogram,
instead of the simple absolute difference. We demonstrate the periodogram’s
performance using both simulated and real-life data. This adaptation makes the
PDC periodogram much more useful, demonstrating it can be helpful in the
exploration of large time-resolved astronomical databases, ranging from Gaia
radial velocity and photometry data releases to those of smaller surveys, such
as APOGEE and LAMOST. We have made a public GitHub repository available, with
a Python implementation of the new tools available to the community.

###### Key Words.:

methods: data analysis – methods: statistical – binaries: general – planets
and satellites: detection

##  1 Introduction

Time-domain astronomy poses distinctive challenges, including sparse and
uneven sampling, heteroscedasticity, high-dimensionality, and large data
volumes. In practice, the astrophysical information obtained from an
astronomical observation is often reduced to some scalar quantity, such as the
radial velocity or photometric magnitude. The temporal behavior of these
observables helps shed light on the system’s physical properties.

One key aspect of astronomical time series analysis is the search for periodic
modulation in the observed data. A common approach to detecting periodicity is
to generate a periodogram, produced by scanning a grid of trial periods (or
frequencies) and ascribing each one with a score. The score quantifies the
plausibility that the data are indeed modulated with a given period. A more
thorough characterization of the observed object is usually initiated only if
a periodic modulation is identified. The appeal of this approach stems from
its conceptual simplicity and reduced computational costs.

The Lomb-Scargle (LS) periodogram is a widely used, weighted, least squares-
based method for periodicity detection in unevenly sampled data (Ferraz-Mello,
1981). It assumes that a sine wave modulation is able to explain most of the
variation (‘energy’) in the data and that the measurement error is additive,
uncorrelated (’white’), and Gaussian. If the underlying signal is sinusoid and
the noise is white and Gaussian, LS would be an optimal detection scheme,
based on the Neyman-Pearson hypothesis-testing theory. The logic behind its
application largely relies on the ability to approximate smooth periodic
functions by Fourier series, assuming that most energy will concentrate on one
fundamental mode (see VanderPlas, 2018).

However, the assumption that the periodic signal is well represented by a
single harmonic function can be invalid in certain cases. This holds
particularly true for light curves exhibiting the majority of photometric
variability patterns. As a consequence, several extensions to LS have been
developed, aiming to broaden the flexibility of the search. For example,
Fawcett (2006) used several harmonics, instead of the pure sinewave used in
the classical LS periodogram; another extension is attributed to Zechmeister &
Kürster (2009), where a constant offset term was incorporated to account for
non-zero baseline values. The latter approach is often referred to as the
generalized LS periodogram (GLS).

Nevertheless, some astrophysical phenomena may give rise to signals that
cannot be approximated with only a few harmonics. Some studies have shown that
GLS and other parametric periodograms tend to underperform in such cases
(Pinamonti et al., 2017; Zucker, 2018). Examples of such signals include the
radial-velocity curves of highly eccentric Keplerian orbits or sawtooth-like
photometric modulations of pulsating stars. Exoplanetary transits are another
well-known example of a type of signal that is poorly approximated by
sinusoids. Instead of increasing the number of harmonics in the template model
used for the search, the common approach is changing its functional shape
altogether. The classical transit detection method BLS (box-fitting least
Squares) looks for box-shaped signals in the data (Kovács et al., 2002; Panahi
& Zucker, 2021; Shahaf et al., 2022). As the photometric accuracy improved,
the assumed underlying model was further refined in various ways to employ
more realistic transit-shaped models (e.g. Mandel & Agol, 2002; Kipping,
2023), further improving the search sensitivity.

The approaches discussed above are parametric, as they rely on the ability to
model or approximate the functional shape of the signal. There are cases,
however, where the functional shape is unknown or requires significant
computational resources for the calculation. Several non-parametric
periodograms were developed to address such situations. These techniques are
designed to be more agnostic to the particular modulation shape, expanding the
reach of the search, sometimes at the cost of reduced statistical efficiency:

String-length techniques are a set of commonly used model-independent methods
(Clarke, 2002). They include, among others, the Lafler-Kinman method (Lafler &
Kinman, 1965), as well as the methods of Renson (1978) and Dworetsky (1983).
In all of those methods, we quantify the dependence between consecutive phase-
folded measurements by estimating the length of an imaginary string connecting
them.

The analysis of variance method (AoV; Schwarzenberg-Czerny, 1989) is another
non-parametric approach, which basically fits a periodic piecewise constant
function to the measurements. After dividing the measurements into phase
subsets (’bins’) based on the trial period, the variance within each subset is
compared to the variance of the entire dataset to determine the significance
of the periodic signal at each period. A method closely related to the AoV
method is the phase dispersion minimization method (PDM; Stellingwerf, 1978),
minimizing the dispersion of the phase-folded data sets to detect the most
probable period.

The phase distance correlation periodogram (PDC; Zucker 2018) quantifies the
statistical dependence between the observable and the phase at which the data
were taken, given the trial period. The statistical dependence is quantified
by the distance correlation (Székely et al., 2007). The PDC periodogram does
not assume a specific shape for periodic variability and outperforms GLS in
detecting some highly non-sinusoidal periodic signals (Zucker, 2018).

One weakness of the non-parametric methods discussed above is their inability
to consider measurement uncertainties. Error estimates are essential
components of all physical measurements (e.g. Barlow, 1995). Uncertainties can
reflect various noise sources, such as inherently unavoidable photon noise or
instrumental glitches, and their estimates are particularly fundamental for
understanding and analyzing astronomical measurements (see Andrae, 2010).

Uncertainties constitute a unique challenge for model-independent periodograms
such as the PDC. Without an assumed functional shape, it is not
straightforward to define, let alone quantify, the deviations of the data from
their expected value. As a result, uncertainties are often omitted from the
calculation in these cases. One notable exception is related to the AoV and
PDM methods mentioned earlier, where an assumed functional form does indeed
exist, specifically one that minimizes dispersion. In this case, the
minimization procedure can involve sampling the assumed distribution of the
noise to incorporate them into the procedure.

In this work, we present an improvement of the PDC periodogram that considers
measurement uncertainties. We do so by viewing each measurement paired with
its corresponding uncertainty as a probability distribution, and then
introducing a new metric that enables the calculation of distance correlation
and therefore the PDC periodogram.

In the next section, we present the details of adapting the PDC periodogram to
account for measurement uncertainties. In Sect. 3, we demonstrate the
performance of the improved periodogram using simulated data and its
application to the HARPS data of a planet-hosting star. We present our
conclusions in Sect. 4 and discuss the method and its potential for future
studies.

##  2 Incorporating measurement errors in PDC

###  2.1 A new metric

As mentioned above, the distance correlation is a way to quantify the
statistical dependence between two random variables (Székely et al., 2007).
One of its merits is that the two random variables in question do not need to
be of the same dimensionality (Lyons, 2013).

Distance correlation is reliant on providing each random variable with an
adequate metric to determine the pairwise distance of the sample (Lyons 2013).
Therefore, to allow the PDC periodogram to account for measurement
uncertainties, we must define a suitable metric to consider both the
measurements and their accompanying error estimates. To do so, we assume that
each measurement is sampled from a Gaussian distribution, with the nominal
measurement value as its mean ($\mu$) and the accompanying error estimate as
its standard deviation ($\sigma$). We define the distance between two
measurements as the ‘Energy Distance’ between their two corresponding
distributions (Székely & Rizzo, 2013). This metric satisfies the required
conditions, outlined by Rizzo & Székely (2016), for the distance correlation
to be applied as a measurement of statistical independence.

The distance between two measurements is taken as the energy distance between
the two probability distributions from which these values were sampled,
$F_{i}$ and $F_{j}$. The explicit expression for the energy distance is
(Székely & Rizzo, 2013):

$\mathcal{D}^{2}(F_{i},F_{j})=2E||X_{i}-X_{j}||-E||X_{i}-X_{i}^{\prime}||-E||X_{j}-X_{j}^{\prime}||\
,$ (1)

where $X_{i}$ and $X_{j}$ are independent random variables distributed
according to $F_{i}$ and $F_{j}$; $X_{i}^{\prime}$ and $X_{j}^{\prime}$ are
independent and identically distributed copies of $X_{i}$ and $X_{j}$;
$||\cdot||$ denotes the Euclidean norm; $E$ represents the expectancy-value.

The difference between two normally distributed random variables is also
normally distributed. Therefore, $X_{i}-X_{j}$ is a Gaussian random variable
with mean $\mu_{i}-\mu_{j}$ and variance $\sigma_{i}^{2}+\sigma_{j}^{2}$.
Similarly, $X_{i}-X_{i}^{\prime}$ and $X_{j}-X_{j}^{\prime}$ are zero-mean
Gaussian distributions with variances given by $2\sigma_{i}^{2}$ and
$2\sigma_{j}^{2}$. Finally, to obtain an expression for the energy distance,
we need to know the distribution of the random variable $||X_{i}-X_{j}||$.
Since we are considering random scalars, the last expression is a random
variable representing the absolute value of the difference between two
normally distributed variables, which follows a ‘folded normal distribution’,
which is the distribution of the absolute value of a normally distributed
random variable.

Leone et al. (1961) showed that the expectancy value of a folded normal
distribution is

$\sigma\sqrt{\frac{2}{\pi}}\exp\left({-\frac{\mu^{2}}{2\sigma^{2}}}\right)+\mu\,\operatorname{erf}\left(\frac{\mu}{\sqrt{2\sigma^{2}}}\right)\
,$ (2)

where $\mu$ and $\sigma^{2}$ represent the mean and variance of the underlying
Gaussian, whose absolute value was taken, and $\operatorname{erf}$ is the
Gauss error function, as follows:

$\operatorname{erf}\left(x\right)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp\left(-t^{2}\right)dt\
.$

With these results in hand, we can explicitly write the energy distance (Eq.
1) for two normally distributed random variables. For convenience, we define
$\sigma_{ij}^{2}\equiv{{\sigma_{i}^{2}+\sigma_{j}^{2}}}$. In these terms, Eq.
1 becomes

$\displaystyle\varepsilon^{2}_{ij}\,=$
$\displaystyle{\sqrt{\frac{8}{\pi}}}\,\sigma_{ij}\,\left[{\exp({-x^{2}})}+x\,\operatorname{erf}\left(x\right)-y\right]\
,$ (3)

where

$\displaystyle
x\equiv\,\frac{\mu_{i}-\mu_{j}}{\sqrt{2}\sigma_{ij}}\quad\text{and}\quad
y\equiv\,\frac{\sigma_{i}+\sigma_{j}}{\sqrt{2}\sigma_{ij}}\,.$

It is illuminating to closely examine the behavior of the expression in Eq. 3.
First, it can be shown to be always non-negative, as is required for a metric.
The part that depends on $x$ is an increasing function of the absolute
difference of the two means (as Fig. 1 demonstrates): it is larger or equal to
$1$ and asymptotically converges to $|x|$. The second term, $y$, reflects the
heteroscedasticity of the measurement pair. It is bounded between $2^{-1/2}$,
for extreme variance ratios, to $1$, in the case of equal variances. In
essence, the heteroscedasticity makes the energy-distance-based metric treat
equal-variance distributions as closer to one another. Finally, the term in
square brackets is scaled by $\sigma_{ij}$, which increases the distance
between distributions with larger variances. The constant multiplicative
factor $\sqrt{8/\pi}$ can be omitted in practical contexts.

Figure 1: Illustration of the functional form attained by the energy distance
between two Gaussians as a function of the normalized difference between their
means, $x$. For simplicity, the diagram shows the value $x$-dependent term in
Eq. 3, illustrating that it is bounded from below by $1$ and asymptotically
approaches $|x|$. The additional term of the energy-distance-based metric
imposes a constant offset to the curve presented in this figure.

Following Lyons (2013, corollary 3.18, therein), we take the square root of
the metric to guarantee it is of ‘strong negative type’. Exploring the
definition and intricacies of strong negative-type metric spaces is outside
the scope of this paper, but it is important to note this attribute is a
sufficient condition for a metric to be applicable for distance correlation
computation (for further details, see Székely et al., 2007; Lyons, 2013).

$\mathcal{D}(F_{i},F_{j})=\sqrt{\varepsilon_{ij}^{2}}\,.$ (4)

Thus we obtained a strong-negative-type metric, which we then used to
formulate an adaptation of the PDC periodogram, to use measurement
uncertainties (as described in the following section).

###  2.2 Phase-distance periodogram

Following Zucker (2018, 2019), we defined a distance matrix based on the
pairwise energy distance introduced above. Each entry of the distance matrix
represents the distance between the corresponding measurements, namely

$a_{ij}=\mathcal{D}(i,j)\ .$ (5)

Once a distance matrix is obtained, we apply $\mathcal{U}$-centering
correction (Székely & Rizzo, 2014),

$A_{ij}=\begin{cases}\begin{split}a_{ij}-\frac{1}{N-2}\sum\limits_{k=1}^{N}a_{ik}-\frac{1}{N-2}\sum\limits_{k=1}^{N}a_{kj}\\\\[8.0pt]
+\frac{1}{(N-1)(N-2)}\sum\limits_{k,l=1}^{N}a_{kl}\end{split}&\text{if $i\neq
j$ ,}\\\ \\\ 0&\text{if $i=j$ .}\end{cases}$ (6)

Here, we use $N$ to represent the number of measurements. The
$\mathcal{U}$-centering procedure is required to produce an unbiased estimator
of the distance correlation.

Thus far, we have discussed the distances between pairs of observables. To
construct a periodogram, we must construct a distance matrix representing the
pairwise phase distance for each trial period, $P$. We used the metric defined
by Zucker (2018) to do so. The phase distance between two measurements is

$\begin{split}b_{ij}&=\phi_{ij}(P-\phi_{ij})\ ,\end{split}$ (7)

where

$\phi_{ij}=(t_{i}-t_{j})\mod P$

is the phase difference between two measurements. The resulting distance
matrix is then $\mathcal{U}$-centered, in a process identical to the one
presented in Eq. 6, producing the $\mathcal{U}$-centered phase distance
matrix, denoted as $B_{ij}$.

The unbiased estimator of the distance correlation is calculated using

$D=\frac{\sum\limits_{ij}A_{ij}B_{ij}}{\sqrt{\sum\limits_{ij}A^{2}_{ij}\sum\limits_{ij}B^{2}_{ij}}}\
.$ (8)

Similar to the original version of the PDC periodogram, the proposed recipe is
computationally intensive and entails an $O(N^{2})$ calculation for each trial
period. This quadratic dependence may be reduced in future applications to
$\mathcal{O}(N\log N)$ by using newly developed fast techniques to compute the
distance correlation (e.g. Huo & Szekely, 2014; Chaudhuri & Hu, 2019).

Following Binnenfeld et al. (2023), we used a $\chi^{2}$ test to assess the
significance of periods detected by the PDC (Shen et al., 2022), producing its
false alarm probability (FAP). The $\chi^{2}$-based FAP is simple to compute
and spares the need for computationally heavy permutation tests, which we
previously used for this purpose (Binnenfeld et al., 2020, 2022). The FAP is
calculated as follows:

$FAP=1-F_{\chi_{1}^{2}-1}(N\cdot D)$ (9)

with $F$ being the $\chi^{2}$ distribution cumulative distribution function
(CDF), $N$ is the number of measurements, and $D$ defined in Eq. 8.

Figure 2: Simulated demonstration of an Earth-analog exoplanet detection. Top:
Simulated RV curve of a planet-hosting star, as described in Sect. 3.1. The
red line represents the orbital model, and the black circles represent the
simulated noisy measurements, together with the simulated error bars. Bottom:
Comparison of the two versions of the PDC periodogram: with (solid blue) and
without (dashed red) considering the error bars. The dashed vertical line
marks the frequency associated with the injected orbital period of $365$ days.
The lower dotted horizontal line corresponds to a FAP level of $10^{-2}$, and
the upper one to a level of $10^{-3}$.

##  3 Demonstration and Verification

After introducing the theoretical foundations of the method, we demonstrate
its performance using simulated and real datasets. To quantify its merit, we
produced receiver operating characteristic (ROC) curves (e.g. Fawcett, 2006).
The ROC is a diagnostic tool for assessing the performance of detection
schemes by comparing the rates of true- and false-positive detection as a
function of the threshold applied to the detection statistic.

###  3.1 Earth-analog exoplanet detection

As a first example, intended to demonstrate the capabilities and advantages of
the new method for instances involving non-homogeneous errors, we generated a
sinusoidal RV curve corresponding to some hypothetical exoplanet with a
circular orbit of $365$ days, and radial velocity semi-amplitude of
$10\,\mathrm{cm\,s}^{-1}$. The radial velocity curve was sampled in randomly
drawn $35$ epochs from a uniform distribution on an interval of $1000$ days.
We added noise to each measurement by randomly drawing from a normal
distribution centered around $0$, with standard deviations that were drawn
from an exponential distribution with a scale parameter of
$1\,\mathrm{m\,s}^{-1}$. The standard deviations were paired with the
corresponding simulated measurements as their ‘estimated’ error bars. We used
this procedure to produce a wide range of noise standard deviations in the
same dataset. The simulated RV curve, with the simulated error bars, can be
seen in Fig. 2 (top panel).

As can be seen in the bottom panel of Fig. 2, the peak corresponding to the
actual periodicity in the original PDC periodogram is not prominent enough
compared to other spurious peaks, while the improved version of PDC does show
a dominant peak at the expected annual frequency, with a FAP of $10^{-3}$.

###  3.2 Detecting eccentric orbits

Periodic RV signals related to eccentric orbits were shown to be more
challenging to detect (Pinamonti et al., 2017). Their diverse variability
patterns pose difficulties for model-dependent methods, and they were
therefore identified as preferable targets for the PDC periodogram (Zucker,
2018).

To demonstrate how the new adaptation of PDC improves the PDC detection
capabilities even further, we chose the most eccentric exoplanet currently
known: HD 20782 b, with $e=0.97\pm 0.01$ and an orbital period of $597$ days
(O’Toole et al., 2008; Kane et al., 2016). We used RVs from the recently
published HARPS radial velocity database (Trifonov et al., 2020), which
consists of fifteen years of HARPS RVs.

HARPS optical fibers were upgraded in May 2015, changing the instrumental
profile and introducing an RV offset between the pre- and post-upgrade RVs.
Trifonov et al. recommended treating the pre- and post-upgrade time series
from the database as taken from two different instruments, a recommendation we
indeed followed in our analysis. Out of $87$ measurements in total, we only
used the $72$ taken prior to the fiber change. The RVs and their error
estimates can be seen in the bottom panels of Fig. 3.

As seen in the top panel of Fig. 3, the GLS periodogram does not exhibit a
real prominent peak in the previously published frequency, while the regular
PDC version does. The significance of the peak is even emphasized in the
improved version of the periodogram. The change translates to a major decrease
in the detection FAP, from $10^{-9}$ to $10^{-13}$.

RV curves of very eccentric orbits can benefit significantly from the new
adaptation of the periodogram. As demonstrated by the case of HD 20782,
samples that occur during periastron passage can easily be mistaken for
outliers, even if their uncertainty estimates are small. In other cases, the
error bars might be more ambiguous. The new and improved version offers an
objective and quantitative way to include the error bars in the analysis

Figure 3: Demonstration of the advantages of the new version of the PDC
periodogram, which takes into account the uncertainties, using RV data of the
planet-host star HD 20782 b data. Top: GLS and PDC periodograms. The dashed
vertical line marks the frequency corresponding to the known orbital period of
$597$ days. The lower dotted horizontal line corresponds to a FAP level of
$10^{-9}$, and the upper one to a level of $10^{-13}$. Middle panels: HARPS
RVs measurements for HD 20782, and their estimated measurement uncertainties
(plotted separately for clarity). Bottom: RVs phase-folded according to the
previously published orbital period.

###  3.3 ROC curves

We ran the tests for the ROC curve using the same procedure described in Sect.
3.1, only on $1000$ sinusoidal RV curve realizations with periods drawn from a
uniform distribution between $7$ and $180$ days. Similarly, we drew a random
number of samples between $20$ and $70$. We added noise to each measurement by
randomly drawing from a normal distribution centered around $0$, with standard
deviations that were drawn from an exponential distribution with a scale
parameter of $1\,\mathrm{m\,s}^{-1}$. Half of the light curves were then
intentionally shuffled to eliminate the periodic signal, so that we could
later estimate the FP rate. We calculated the two periodograms for each
realization on the same frequency grid with 200 frequencies between $10^{-4}$
and $0.2\,\textrm{day}^{-1}$.

We then used the resulting periodograms to generate the ROC curve by sorting
the FAP thresholds of each periodogram in ascending order. By serially
computing the pairs of FP and TP rates, we created the ROC curves presented in
Fig. 4.

As can be seen in the figure, the new version of the PDC (solid blue line), as
it takes into account the uncertainties, outperforms the old version (dashed
red line).

As noted in the previous subsection, the PDC periodogram outperformed other
methods in detecting signals of eccentric orbits. To demonstrate this, we
repeated the ROC curve test, this time simulating $1000$ realizations of
eccentric Keplerian orbits (uniformly distributed $e>0.35$). As shown in Fig.
5, the PDC periodogram in its new version (using the information contained in
the error bars) does indeed outperform both the GLS and the regular version of
the PDC in this detection test.

The area under the ROC curve (AUC) metric is often used for analyzing ROC
curves, offering an aggregated performance measure, invariant to both scale
and detection threshold. In the context of the PDC periodogram, it can be
effectively employed to measure the method contribution, particularly when
dealing with data exhibiting diverse characteristics such as variance and
signal-to-noise ratio (S/N). Future studies, focusing on defining and
exploring the physical parameter space for different test cases, could
leverage the AUC metric to gain a more nuanced understanding of the
periodogram efficacy.

####  3.3.1 Contribution of the uncertainties information

To further illustrate the capabilities of the updated version of the PDC and
establish the importance of the contribution of measurement uncertainties to
this detection scheme, we created an additional ROC curve. We repeated the
$1000$ sinusoidal realization test described in the previous subsections. This
time, we evaluated the performance of the periodogram under two distinct
settings: the standard configuration and a scenario in which we randomly
shuffled the uncertainti estimated, highlighting their critical role in the
enhanced performance (presented in Fig. 6).

As is evident from the figure, the shuffled version is outperformed by the
regular PDC, while the newly developed method in its standard configuration is
better than both. Nonetheless, the fact that false error estimates do not
significantly impair the performance of the PDC periodogram is encouraging and
attests to the robustness of the PDC.

##  4 Conclusion

In this paper, we present a modification of the PDC periodogram that is aimed
at accounting for measurement uncertainties. We have demonstrated its
performance using both simulations and data from the recently published HARPS-
RVBank archive.

The original version of the PDC periodogram ignored error bars, which was a
real disadvantage of the PDC and is likely the reason the method was not
widely adopted. Now that we have found a way to incorporate into PDC the
information contained in the error bars, we believe we have made it much more
useful.

Although this study emphasizes RV curve illustrations, it is important to note
that the PDC periodogram and its current modification presented in this paper
are well suited for photometric light curves as well. As we have shown
(Zucker, 2018), the PDC periodogram outperforms other methods in various cases
of non-sinusoidal variability patterns, such as the sawtooth-like variability
shapes that are typical of pulsating stars.

The adaptation presented in this paper relies on a Gaussian model for the
uncertainty, which we used in order to construct the energy distance metric.
The Gaussian uncertainty model is used very often in many contexts (almost by
default) due to its simplifying properties. However, somewhat less often,
other uncertainty models are used (e.g., Poisson distribution). The energy
distance expression shown in Eqs. 3-4 should be replaced in those cases with
the relevant expressions suitable for the assumed noise distribution. This
approach may be useful in a range of applications, including the potential to
extend the definition of distance to encompass and mitigate correlated noise.

In recent papers (Binnenfeld et al., 2020, 2023), we introduced USuRPER, which
is essentially a PDC periodogram for spectral observations and also a version
of the PDC periodogram suitable for scanning astrometry. Despite it not being
suited to the exact modification presented in this paper (i.e. scalar
observations), we believe similar approaches of viewing measurements and their
errors as probability distributions can be applied for those periodograms as
well, to make them account for measurement error bars.

Additional relevant PDC extension is the ‘partial phase distance correlation
periodograms’, which allow for ‘nuisance’ parameters to be accounted for with
the aim of eliminating spurious peaks related to those nuisance parameters
(Binnenfeld et al., 2022). In the case of two simultaneously measured scalar
quantities, such as RVs and some activity indicator, the latter may be
considered a nuisance; thus, it can be used to eliminate periodogram peaks
originating from stellar activity. The formulation presented in this work can
also be used to improve the partial periodogram by using the relevant
uncertainty estimates.

Finally, we provide our Python implementation of the periodogram in the form
of a public GitHub repository111PDC and its extensions, including USuRPER,
partial distance correlation periodograms, PDC periodogram for scanning
astrometry, and the periodogram presented in this work, are all available as
part of the SPARTA package (Shahaf et al., 2020), at
https://github.com/SPARTA-dev..

###### Acknowledgements.

We thank the editor and the referee, Andrej Prša for reviewing this paper and
for their enlightening and useful comments. This research was supported by the
Ministry of Innovation, Science & Technology, Israel (grant 3-18143) and by
the ISRAEL SCIENCE FOUNDATION (grant No. 1404/22). The research of SS is
supported by a Benoziyo prize postdoctoral fellowship. The analyses done for
this paper made use of the code packages: NumPy (Harris et al., 2020), SciPy
(Virtanen et al., 2020) and SPARTA (Shahaf et al., 2020).

Figure 4: ROC curves comparing the regular PDC periodogram with its newly
developed version. The dashed red line represents the regular PDC, and the
thick blue one represents the new version that considers uncertainty
estimates. The thin black line stands for a $1$:$1$ ratio between the true and
false positive rates, i.e., random guess. Figure 5: ROC curves for the case of
eccentric Keplerian RV curves ($e>0.35$). The dashed red line represents the
regular PDC, and the solid blue line represents the new version presented in
this paper. They both outperform the GLS periodogram, represented by the
dotted grey line. The thin black line stands for a $1$:$1$ random guess.
Figure 6: ROC curves comparing our newly developed method in its standard
configuration against a scenario where the uncertainty estimates were randomly
shuffled. The dashed red line represents the regular PDC, and the solid blue
line represents the new version. The green dotted line represents the
performance with the randomly shuffled uncertainty estimates. The thin black
line represents a random guess.

## References

  * Andrae (2010) Andrae, R. 2010, arXiv e-prints, arXiv:1009.2755
  * Barlow (1995) Barlow, R. 1995, Statistics: A Guide to the Use of Statistical Methods in the Physical Sciences, Manchester physics series (Wiley)
  * Binnenfeld et al. (2022) Binnenfeld, A., Shahaf, S., Anderson, R. I., & Zucker, S. 2022, A&A, 659, A189
  * Binnenfeld et al. (2023) Binnenfeld, A., Shahaf, S., Anderson, R. I., & Zucker, S. 2023, A&A, 659, A189
  * Binnenfeld et al. (2020) Binnenfeld, A., Shahaf, S., & Zucker, S. 2020, A&A, 642, A146
  * Chaudhuri & Hu (2019) Chaudhuri, A. & Hu, W. 2019, Computational Statistics & Data Analysis, 135, 15–24
  * Clarke (2002) Clarke, D. 2002, A&A, 386, 763
  * Dworetsky (1983) Dworetsky, M. M. 1983, MNRAS, 203, 917
  * Fawcett (2006) Fawcett, T. 2006, Pattern Recognition Letters, 27, 861, rOC Analysis in Pattern Recognition
  * Ferraz-Mello (1981) Ferraz-Mello, S. 1981, AJ, 86, 619
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Huo & Szekely (2014) Huo, X. & Szekely, G. J. 2014, Fast Computing for Distance Covariance
  * Kane et al. (2016) Kane, S. R., Wittenmyer, R. A., Hinkel, N. R., et al. 2016, The Astrophysical Journal, 821, 65
  * Kipping (2023) Kipping, D. 2023, MNRAS, 523, 1182
  * Kovács et al. (2002) Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369
  * Lafler & Kinman (1965) Lafler, J. & Kinman, T. D. 1965, ApJS, 11, 216
  * Leone et al. (1961) Leone, F. C., Nelson, L. S., & Nottingham, R. B. 1961, Technometrics, 3, 543
  * Lyons (2013) Lyons, R. 2013, Ann. Probab., 41, 3284
  * Mandel & Agol (2002) Mandel, K. & Agol, E. 2002, ApJ, 580, L171
  * O’Toole et al. (2008) O’Toole, S. J., Tinney, C. G., Jones, H. R. A., et al. 2008, Monthly Notices of the Royal Astronomical Society, 392, 641
  * Panahi & Zucker (2021) Panahi, A. & Zucker, S. 2021, PASP, 133, 024502
  * Pinamonti et al. (2017) Pinamonti, M., Sozzetti, A., Bonomo, A. S., & Damasso, M. 2017, MNRAS, 468, 3775
  * Renson (1978) Renson, P. 1978, A&A, 63, 125
  * Rizzo & Székely (2016) Rizzo, M. L. & Székely, G. J. 2016, Wiley Interdiscip. Rev. Comput. Stat., 8, 27
  * Schwarzenberg-Czerny (1989) Schwarzenberg-Czerny, A. 1989, Monthly Notices of the Royal Astronomical Society, 241, 153
  * Shahaf et al. (2020) Shahaf, S., Binnenfeld, A., Mazeh, T., & Zucker, S. 2020, Astrophysics Source Code Library, ascl:2007.022
  * Shahaf et al. (2022) Shahaf, S., Zackay, B., Mazeh, T., Faigler, S., & Ivashtenko, O. 2022, MNRAS, 513, 2732
  * Shen et al. (2022) Shen, C., Panda, S., & Vogelstein, J. T. 2022, J. Comput. Graph. Stat., 31, 254
  * Stellingwerf (1978) Stellingwerf, R. F. 1978, ApJ, 224, 953
  * Székely & Rizzo (2013) Székely, G. J. & Rizzo, M. L. 2013, J. Stat. Plan. Inference, 143, 1249
  * Székely & Rizzo (2014) Székely, G. J. & Rizzo, M. L. 2014, Ann. Stat., 42, 2382
  * Székely et al. (2007) Székely, G. J., Rizzo, M. L., & Bakirov, N. K. 2007, Ann. Stat., 35, 2769
  * Trifonov et al. (2020) Trifonov, T., Tal-Or, L., Zechmeister, M., et al. 2020, A&A, 636, A74
  * VanderPlas (2018) VanderPlas, J. T. 2018, ApJS, 236, 16
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261
  * Zechmeister & Kürster (2009) Zechmeister, M. & Kürster, M. 2009, A&A, 496, 577
  * Zucker (2018) Zucker, S. 2018, MNRAS, 474, L86
  * Zucker (2019) Zucker, S. 2019, MNRAS, 484, L14


# Joint Stochastic Model for Electric Load, Solar and Wind Power at Asset
Level and Monte Carlo Scenario Generation

René Carmona & Xinshuo Yang  Department of Operations Research & Financial
Engineering
Princeton University
Princeton, NJ 08544, USA

(Date: February 20, 2021)

###### Abstract.

For the purpose of Monte Carlo scenario generation, we propose a graphical
model for the joint distribution of wind power and electricity demand in a
given region. To conform with the practice in the electric power industry, we
assume that point forecasts are provided exogenously, and concentrate on the
modeling of the deviations from these forecasts instead of modeling the actual
quantities of interest. We find that the marginal distributions of these
deviations can have heavy tails, feature which we need to handle before
fitting a graphical Gaussian model to the data. We estimate covariance and
precision matrices using an extension of the graphical LASSO procedure which
allows us to identify temporal and geographical (conditional) dependencies in
the form of separate dependence graphs. We implement our algorithm on data
made available by NREL, and we confirm that the dependencies identified by the
algorithm are consistent with the relative locations of the production assets
and the geographical load zones over which the data were collected.

_Keywords:_ Load, Solar Power, Wind power, Monte Carlo simulations, Graphical
LASSO, Conditional Simulation.

##  1\. Introduction

Modern operation of the electric grid has become highly technical but despite
clear evidence of the stochasticity of electricity demand the major unit
commitment and economic dispatch computer algorithms remained intrinsically
deterministic. This is the more problematic given the increased penetration of
renewable energy production sources like solar and wind whose uncertainly and
volatility have been universally documented. For this reason, it became
commonplace to have engineers rely on probabilistic forecasts of wind, solar
and even loads in lieu of deterministic point forecasts of these quantities.
Even though these new classes of forecasts provide a form of quantification of
the variability of the variables, and hence of the reliability of the point
forecasts, these forecasts remain very limited in the sense that they have no
ability to capture temporal and spatial dependencies between the random
quantities they try to predict. Still they have been the object of intense
interest during the last decades and the sources of a very large number of
publications. For the methodologies behind these probabilistic forecasts we
refer to the following recent surveys: [7] for load, [21] for wind power, and
[19] for solar power generation. While the application domain is the same, our
approach differs significantly from those papers. Indeed, not only do we not
produce forecasts, but we use existing forecasts and we model the deviations
of the actual values from their forecasts. This choice is based on the fact
that many users have their own forecasts which they do not want to part from,
whether they purchase them or they develop them in house,

Beyond the interest in the production of forecasts, a large number of studies
were devoted to the development of algorithms to generate statistical
scenarios from probabilistic forecasts. For wind power, we refer to the
fundamental paper [14] and the references therein, and to [16] for a recent
_state of the art_ on the subject with a clear goal of providing inputs for
unit commitment and economic dispatch platforms. Evaluation of the quality of
scenarios produced this way was a natural development. See for example [13]
for the analysis of wind scenarios.

Unfortunately, spatial and temporal correlations were very rarely included in
the simulation models, and when it was done, it was often for one single
variable only. See [17] for an example dealing with wind or [18] for a spatio-
temporal model for wind power. One of our motivations, and possibly the main
thrust of our paper, is to propose algorithms to remedy this shortcoming by
allowing temporal and spacial dependencies to be captured for load, wind power
and solar power simultaneously, even at the most granular level of the
individual power plants. Modeling the random behavior of the future values of
a variable, say wind, is commonly done on the basis of existing forecasts
obtained exogenously. See for example [1]. This trend has become a standard in
the industry, and we shall follow this time honored convention in the present
paper.

###  1.1. Contributions of the paper.

This paper is a contribution to the literature on the modeling of the
stochastic behavior of load, wind and solar power at their most granular level
available, and the development of simulation algorithms capable to deliver
Monte Carlo scenarios for their future time evolutions. While often touted in
the published literature, the analysis of fully correlated models is very
rarely pushed beyond the level of ideas and discussions, and even more rarely
implemented by means of simulation algorithms. See nevertheless the inclusion
of time correlation for wind power in [10, 18]. We do not know of any prior
work involving joint spatio-temporal models for load, wind and solar.

Before delving into the gory details of our modeling approach and the inner
working of our simulation engines, we want to emphasize that for all the time
series of interest, we work with the deviations from forecasts. This is a time
honored approach which is widely accepted in the engineering community, and
for this reason, which we follow. The down side is that the historical data
needed to fit our model needs to include actual measurements of load, wind and
solar power as well as their respective forecasts. Our models and their
implementations as simulation engines are agnostic to the source of the
forecasts. However, the quality of the latter does impact the goodness of fit
of our models as well as the predictive power of the Monte Carlo simulator.

We now describe in broad brushstrokes the several prongs of our approach, in
no small part, to emphasize the originality of our contribution. First, we
recognize that the marginal distributions of the hourly zonal load variables
have heavy tails. So our first step is to check for the presence of heavy tail
in marginal hourly measurements, and when found, fit a Generalized Pareto
Distribution (GPD) to the marginals. We use them to process the data and give
them Gaussian marginals. We explain and illustrate the pros and cons of this
step with simulation illustrations on days for which the forecasts
overestimate the actual future values and days for which they underestimate
the actuals. In some respect, our approach is reminiscent of the Gaussian
copula approach as advocated for example in [14] and [18]. However, our
approach differs in significant ways from theirs. First because we fit heavy
tail distributions to the marginals, and then we use non-parametric estimation
techniques based on the LASSO $L^{1}$ regularization to estimate large
covariance matrix instead of the parametric forms used in [18]. Another major
feature of our model and its implementations is to be usable when the number
of historical observations used to fit the model is small compared to the
dimension of the state variables we model and simulate. This issue is often
encountered in modern high dimensional statistics. One of the major issues
with this feature of the data is that it implies a systematic singularity of
the empirical covariance matrix estimates. In order to overcome this
shortcoming, we use an approach based on graphical Gaussian models estimated
by a form of the graphical lasso algorithm glasso. In this way, we obtain
reasonable sparse proxies for the precision matrices, and work with Gaussian
models which can be simulated. As a further step to reduce the number of
degrees of freedom we use the gemini algorithm introduced in [22], to search
for precision matrices in the form of tensor products, allowing us to
disentangle the contributions of the spatial and temporal components to the
correlation structures.

The model introduced in this paper is built from the ground up. We first
introduce marginal models separately for the three features: load, wind power
and solar power. Historical data for load actuals and forecasts being only
available at the zone level, possible correlations between load and wind on
one hand, and between load and solar on the other, need to be detected and
implemented at the zonal level, even if eventually, wind and solar power need
to be modeled and simulated at the individual asset level.

The marginal model for load at the zone level is inspired by our earlier work
[15] on the ERCOT data publicly available from the Texas ISO web site. While
the load data used in the present paper is very different from the data used
in [15], load data is still at the zone level, and even though the zones are
different in their number and coverage, the modeling approach used in the
present paper is similar to the one of our previous work. The stochastic
models we developed for load, and wind and solar power are the fulcrum for a
Monte Carlo simulation engine providing scenarios on demand. Given a date and
a set of historical observations (including forecasts) prior to that date, the
simulation engine produces as many Monte Carlo scenarios as requested for the
next $24$ hourly values of load, wind and solar,

The simulation engine outputs scenarios in the format predicated by the format
of the historical data used to fit the models. A first version of the
simulation engine was implemented in a Python package by Michal Grzadkowski
and made available with extensive README files and examples which can be run
from Jupyter Notebooks. It can be found at [12].

###  1.2. Organization of the paper.

The paper is organized as follows: in Section 2 we describe the data used for
the analysis. The following Section 3 is devoted to the analysis of the
aggregate demand for electric power at the zone level. The fitted model is the
basis for the Monte Carlo simulation engine whose performance is demonstrated
by examples of batches of scenarios. Our asset level model is developed and
fitted to the wind power generation historical data in Section 4. It is shown
there that, when aggregated at the zone level, no significant dependence with
load is detected. The situation is very different for solar power. The first
part of Section 5 is devoted to the development of a marginal model for the
$24$ hours of solar power production for each of the $222$ solar stations for
which we have historical data. The second half of the section is devoted to
the aggregation of the solar production at the level of the load zones and to
the quantification of the dependencies between load and solar at these levels.
Load and solar generation are modeled together in Section 6. The main take-
away from this section is that, when designing Monte Carlo simulation engines
from the model, we first simulate load and solar power at the aggregate zone
level and then, from the zone solar production scenarios generated jointly
with load, we simulate $24$-dimensional solar generation for each asset,
conditioned by the value of the vector of aggregate production at the zone
level already simulated.

We conclude with a short recall of the main contributions of the paper in
Section 7.

##  2\. Problem Set-Up and Description of the Data

###  2.1. Case Study Data

For the purpose of illustration and demonstration, we use the first of the
data sets assembled by NREL for the Department of Energy program ARPA-E
PERFORM, [11]. The data covers the Texas region. This synthetic data set
contains time series for electricity demand at the regional level, and actuals
and forecasts of wind and solar power generation at the asset level. Several
time frequencies are provided, but for the sake of definiteness, we shall
limit ourselves to hourly data. The geographical resolution is the highest in
the sense that data is provided at the site-level.

The solar data comprises $226$ sites (including $22$ existing and $204$
proposed solar sites) in the Texas region. The solar power forecasts are
provided for $2018$ as probabilistic forecasts in the form of $1$ to $99$
percentiles. These probabilistic forecasts are generated using a Bayesian
Model Averaging (BMA) [11] procedure to the ensembles of $51$ deterministic
forecasts provided by the European Centre for Medium-Range Weather Forecasts
(ECMWF). However, since the ECMWF forecasts included in the NREL data set are
for $2017$ as well as $2018$, in order to maximize the size of historical data
to which we fit our model, we chose to ignore the NREL probabilistic forecasts
and use instead, a plain average of the $51$ ECMWF deterministic forecasts as
a point forecast. In summary, the day-ahead forecasts we use have a $12$-hour-
ahead lead time, a $24$-hour horizon, an hourly resolution, and a daily update
rate. For each day, at $12:00$ (issue time), we have a set of hourly forecasts
for the next day starting from $00:00$ ending at $23:00$ (forecast time). For
the actual solar power, the data is given with a $5$-min resolution for $2017$
and $2018$, so we computed an average for each hour to obtain hourly solar
power.

The wind data comprises $264$ sites (including $125$ existing and $139$
proposed sites) in the Texas region. The wind day-ahead forecasts are provided
by NREL in the same format as the solar data – probabilistic forecasts for
$2018$ and $51$ ECMWF forecasts for $2017$ and $2018$. However, the ECMWF
forecasts are highly biased and point forecasts obtained by plain averaging
yield much larger deviations. For this reason, we prefer to use the median of
the probabilistic forecasts as a point forecast. However, in order to
demonstrate the impact of the quality (or lack thereof) of the forecasts on
the performance of the output scenario generation, we also use the plain
average of the $51$ ECMWF forecasts. Note that for the median of the
probabilistic forecast, the available historical data to fit and test the wind
model is limited to $2018$ only, while for the average of the ECMWF forecast,
the available historical data covers $2017$ and $2018$. For the actual wind
power, the data is given with a $5$-min resolution for $2017$ and $2018$.
Again, we computed an average for each hour to obtain hourly wind power data.

The NREL data set also contains load data at the zone level. In fact, the load
actual and forecast data were collected from the ERCOT official website at
hourly resolution for $8$ load zones. The load forecasts were obtained from
ERCOT’s ”Seven-Day Load Forecast by Weather Zone”, see [5], by selecting the
data issue and forecast times to match those of the solar day-ahead forecasts.

######  Remark 1.

While the examples given in this paper are based on the NREL synthetic data
set described above, we also implemented our model on the IEEE’s Reliability
Test System (RTS) whose update is presented in [2]. The interested reader can
obtain simulation results for RTS by running the examples provided in the
Jupyter Notebooks which can be found on the GitHub repository of the
simulation platform PGscen, [12].

##  3\. Marginal Model for Electricity Loads

We first tackle the challenge of the design of the Monte Carlo scenario
simulation engine for the loads over the region.

###  3.1. Strategy

As explained in the introduction, we model the deviations (which are often
called errors in the literature on the subject) defined as

_load deviations_ = _load actuals_ – _load forecasts_.

Recall that the historical data we secured covers the period 2017-01-01 to
2018-12-31. For each day $d$ and for each load zone

$z\in\\{\text{Coast},\text{East},\text{Far West},\text{North},\text{North
Central},\text{South},\text{South Central},\text{West}\\},$

we read the vector of the $N_{lag}=24$ load point forecasts from the load
forecast data set, and for each lag $\ell\in\\{0,1,\cdots,N_{lag}-1\\}$, we
subtract the point forecast from the actual load found in the load data set.
This gives what we call the load deviation. The data set of all deviations so
constructed has the same structure as the data set of load forecasts.

1) For each load zone $z$ and for each lag $\ell\in\\{0,1,\cdots,N_{lag}-1\\}$
the load deviations form a daily time series of length $N$, the number of days
prior to the day we generate scenarios for. The first step of our analysis is
to analyze the marginal distributions of these $N^{L}_{zone}\times
N_{lag}=8\times 24=192$ stationary univariate time series.

2) For each load zone $z$ and for each lag $\ell$ we produce a Q-Q plot of the
empirical quantiles of the load deviations against the theoretical quantiles
of the standard Gaussian distribution $N(0,1)$. Figure 1 give a sample of $4$
of these Q-Q plots. They are typical of what we found throughout. This is
clear evidence that the marginal distributions of the load deviations have
_heavy tails_ , justifying our next step.

Figure 1. Normal Q-Q plots for load deviation remainders in the
$z=\text{Coast}$ region and lags $\ell=0,1,2$ and $\ell=3$ hours ahead.

3) For each load zone $z$ and for each lag $\ell$ we fit a Generalized Pareto
Distribution (GPD), see for example [4, Chapter 2], to the load deviations. We
denote such a distribution as $G^{L}_{z,\ell}$ for the sake of later
reference. Throughout this work, in order to manipulate GPD distributions, we
use functions of the R library Rsafd [3] developed for the analysis of the
examples and problem sets of the book [4], which we wrapped into a Python
package.

4) For each load zone $z$ and for each lag $\ell$ we transform the load
deviation time series into a uniform time series by computing the cumulative
distribution function of the GPD just fitted on each of the entries of the
time series. If the fitted GPD is a good fit, we should expect the marginal
distribution of the new time series to be uniform. The histogram on the left
pane of Figure 2 is a testament to this claim.

Figure 2. Transformations of the load deviations in the Coast region $4$ hours
ahead.

5) For each load zone $z$ and for each lag $\ell$ we transform the uniform
time series just obtained into a time series with a standard Gaussian
distribution $N(0,1)$. In order to do that, we just compute the quantile
function of the standard Gaussian distribution on each of the entries of the
uniform time series. The normal Q-Q plot on the right pane of Figure 2 gives a
graphical evidence of the fact that the marginal distribution of the resulting
time series is indeed a standard Gaussian distribution.

6) At this stage of the analysis, for each day $d$, we have a
$N^{L}_{zone}\times N_{lag}=8\times 24=192$ dimensional numerical vector, each
component of this vector having a standard Gaussian distribution. While we do
not know if this $192$-dimensional vector is actually _jointly Gaussian_ , we
shall behave as if that was. Given the high dimensionality of the vectors when
compared to the small amount of historical data available to estimate the
covariance structure of such a Gaussian vector, we expect that any empirical
estimate of the covariance matrix will be degenerate and hence not invertible.
So in order to estimate the precision matrix capturing the dependencies
between the various components of the vector, we use a regularization
procedure and we proceed to fit a _Gaussian Graphical Models_ using LASSO. The
interested reader will find a readable introduction to Gaussian graphical
models, including LASSO estimation in [20, Chapter 11]. Here we could use the
R package glasso to fit the graphical model to our modified load data.
However, we believe that we would miss the special structure of our $192$
dimensional vector which has a $8$-dimensional spatial component and a
$24$-dimensional temporal component. Taking one more step to reduce the
mismatch between the high dimension of the vector and the short historical
data, and relying on the intuition developed for classical statistical
analysis of variance with multiple factors, we would like to think of the
$192\times 192$ covariance and precision matrices as the tensor product (also
known as Kronecker products) of a $8\times 8$ covariance and precision
matrices capturing the dependencies between the spatial components, and a
$24\times 24$ covariance and precision matrices for the temporal component.
Searching for a precision matrix and a covariance matrix in this form
dramatically reduces the number of parameters to be estimated. In the
statistical literature, these special covariance matrices are called
_separable_. Testing for the presence of such a special structure is sometime
possible. See for example [9, 8] and the more recent paper [6], even though
none of these tests have enough power in our case because of the lack of
historical data. In any case, we shall work under this assumption here, and we
shall use the variation on the glasso algorithm called GEMINI which was
proposed in [22] to take advantage of the separable structure.

Figure 3 shows the graph of the conditional temporal dependence in the left
pane, as well as the spatial conditional dependence component provided by the
GEMINI algorithm.

Figure 3. Temporal component of the conditional correlation (left), spatial
component of the conditional correlation (middle), geographical map of the 8
load zones in ERCOT.

The graph of the temporal dependencies comprises two types of links. The thick
links form a perfectly ordered linear chain. It points clearly at some form of
Markovian structure. However, the thinner links show a residual dependence for
hours two lags apart, which is quite intuitive. The edges, as well as their
strengths in the graph of the spatial dependencies form a pattern consistent
with the contiguity of the load regions as shown in the right pane of Figure
3.

7) The next step of the analysis is to generation of Monte Carlo samples from
the model. Given the covariance matrix estimated above, we can easily generate
as many Monte Carlo samples from the $192$-dimensional mean zero Gaussian
distribution with this covariance matrix. Computing the cumulative
distribution function of the standard Gaussian distribution on each of the
entries provides sample of $192$-dimensional random vectors with uniform
marginals. Each of the $192$ components corresponds to a specific zone $z$ and
a specific lag $\ell$ so if we now compute the quantile function of the GPD
distribution $G^{L}_{z,\ell}$ fitted to the original load deviations for zone
$z$ and lag $\ell$, we obtain samples of a $192$ dimensional random vector
with marginal distributions originally fitted to the load deviation time
series. At this stage, it is important to emphasize that these $192$
components depend upon each other in the way dictated by the covariance
structure found in the estimate of the graphical Gaussian model. The astute
reader will have certainly noticed that we are actually using a
$192$-dimensional copula to build the dependencies between the
$G^{L}_{z,\ell}$ marginals of this $192$-dimensional vector.

At this stage, it is very easy to complete the simulation cycle by adding back
the forecasts, in order to end up with Monte Carlo scenarios of the actual
loads.

Figure 4 shows what the scenarios look like vis-a-vis the actual and forecasts
loads for May 21, 2018.

Figure 4. $1000$ Monte Carlo scenarios generated for the next $24$ hours for
May 21, 2018. The red line gives the actual load demand observed over these
$24$ hours while the blue line shows the point forecasts. The dashed lines are
$10$ randomly selected scenarios. The light blue bands are the traces left by
the Monte Carlo scenarios after we remove the smaller $1\%$ and the larger
$1\%$ of the bunch.

######  Remark 2.

Notice that we assumed some form of stationarity of the load time series as we
worked under the assumption that the load deviations have a single maarginal
distribution, irrespective of the day of the year (whether $d$ is in the
summer or the winter, or $\ldots$) or the lag $\ell$. Our reason for doing so
is that we expect that the major trends and seasonal effects are captured by
the exogenously provided forecasts which are usually constructed from a large
number of meteorological variables and Bayesian averaging of many forecasting
models.

###  3.2. The Impact of Fitting Heavy Tail Distributions

Here, we discuss the impact of our decision to check for the presence of heavy
tails and to possibly fit generalized Pareto distributions to the marginal
laws of the deviations. In order to do so, we provide a detailed examination
of the smoothness of the Monte Carlo scenarios and the coverage provided by
these scenarios.

Repeated overestimation by the forecasts create instances of negative values
for the deviations, and if these instances are significant because of their
sizes, the empirical distribution of the deviations can exhibit a heavy lower
tail which cannot be accounted for by Gaussian distributions. We illustrate
such an instance in Figure 5 with the example of the load zone Coast on May
21, 2018. The left pane of the figure shows that the actual loads (red curve)
hover around the lower boundary of the band created by the Monte Carlo
scenarios, especially during the later hours. Recall that these scenarios are
obtained by adding sample scenarios of the deviations to the forecasts (blue
curve). Fitting a distribution with a heavy lower tail guarantees that the
corresponding scenarios will cover the actual values. The latter should be
viewed as extremes, possibly rare events, but still in the realm of the model.
The right pane shows the same plot produced from Monte Carlo scenarios
obtained from a Gaussian model ignoring the presence of heavy tails. While one
could think that having a narrower coverage band is desirable, the fact that
the actual values are so far from what the model considers to be reasonable is
a major flaw of this form of model.

Figure 5. Actual loads (red curves) and forecasts (blue curve) as well as a
gray band covering $10,000$ Monte Carlo scenarios, as well as the sample
scenario (dashed lines) with the minimum mean absolute percentage error (MAPE)
for the $24$ hours following May 21, 2018. The left pane was produced using
the model described in the paper, while the right pane was produced from a
Gaussian model ignoring the possible presence of heavy tails in the
distributions of the deviations.

In a similar fashion, we can discuss the case of underestimation by the
forecasts. Clearly such instances can have dire consequences for the
management of the electric grid, and they are feared by system operators.
Figure 6 provides an illustration which makes a resounding case for the use of
heavy tail distributions in the model. We use the example of the load in the
West region on February 11, 2018. With the same convention as before, the left
pane of Figure 6 shows that the actual loads for the next $24$ hours (red
curve), while much higher than the forecasts (blue curve) are mostly covered
by the scenarios. So even on the tail end of the possibilities given by the
Monte Carlo simulations, it can still be expected to occur with significant
probability and the operator of the grid should not dismiss such a pattern. In
contrast, the right pane of Figure 6 shows the result of the same Monte Carlo
simulations when the model does not attempt to detect and fit heavy tail
distributions. Clearly the curve of actual loads is far beyond the simulation
band, and an operator relying on such a model will consider the $24$ actual
values as a possibility with zero likelihood.

Figure 6. Same plots as in Figure 5 for February 28, 2018 to illustrate the
possible consequences of an underestimation by the forecasts.

The two instances discussed above clearly illustrate the pros of checking for
the presence of heavy tails in the marginal distributions of the deviations,
and when detected, of using generalized Pareto distributions to process the
data.

##  4\. Marginal Model for Wind Power

This section will be shorter than the previous one. Indeed, several steps of
the analysis are the same as for the loads, so we avoid most repetitions by
emphasizing only the significant differences when they are relevant.

###  4.1. Modeling Steps

As before, we introduce the _deviations_ :

_wind power deviations_ = _wind power actuals_ – _wind power forecasts_ ,

and prepare the data set of wind power deviations exactly as for the loads,
the only noticeable difference being that we have now $N_{w}=264$ production
assets instead of $N_{z}=8$ zones. Obviously, we still have the same
$N_{lag}=24$ lags.

As for the step 2) above, the marginal distributions do not exhibit heavy
tails. Indeed, the power produced by a single wind station is a non-negative
number bounded from above by the capacity of the wind farm. So for a specific
wind farm, the deviation should have a distribution with support contained in
an interval $[-c,+c]$ where $c$ is the capacity of the generation asset. So
for each of the $6336$ $(w,\ell)$ asset-lag couples, we use the empirical
cumulative distribution function (cdf) to fit the marginal distribution of the
_stationary_ time series, and we apply this cdf to each of the entries to
produce a series with uniform marginals, and eventually standard Gaussian
after applying the quantile function of the Gaussian distribution. This takes
care of the steps 3) - 5) of the strategy used for the load data. Fitting a
Gaussian graphical model like we did in step 6) above with the GEMINI
algorithm to capture the precision and covariance matrices leads to the analog
of Figure 3 which we reproduce as Figure 7. In the estimation of the spatial
components of the covariance matrix using LASSO, we chose the regularization
parameter to be a $N_{w}$-by-$N_{w}$ symmetric matrix so that the $ij$-th
element is proportional to the Euclidean distance between the $i$-th and
$j$-th wind farms.

Figure 7. Temporal component of the conditional correlation (left), spatial
component of the conditional correlation (right) using the plain average of
$51$ wind ECMWF forecast.

###  4.2. Impact of the Quality of the Forecasts.

As mentioned in Section 2, we implemented simulation models for wind using two
different sets of forecasts: the median of the probability forecasts and a
plain average of the $51$ ECMWF forecasts. We show the temporal and spatial
components of the conditional correlations using the two forecasts in Figures
7 and 8. For both cases, the spatial dependency graphs are consistent with the
geographical locations of the wind power production assets. Figure 7 shows
that using point forecasts obtained from the ECMWF forecasts, the graphical
model shows the usual _Markovian-like_ temporal dependence structure given by
the clear chain pattern of the temporal dependency graph. However, the
temporal dependency graph derived from the probabilistic forecasts doesn’t
indicate a clear pattern. Moreover, it depends strongly on the choice of the
LASSO regularization parameters. For example, the graph shown in Figure 8
exhibits many isolated nodes. Moreover, although the temporal dependency graph
contains a chain pattern for the hours from $8$ to $14$, the largest
correlation coefficient in the graphical model is less than $20\%$.

Figure 8. Temporal component of the conditional correlation (left), spatial
component of the conditional correlation (right) using the median of the
probabilistic forecast.

###  4.3. Monte Carlo Generation Specifics

Exactly as in the case of the load simulations, we can generate Monte Carlo
scenarios for the Gaussian vectors whose distributions were identified from
the graphical Gaussian model, and remove the Gaussian marginals by computing
the Gaussian cumulative distribution function on the entries of these vectors.
This way of building dependencies in the Monte Carlo simulation is reminiscent
of the use of the Gaussian copula. However, transitioning to the actual power
productions (recall step 7 of our discussion of the load) requires a modicum
of care. Indeed, in the case of wind power, the marginal distribution of the
deviations strongly depends upon the actual value of the forecast. See Figure
9 for empirical evidence.

Figure 9. Scatterplots of the deviations (vertical axis) against the
corresponding forecasts (horizontal axis).

So in order to account for this phenomenon, we restore the marginal
distributions of the scenarios by using the quantile function of conditional
distribution of the deviations given a forecast value. We estimate this
conditional distribution by aggregating all the couples (forecast, deviation)
in the data, grouping the values of the forecasts in a finite set of bins, and
estimating the distribution of the actual generation above each bin. Once this
is done, for each couple (asset, lag), we generate the value of each scenario
by adding to the forecast, the value of the conditional quantile function for
the bin in which the forecast landed on, computed on the Monte Carlo value
obtained by the Gaussian simulation made uniform.

For illustration purposes, we show the impact at the asset-level of the
quality of the forecasts by the plots of wind power scenarios for $8$ wind
farms randomly selected from the $264$ wind farms. The Monte Carlo scenarios
are generated from our model using the mean of the ECMWF forecasts in Figure
10, and the median of the probabilistic forecasts in Figure 11. Because the
probabilistic forecasts were constructed from a bias correction procedure
applied to the ECMWF point forecasts, the forecasts seen in Figure 11 are much
closer to the actuals than in the case of Figure 10. The precision with which
these forecasts predict the actual values of the wind generation is
remarkable. This being also true for the historical data used to fit the
models, the bands covered by the scenarios in Figure 11 are narrower.

Figure 10. $1000$ Monte Carlo scenarios generated for the next 24 hours for
May 21, 2018, using the mean of the ECMWF forecasts as point forecasts. The
red line gives the actual wind power observed over these $24$ hours, the blue
line shows the point forecasts for these $24$ hours, while the dashed lines
are 10 randomly selected scenarios. The light blue bands are the traces left
by the Monte Carlo scenarios after we remove the smaller $1\%$ and the larger
$1\%$ of the bunch.

Figure 11. Same caption as for Figure 10 except for the fact that the
scenarios are now generated from a model using the median of the probabilistic
forecasts as point forecasts.

##  5\. Marginal Model for Solar Power

The fact that solar power production follows a very different pattern than
wind power production, is in no small part due to the stark differences
between nocturnal and diurnal productions, so we adopt a different strategy to
fit a model to the data.

In line with our analysis of load demand and wind power production, we work
with the deviations between the solar power production actuals and forecasts,
i.e.,

_solar power deviations_ = _solar power actuals_ – _solar power forecasts_.

The dataset consists $N_{s}=226$ production assets and $N_{lag}=24$-hour time
horizon for the period 2017-01-01 to 2018-12-31. For each of the $5424$
$(s,l)$ asset-lag couples, the data is in the form of a time series of length
$730$.

In order to make sure that our model (and the subsequent Monte Carlo
simulations) are non-anticipative (i.e. do not use a crystal ball), for each
given day $t$ of 2018, we collect the historical data comprising the $n$ days
preceding $t$, and the days of 2017 which fall in the $2n+1$ days window
centered around day $t$ in 2017. We use this approach for any day $t$ in 2018,
and rely on the historical data identified in this way to fit the model from
which we generate Monte Carlo scenarios for day $t$. In typical experiments we
use the default value $n=50$. This is reasonable, even when $t$ is early in
the year. For instance, if we are interested in creating scenarios for
2018-07-01, the historical data points to fit the model would be 2017-05-12 to
2017-08-20 and from 2018-05-12 to 2018-06-30.

To transform the data so that the marginal distributions become Gaussian, we
first use the empirical cumulative distribution function to produce time
series with uniform marginals, and then the standard Gaussian quantile
function to end up with Gaussian marginal distributions.

Next, we perform a Principal Component Analysis (PCA) of the transformed data.
Denoting by $N$ the number of days used to fit the model, we aggregate the
transformed data over all historical days and the assets and end up with a
$226N\times 24$ matrix, recall that we have $226$ solar stations, one column
for each hour of the day, and one row for each couple (asset, day). We then
compute the PCA orthonormal basis of $24$ loadings, and a $226N\times 24$
matrix of components of the original data onto the loadings basis. Given the
proportions of the variation in the data explained by the successive loadings,
we chose to keep, for each of the $226N$ couples (asset, day), only the first
$k$ entries. The switch from the original data to their first $k$ components
is intended to resolve the issue of the diurnal zero production, more than for
its _data compression_ effect. We now reshape the data into $k\times 24$ time
series of length $N$, and standardize the data to have zero mean and unit
variance. We then fit a graphical Gaussian model using GEMINI to obtain a
$226$-by-$226$ and a $k$-by-$k$ covariance matrices that capture the
dependencies between the spatial components and dependencies between the PCA
loadings. To generate scenarios, we first draw Monte Carlo samples from the
mean zero Gaussian distribution with covariance matrix defined as the tensor
product of the two covariance matrices computed from the graphical model. We
then reconstruct $24$-dimensional vectors using the orthonormal basis of the
PCA loadings. We complete the simulation by adding back the forecasts. In
Figure 13, we show the graph of the spatial conditional temporal dependence,
which is again consistent with the geographical locations of the solar power
production assets.

Figure 12 shows an example of asset-level solar production scenarios for $8$
solar farms randomly selected from the $226$ solar farms. As we can see, the
values of the solar production given by the scenarios are only non-zero during
the day. This is the result of our use of the projection of the data on the
first $k$ loadings of the PCA.

Figure 12. $1000$ Monte Carlo scenarios generated for the next 24 hours for
May 21, 2018, using the plain average of the ECMWF forecasts as point
forecasts. The red line gives the actual solar power observed over these $24$
hours, the blue line shows the point forecasts for these $24$ hours, while the
dashed lines are $10$ randomly selected scenarios. The light blue bands are
the traces left by the Monte Carlo scenarios after we remove the smaller $1\%$
and the larger $1\%$ of the bunch.

Figure 13. Spatial component of the conditional correlation for the analysis
of the solar deviations using PCA.

##  6\. Joint Model and Consolidated Simulation Engine

Once the three marginal models discussed earlier are fitted, and the
corresponding simulation engines have been developed, a joint model and a
joint simulation engine can be obtained following the steps described in this
section.

###  6.1. In-sample vs Out-of-sample Modeling and Simulation

To guarantee that our model (and the subsequent Monte Carlo simulations) are
non-anticipative, for each given day $t$ of 2018, we fit our model to the
bunch of historical data identified earlier in the case of solar power
generation. However, if we want to use the wind point forecasts given by the
median of the probabilistic forecasts provided by NREL, we do not have
historical data for 2017. So when the day $t$ is early in the year, we fit the
model with data for days in a window covering $t$, extending in the future by
$n$ days and excluding $t$. So, for those values of $t$, the estimation of the
model is _in-sample_. This is unfortunate as it should not have to be if we
had enough historical data.

###  6.2. Quantifying and Implementing the Dependencies

Load data being available at the zone level only, the strategy for the joint
modeling of load, wind power and solar power is based on the following steps.

  1. (1)

The first step of the analysis is again a form of _Gaussianizeation_ of the
historical data of load, wind and solar power productions. We do that
separately, just as we did for the marginal models;

  2. (2)

For each of the load zones and for each lag $\ell$, we aggregate the wind and
solar Gaussianized power production historical data from the individual asset
levels to the zone levels. So for each of the $8$ load zones and for each lag
$\ell$, we sum up all the Gaussianized time series of the individual wind
asset deviations to create one single wind power deviation series for this
zone. Similarly, for each of the $8$ load zones and for each lag $\ell$, we
sum up all the Gaussianized time series of the individual solar asset
deviations to create one single solar power deviation series for this zone;

  3. (3)

Next, for each of the load zones, we aggregate the zonal load, wind and solar
Gaussianized power productions across time. As solar plants only generate
power during the day, we only aggregate across those hours. In particular, if
we denote the lags during sunrise and sunset as $\ell_{sunrise}$ and
$\ell_{sunset}$ respectively, then we sum up all the Gaussianzed zonal load
time series to create one single load series for that zone during that day.
Similarly, we aggregate wind (and solar) to create one single series for that
zone during that day.

  4. (4)

For each load zone $z$, we concatenate the three time series into one three
dimensional time series and fit a joint Gaussian graphical model at the zone
level, and using a plain glasso algorithm. For this step, we do not use the
gemini variation on the glasso algorithm because the temporal dependency
graphs of load, wind and solar are sufficiently different for us to be
convince that the assumption of _joint separability_ of the covariance matrix
is unreasonable.

  5. (5)

When running the glasso algorithm with load. solar power and wind power, we
find that the nodes of the dependency graph corresponding to wind are
disconnected from the nodes corresponding to load and solar. Since we are
still in the _Gaussian world_ when we perform this step of the analysis, our
conclusion is that, at least at the zonal level, wind power generation is
statistically independent of the electricity demand and the solar power
generation. As a.result, wind is modelled separately and Monte Carlo scenarios
are generated for wind power from the marginal model described in Section 4

  6. (6)

As for diurnal load and solar power, Figure 14 provides the spatial dependency
graph as given by the precision matrix, while Figure 15 provides a heat-map
for the $16\times 16$ correlation matrix of the diurnal aggregates of load and
solar power over the $8$ load zones. The dependence highlighted in Figure 14
by the links between the load zones is consistent with the geographical
locations of these zone, and the fact that the node corresponding to load in
the ”Far-West” zone is isolated is not a surprise. Links between different
solar zones are also natural. Strong links between load and solar power for
the same zone are also significant and expected, for example load and solar in
zone South Central, or zones East and Coast, and to a lesser degree zone
South.

Figure 14. Graphical LASSO estimate of the conditional correlation of load and
solar power in the joint model described in the text.

As expected, the corresponding correlation matrix whose heat-map is given in
Figure 15 is not sparse. The bottom-right fourth of the matrix confirms the
strong correlations between the zone level diurnal solar power productions.
These correlations are weaker among the zonal loads. The green blocks
appearing in the bottom-left quadrant of the matrix confirm the remarks we
made above concerning the correlations between load and power production in
some regions.

Figure 15. The correlation heatmap of the load-solar joint model.

###  6.3. Putting Everything Together with Conditional Simulation

Giving the above discussion, in order to generate Monte Carlo scenarios for
the entire system on a given day $t$, we can first simulate samples of $24$
hours wind power generation for the $264$ wind farms using the marginal model
introduced and fitted in Section 4. Next, and independently of the scenarios
generated for wind power, we generate jointly correlated Monte Caro scenarios
for load and solar power according to the following steps based on conditional
simulation of Gaussian processes.

  1. (1)

We generate jointly Gaussian scenarios for load and zonal solar power as their
daily sum of the lags from $\ell_{sunrise}$ to $\ell_{sunset}$. We do just
that using the joint Gaussian model estimated from the glasso algorithm
implemented as described above.

  2. (2)

Next for each such scenario, we generate a $8\times 24$ ($8$ for each load
zone and $24$ for the $24$ time lags) Gaussian load Monte Carlo scenario from
the marginal distribution fitted using the gemini algorithm as described in
Section 3, conditional on the fact that the sum of the lags in the range
$[\ell_{sunrise},\ell_{sunset}]$ are exactly equal to the values of the load
scenario generated in step 1. We do that for the eight load zones
_simultaneously_ to make sure that we do not wipe out the dependencies between
the zones!

  3. (3)

Next for each scenario generated in step 1, we generate a $226\times 24$
($226$ for each of the solar assets, and $24$ for the $24$ time lags) Gaussian
load Monte Carlo scenario from the marginal distribution fitted using the PCA
based algorithm described in Section 5, conditional on the fact that the zonal
aggregates of the lags in the range $[\ell_{sunrise},\ell_{sunset}]$ are
exactly equal to the values of the particular solar scenario generated in step
1. We do that for the $226$ solar assets _simultaneously_ to make sure that we
do not wipe out the dependencies between the individual asset!

  4. (4)

Once we have Monte Carlo scenarios for the Gaussian versions of zonal loads
and individual solar assets, we can undo the Gaussian nature of the scenarios
as explained in Section 3 and Section 5 and reinstate the non-necessarily
Gaussian marginal distributions.

##  7\. Conclusion

In this paper, we propose an original model for electricity demand and wind
and solar power generation together with an implementation producing joint
Monte Carlo simulations for the purpose of unit commitment and economic
dispatch performed by highly sophisticated optimization programs. For this,
while electricity load can remain at the zonal level, wind and solar
productions need to be modeled and simulated at the individual asset level.
The main thrust of our model is to capture the possible and intricate
dependencies both temporal and geographical. We consider hourly data, but our
model and our simulation algorithms are agnostic to the time frequency and can
be ported to higher frequencies as long as the historical data is available.

The main novelty of our approach resides on 1) checking for the presence of
heavy tails in the marginal distributions and the use of generalized Pareto
distributions to _Gaussianize_ the data in the spirit of the use of Gaussian
copulas in portfolio theory; 2) using graphical Gaussian model and producing
precision matrix proxies with the gemini variation on the glasso algorithm
when appropriate, disentangling the temporal and spatial contributions to the
correlations, and allowing for implementations with high dimensional state
variables. We illustrate the potential of the model and its implementation on
synthetic data produced by NREL. Our implementation, called PGscen, is part of
the GitHub

https://github.com/PrincetonUniversity/PGscen

##  8\. Acknowledgments

Both authors were partially supported by ARPA-E grants DE-AR0001289 and DE-
AR0001390 under the PERFORM program of the US Department of Energy. A
preliminary version of the paper was presented in November 2020 during the
INFORMS Annual Meeting, and on December 21, 2020 at a DOE Review Meeting. We
would like to thank Mike Ludkovski (University of California at Santa
Barbara), Ronnie Sircar (Princeton University) and Glen Swindle (Scoville Risk
Partners) for enlightening conversations on the content of the paper, and
Michal Grzadkowski for the development of the GitHub package PGscen
implementing the model presented in the paper.

## References

  * [1] A.T. Al-Awami and M.A.El-Sharkawi.  Statistical characterization of wind power output for a given wind power forecast.  In North American Power Symposium (NAPS), pages 1–4, 2009.
  * [2] A. Ehlen J.Ikäheimo J. Jorgensen D. Krishnamurphy J. Lau B. McBennett M. O’Connell E. Preston A. Staid G. Stephen J.P. Watson C. Barrows, A. Bloom.  The ieee reliability test system: A proposed 2019 update.  IEEE Transactions on Power Systems, 35(1):119 – 127, January 2020\.
  * [3] R. Carmona.  Rsafd.  https://rdrr.io/rforge/Rsafd/, 2014.
  * [4] R. Carmona.  Statistical Analysis of Financial Data in R, Second Edition.  Springer Texts in Statistics. Springer Verlag, 2014.
  * [5] ERCOT.  Ercot load.  GitHib repository, 2021.
  * [6] P. Guggenberger, F. Kleibergen, and S. Mavroeidis.  A test for Kronnecker product structure covariance matrix.  Journal of Econometrics, available on line, 2022.
  * [7] T. Hong and S. Fan.  Probabilistic electric load forecasting: A tutorial review.  International Journal of Forecasting, 32:914 – 938, 2016.
  * [8] B. Li, M. Genton, and M. Sherman.  Testing the covariance structure of multivariate random fields.  Biometrika, 95:813 – 829, 2008.
  * [9] M. Mitchell, M. Genton, and M. Gumpertz.  A likelihood ration test for separability of covariances.  Journal of Multivariate Analysis, 97:1025 – 1043, 2006.
  * [10] J.M. Morales, R. Minguez, and A.J. Conejo.  A methodology to generate statistically dependent wind speed scenarios.  Applied Energy, 87:843–855, 2010.
  * [11] NREL.  Arpa-e perform ercot dataset produced by nrel.  GitHib repository, 2021.
  * [12] ORFEUS.  Pgscen, a python package for load, wind and solar power simulation.  GitHib repository, 2021.
  * [13] P. Pinson and R. Girard.  Evaluating the quality of scenarios of short-term wind power generation.  Applied Energy, 96:12–20, 2012.
  * [14] P. Pinson, H. Madsen, H. Nielsen, G. Papaefthymiou, and B. Klöckl.  From probabilistic forecasts to statistical scenarios of short-term wind power production.  Wind Energy, 12:51–62, 2009.
  * [15] X. Yang R. Carmona.  Glasso model for electric load and wind power and monte carlo scenario generation.  Technical report, Princeton University, 2021.
  * [16] A. Staid, J.P. Watson, R.J.B. Wets, and D.L. Woodruff.  Generating short-term probabilistic wind power scenarios via nonparametric forecast error density estimators.  Wind Energy, 20:1911–1925, 2017.
  * [17] J. Tastu, P. Pinson, E. Kotwa, H. Madsen, and H.A. Nielsen.  Spatio-temporal analysis and modeling of short-term wind power forecast errors.  Wind Energy, 14:43–60, 2011.
  * [18] J. Tastu, P. Pinson, and H. Madsen.  Space-time trajectories of wind power generation: Parameterized precision matrices under a gaussian copula.  In Modeling and Stochastic Learning for Forecasting in High Dimensions, volume 217 of Lecture Notes in Statistics, pages 267 – 296\. Springer Verlag, 2015.
  * [19] D.W. van der Meer, J. Walfen, and J. Munkhammar.  Review on probabilistic forecasting of photovoltaic power production and electricity consumption.  Renewable & Sustainable Energy Review, 81:1484 – 1512, 2018.
  * [20] M.J. Wainwright.  High-Dimensional Statistics: A Non-asymptotic Viewpoint.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2019.
  * [21] Y. Zhang, J. Wang, and X. Wang.  Review on probabilistic forecasting of wind power generation.  Renewable & Sustainable Energy Review, 32:255 – 270, 2014.
  * [22] S. Zhou.  GEMINI: Graph estimation with matrix variate normal instances.  Annals of Statistics, 42:532–562, 2014.