Large-scale Probabilistic Forecasting in Energy Systems using Sparse
Gaussian Conditional Random Fields
Matt Wytock and J. Zico Kolter
Abstract— Short-term forecasting is a ubiquitous practice ina wide range of energy systems, including forecasting demand,renewable generation, and electricity pricing. Although it is
known that probabilistic forecasts (which give a distribution
over possible future outcomes) can improve planning andcontrol, many forecasting systems in practice are just used as“point forecast” tools, as it is challenging to represent highdimensional non-Gaussian distributions over multiple spatial
and temporal points. In this paper, we apply a recentlyproposed algorithm for modeling high-dimensional conditional
Gaussian distributions to forecasting wind power and extend
it to the non-Gaussian case using the copula transform. On awind power forecasting task, we show that this probabilistic
model greatly outperforms other methods on the task of
accurately modeling potential distributions of power (as would
be necessary in a stochastic dispatch problem, for example).

I. I NTRODUCTION
Forecasting, the task of predicting future time series
from past observations, is ubiquitous in energy systems. As
well-known examples, electricity system operators routinely
forecast upcoming electrical load and use these forecasts inmarket planning [17], [19]; wind farms forecast future power
production when offering bids into these markets [7], [11],[15]; and there is a growing use of forecasting at the microscale for coordinating smart grid operations [1]. Despite theirubiquity and the complexity of many forecasting methods,
most methods are ultimately employed as “point forecast”
strategies; users train a system to output point predictions ofupcoming values, typically to minimize a metric such as root
mean squared error. However, for many complex control andplanning tasks, such point forecasts are severely limited: theprocesses that make up electrical demand, wind power, etc,
are stochastic systems and the notion of a “perfect forecast”
is unattainable. Thus, probabilistic forecasts, which output
a distribution over potential future outcomes instead of asingle prediction, are of substantial practical interest. Indeed,studies have demonstrated that in the context of electricaldemand and wind power, probabilistic forecasts can offer
substantial benefits over point predictions [14].
Unfortunately, the challenge of probabilistic forecasts isthat it is often very hard to describe the joint distribution
over all predicted values because many variables of interest
are highly non-Gaussian and it can be difficult to accurately
model correlations in a high-dimensional output space. For
this reason, most of the literature on probabilistic forecasting has often made simplifying assumptions, for example
M. Wytock and J. Z. Kolter are with the School of Computer
Science, Carnegie Mellon University, Pittsburgh, Pennsylvania, USA{mwytock,zkolter}@cs.cmu.edu

using specific forms of Gaussian linear models, such as
autoregressive moving average (ARMA) models (e.g. [10]);or only predicting non-Gaussian marginal distributions for
single output variables (e.g. [9]). Indeed, past work has
explicitly highlighted the challenge of developing models
that can capture joint distributions over future values.
In this paper, we apply and extend a recently-proposed
method for learning high-dimensional conditional Gaussian
distributions, called sparse Gaussian conditional random
fields, resulting in a method that is able to both capture
high-dimensional correlations between a very large number of output variables and model non-Gaussian marginal
probabilities. We accomplish this by 1) applying recent
techniques in machine learning and statistics, which model
high-dimensional correlations by exploiting sparsity in the
inverse covariance matrix, and by 2) using copula transforms
to capture non-Gaussian marginal distributions. We apply
this model to the task of wind power forecasting (a setting
where capturing non-Gaussian distribution and temporal andspatial correlation is critical), and show that it is able to makeprobabilistic predictions far better than forecasts that do not
explicitly model spatial and temporal correlation. We also
highlight the applicability of the method to challenging open
problems in this area such as wind power ramp prediction.
II. T HE PROBABILISTIC FORECASTING SETTING
We consider the following setting: let zt ∈ Rn denotea vector-valued observation at time t; for example, the ithelement of zt , denoted (zt )i could denote the power output
by a particular wind farm at time t, and i could range overa collection of wind farms. We let wt ∈ Rm denote a set of(known) exogenous variables that may affect the evolution of
the sequence; for example, wt may include the current timeof day, day of the year, and even external variables suchas wind forecasts at upcoming time points. The goal of ourforecasting setting is to predict Hf future values given Hp
past values and the exogenous variables:
Given zt−Hp +1:t , wt , predict zt+1:t+Hf .

(1)

This setting can be referred to as the vector autoregressive
exogenous (VARX) setting [13]; however, this terminology
is also used to describe a particular form of probabilistic
model for the sequence, so we just refer to it generally as amultivariate forecasting problem.
Although predicting a single estimate of future observations from past observations can be very useful in many
situations, we often want to understand more broadly thedistribution of possible future observations given past observations and exogenous factors, denoted asp(zt+1:t+Hf |zt+1:t+Hf , wt ).

(2)

We are particularly interested in the case where individual
observations zt are high-dimensional and we want to predict
their evolution over a relatively long time horizon, resulting
in a high-dimensional probability distribution. The task ismade more challenging by the fact that the observations may
not have Gaussian distributions (indeed, in the case of ourwind power setting, they typically do not) and by the fact that
we may have relatively few past observations upon which tobuild our high-dimensional model.
A. Relation to existing settings and models
A common method for handling such settings is a vector
autoregressive exogenous model in which we model future
observations as a linear combination of past observations,exogenous variables, and a Gaussian noise term
Hp −1

zt+1 =

X

Θi zt−i + Ψwt + t

(3)

i=0

where Θi ∈ Rn×n and Ψ ∈ Rn×m are model parameters,and t ∼ N (0, Σ) is a zero-mean Gaussian random variablewith covariance Σ. If we want to make forecasts over a
multiple future time points, we can iteratively apply thismodel or explicitly build an additional VARX model topredict zt+2 from past values.
A common extension to the autoregressive framework is
to add a moving sum of noise variables, resulting in the
autoregressive moving average (ARMA) model (these andsimilar extensions are described in several sources, e.g. [3]).
Additionally, we may choose Θ such that the overall system
has at least one unit root (the ARIMA model) or consider
only a certain periodic set of past observations (the seasonal
ARIMA model). These and other approaches have been used
extensively in the literature, and while they do impose a joint
probabilistic model over future observations, they are very
limited in the form of this distribution (multivariate Gaussianwith a particular covariance matrix).
In the sequel we will consider forecasting using a modelthat allows for more general dependencies between predictedvariables, can capture non-Gaussian marginal distributionof the variables via a copula transform, and which can belearned from very little data by exploiting sparsity. We focus
largely on extensions of the pure autoregressive setting, but
the model can also be extended to the ARMA setting by
introducing additional latent variables.
III. F ORECASTING WITH THE SPARSE G AUSSIAN
CONDITIONAL RANDOM FIELD

Here we describe the sparse Gaussian conditional random
field (SGCRF), a model that has has recently been proposedby several authors [16], [21], including our own work [20]which develops the algorithm we apply here. In this paper,
we focus on the extension of this method to the case of nonGaussian output variables and the application to probabilistic

forecasting in energy systems. For a full description of thealgorithm and theoretical analysis, see [20] which focuses on
the pure machine learning aspects of the model.
For simplicity of notation, we refer to the set of all known
variables as a single vector x ∈ Rn , while the unknown
variables we are attempting to predict are given by y ∈ Rp




zt
zt+1
 zt−1  zt+2 






..
,
y
=
(4)
x=


.
..
.




.
 zt−Hp +1 
zt+Hfwt
In this section we first present the pure SGCRF formulation,
which models y|x as a Gaussian distribution, exploiting
sparsity to capture high-dimensional correlation using very
few model parameters. We then discuss the extension to
non-Gaussian distributions via the copula transform andsummarize the algorithm.
A. The pure SGCRF model
The SGCRF models y|x as a multivariate Gaussian distribution in terms of its exponential family form, via precision(inverse covariance) and mean parameters rather than acovariance and mean. The advantage of this formulation issparsity: in many forecasting domains, virtually all variables
share some degree of correlation (for instance, the windten hours from now may be correlated with the wind one
hour from now, in a purely statistical sense) and thus
modeling the full joint distribution requires a large numberof parameters. However, it is well-known that the inverseof the covariance matrix for a Gaussian captures just theconditional independencies of the distribution and can behighly sparse even when the covariance matrix is dense. In
the non-conditional case, this fact has been exploited byseveral authors to formulate efficient convex methods forsparse inverse covariance estimation using the `1 penalty [2],also known as the graphical lasso [5].
For prediction tasks such as probabilistic forecasting, it
has repeatedly been observed that a discriminative approachmodeling the conditional distribution (y|x) can be superior
to the generative approach of modeling the full distribution (y, x) [12]. Thus, we extend sparse inverse covarianceestimation to the discriminative case, jointly modeling thecorrelations in the output variables along with their conditional dependence on input features. Since the discriminative
setting coincides with the standard notion of a conditionalrandom field [18], we refer to this model as the sparse
Gaussian conditional random field (SGCRF).
Formally, the SGCRF models the distribution y|x as


1 T
1
T
exp − y Λy − x Θy(5)p(y|x; Θ, Λ) =Z(x)2where Λ ∈ Rp×p and Θ ∈ Rn×p are the parameters of themodel, and Z(x) is the partition function (used to ensure the
distribution over y integrates to one) given by


1
1 T
1/2
−1
= c|Λ| exp − x ΘΛ Θx
(6)
Z(x)
2

where c is a constant term that is independent of x, Λ, andΘ. Critically, we can express models with high correlation
between variables even though Λ and Θ are sparse. To
see this, note that the model can easily be transformed tomean/covariance formp(y|x) ∼ N (−Λ−1 ΘT x, Λ−1 )(7)but Λ−1 Θ and Λ−1 are likely dense even when Λ and Θ
are sparse. Thus, in the forecasting setting, each elementof our prediction zt+1:t+Hf can depend on every element of
zt−Hp +1:t and wt . By exploiting sparsity, we learn the model
efficiently from much less data than would be required toestimate the mean and covariance directly.
B. A second-order optimization method for the SGCRF
Next, we describe an efficient second-order optimization
method for the SGCRF which exploits the sparse structure of
the solution with an active set approach. In particular, given
a collection of input-output pairs {xi , yi }, i = 1, . . . , m, weestimate the parameters Λ and Θ using maximum likelihoodand encourage sparsity by adding an additional `1 penalty
λ(kΛk1 + kΘk1 ), where λ is a regularization parameter andk · k1 denotes the elementwise `1 norm of a matrix. This
results in the optimization problem

second-order expansion along with the `1 penalty is givenby the solution to the regularized quadratic programd = arg min g(∆) + λkx + ∆k1∆(11)which we solve using coordinate descent. Given this direction, we choose our step size using backtracking line
search and iterate until convergence. Note that our model isparameterized by two matrices Λ and Θ and thus the descentdirection is pair of matrices, ∆Λ and ∆Θ .
Using coordinate descent for computing the descent direction is particularly appealing since it allows us to exploit
sparsity in the solution by maintaining a small active set.
In particular, even though Λ and Θ are comprised of p(p +1)/2 + np different variables, due to the sparsity inducedby the `1 penalty and the nature of forecasting applications,
we expect the majority of these coordinates to be zero. We
exploit this fact by considering only the set of variables
that violate the optimality conditions; optimizing over Λij(respectively Θij ) only if|(∇Λ f (Λ, Θ))i,j | > λ or Λij 6= 0|(∇Θ f (Λ, Θ))i,j | > λ or Θij 6= 0.

(12)

It can be shown that this heuristic produces an algorithm
that converges to the optimal solution even though theminimize log |Λ| + tr ΛSyy + 2 tr ΘSyx +descent direction at each step is not necessarily optimal. InΛ,Θ(8) practice, for problems with a sparse solution, this resultstr Λ−1 ΘT Sxx Θ + λ(kΛk1 + kΘk1 )in a substantial speed increase over the naive approach of
considering every variable.
where Syy , Syx and Sxx are sample covariance matrices
Finally, in order to make coordinate descent efficient,
m
m
m
X
X
it
is critical to cache and iteratively update certain matrix
11
1 X
yi yiT , Syx =yi xTi , Sxx =xi xTi . products. The gradient and Hessian of our objective are
Syy =
m i=1
m i=1
m i=1
−1 T
(9) quite involved due to the tr Λ Θ Sxx Θ term, however
Although this is a convex problem, existing solvers suffer performing the coordinatewise updates efficiently reduces−1and ∆Θ Λ−1 . At
from slow convergence due in part to the coupling of the to maintaining the product of ∆Λ Λparameters in the matrix-fractional term tr Λ−1 ΘT Sxx Θ; in each step in the coordinate descent inner loop, we reusepractice, forecasting over long time horizons across multiple these products and when we update a single coordinate,
locations quickly becomes intractable. Therefore, we develop (∆Λ )ij ((∆Θ )ij ), we must update the corresponding ith−1(respectively ∆Λ Λ−1 ). It can be showna custom second-order method which is several orders of row of ∆Λ Λ
magnitude faster than standard solvers, enabling us to apply that the resulting algorithm achieves superlinear convergence,
this method to large-scale probabilistic forecasting tasks such which in practice allows this model to be applied to manyas forecasting wind power production over several adjacent previously intractable probabilistic forecasting problems.
wind farms and multiple days. Here we give the motivationand highlight the key technical points of the algorithm, but C. Non-Gaussian distributions via copula transforms
for complete details and theoretical analysis see [20].
The above SGCRF is limited in that it can only model the
The main idea is to decompose the objective function distribution over y as a multivariate Gaussian. To overcome
into a smooth term plus regularization term and find the this limitation, we employ a (Gaussian) copula transform
descent direction by iteratively optimizing a regularized ver- [22], a method for converting multivariate Gaussian distribusion of the second-order approximation to the smooth term. tions into multivariate distributions with arbitrary marginalIn general, for an optimization problem with the objective distributions. Previous work has applied the copula transformf (x) + λkxk1 , the second-order Taylor expansion of f is
to extend the sparse Gaussian MRF to non-Gaussian distributions [8] and here we extend this to the SGCRF, forming
1
f (x + ∆) ≈ g(∆) ≡ f (x) + ∇x f (x)T ∆ + ∆T ∇2x f (x)∆a model which is well-suited for probabilistic forecasting in2(10) a wide variety of energy systems.
where ∇x f (x) and ∇2x f (x) denote the gradient and Hessian
Formally, suppose u ∈ R is a univariate random variable
respectively. For a fixed x, finding ∆ that minimizes this with cumulative distribution function (CDF) F ; when weTABLE IC OMPARISON OF PREDICTION ERROR

only have samples of u, we use the empirical CDF
mF̂ (u) =1 X
1{u < ui }.
m i=1

(13)

In the case that we expect the variables to come from known
distribution (e.g. the Weibull distribution for modeling windspeeds), we could use the analytical CDF of this distribution
directly. The copula transform simply converts the sample
distribution to a uniform [0, 1] random variable by the CDF
F , then applies the inverse normal CDF Φ−1 to transformthe [0, 1] random variable into a Gaussian random variable.
Our algorithm models the variables using a SGCRF in thistransformed Gaussian space, and then transforms back to theoriginal distribution by applying the inverse copula transform(the normal CDF Φ followed by the inverse CDF F −1 ).
D. Final Algorithm
As with all learning methods, using SGCRFs consists of a
training stage where we learn the parameters that maximize
the model’s likelihood on past observations. Then, for a new
scenario (denoted x0 ∈ Rn ), we use the model to makepredictions about the future observations y 0 . The training
stage consists of the following elements:
1) Given data (xi , yi ), for i = 1, . . . , m (recall that eachxi consists of Hp past observations and external inputs
wt , and each yi consists of Hf future observations),
first estimate the univariate marginal distributions ofeach (yi )j , denoted Fj .
2) Transform each the yi variables to a variable with
marginal Gaussian distributions ỹi by applying theelementwise copula transform(ỹi )j = Φ−1 (Fj ((yi )j ))(14)3) Train a SGCRF model (i.e., estimate the Θ and Λ
parameters) on (xi , ŷi ), i = 1, . . . , m.
With a model, we can perform any of the following tasks:•

Compute the most likely output: Compute the meanin the Gaussian space ỹ 0 = −Λ−1 ΘT x0 ; then transform
each element of ỹ 0 using the inverse copula transform(ŷ 0 )j = F −1 (Φ((ỹ 0 )j ))

•

(15)

Compute the probability of a given output y 0 : Convert y 0 to the Gaussian space using (14) and computep(y 0 |x0 ) =p(ỹ 0 |x0 ; Θ, Λ)
1 0 T 0
1
0 T
0
exp
−
(ỹ
)
Λỹ
−
(x
)
Θỹ
=
Z(x0 )
2(16)•

Draw a random sample of future observations:
Sampleỹ 0 ∼ N (−Λ−1 ΘT x0 , Λ−1 )(17)
and then apply the inverse copula transform (15) to ỹ 0 .

Algorithm
Linear Regression
Linear Regression + copula
ARMAX
SGCRF
SGCRF + copula

RMSE
0.15600.16360.17140.14880.1584

IV. A PPLICATION TO WIND FORECASTING
In this section, we describe the primary applied result of
this paper, an application of the above probabilistic forecasting method to a real-world wind power prediction task. We
use data from the GEFCom 2012 forecasting challenge, awind power forecasting competition that was recently held on
Kaggle [6], where the goal was to predict power output from
7 nearby wind farms over the next 48 hours using forecasted
wind speed and direction as input variables.
In our setup, we model wind power production jointly
across all wind farms as zt ∈ R7 and include the forecasted
wind at each farm as exogenous variables. We model thenon-linear dependence of the wind power using radial basis
functions (RBFs) with centers and variances tuned using
cross-validation, resulting in 10 RBFs for each time point
and location and wt ∈ R3360 . We also include autoregressive
features for past wind power over the previous 8 hours
which we found experimentally to be sufficient to capture
the autoregressive behavior of wind power in this dataset. In
our framework, the input and output variables (xt , yt ) arecompromised of wt and zt ranging over past and future time
points




zt
zt+1
 .. 

(18)
xt =  .  , yt =  ...  zt−7 
zt+48
wt
resulting in xt ∈ R3416 and yt ∈ R336 .
We fit the model using 80% of the provided data (874training examples) and report results on the held out set. As
baselines, we consider a linear regression model (LR) which
predicts each output independently and an ARMAX modelwith AR(3) and MA(2) components, both using the same
input features as the SGCRF.
A. Probabilistic predictions
Typically, forecasting systems are evaluated solely on thequality of the point forecasts produced and we see in Table I
that, on this basis, the SGCRF method performs significantly
better than linear regression and ARMAX. Due to the highdimensionality of the feature space relative to the numberof training examples, we expect the `1 penalty employedby the SGCRF to be statistically efficient in identifying
the underlying structure of the correlations in wind power
production and the dependence of wind power on windforecasts; indeed in Figure 1, we see that the estimatedparameters exhibit a high degree of sparsity.

Prob(Energy = x)

0.2
0.15
0.1
0.05

Fig. 1. Sparsity patterns Λ and Θ from the SGCRF model. Λ is estimatedto have 1412 nonzero entries (1.2% sparse) and Θ is estimated to have 7714
nonzero entries (0.67% sparse). White denotes zero values and wind farms
are grouped together in blocks.

0.25
Prob(Energy = x)LR
LR + copula
SGCRF + copula

0.25

0.2
0.15
0.1
0.05

0

0
0

50
x

100

0

50
x

100

Fig. 2. Examples of predictive distributions for total energy output fromall wind farms over a single day.

TABLE IIC OVERAGE OF CONFIDENCE INTERVALS
MethodLR

LR + copula

ARMAX

SGCRF

SGCRF + copula

Task
Aggregate farms
Aggregate times
Both
Aggregate farms
Aggregate times
Both
Aggregate farms
Aggregate times
Both
Aggregate farms
Aggregate times
Both
Aggregate farms
Aggregate times
Both

90%
0.69430.37900.19440.7256
0.40280.2176
0.56820.67790.24540.82670.61040.4306
0.89810.87430.879695%
0.76530.44510.25000.80510.46630.26390.64730.81880.31020.87910.68850.53700.94680.92660.925999%0.86000.53640.33330.90400.57280.33800.75700.95440.42130.94430.7976
0.63890.98300.97220.9676

However, we are primarily interested not in the accuracyof point forecasts, but in the ability of the models to capture
the distribution of future power production. Indeed, as shownin the same Table I, the inclusion of the copula transformdegrades the performance of the models in terms of RMSE;
this is expected since by assuming a Gaussian distribution
over the noise, the untransformed models are explicitly
minimizing mean squared error. RMSE alone is a poor
measure of how well an algorithm can actually predict future
observations: if we judge the algorithms by the ability to accurately predict the range of possibilities for future outcomes,a different picture emerges. For example, a natural task fora wind farm operation would be to generate a distribution
over total power produced by all seven wind farms in thenext 24 hours, in order to establish 95% confidence intervals
about the power to be produced; this could in turn by used bystochastic optimal dispatch method, to determine how much
power to generate from other sources.
Table II illustrates the coverage of the confidence intervals
generated by different approaches, evaluated on a held-out
test set of the wind power data. For each example in thetest set, we used each method to generate many samples
of upcoming wind power and for each of these samples,
we computed the total power aggregated over all the farms,
all times, or both, and used these to generate histogramsof the aggregated power. Finally, we used these histograms

to estimate 90%, 95%, and 99% confidence intervals of theaggregate power, and evaluated how often the true total wind
power fell into that interval.
As seen in Table II, the SGCRF + copula model produces
intervals that map very closely to their desired coverage level,
whereas linear Gaussian and ARMAX models perform muchworse. To see why this occurs, we show in Figure 2 severalof these estimated distributions of aggregate power, sampled from the different models. The independent Gaussianmodels are substantially overconfident in their predictions,as multiple i.i.d. random variables will tend to tighten thevariance, leading to vastly inaccurate predictions when thosevariables are in fact highly correlated. We note that we
could also consider a joint linear model by forming theunregularized MLE estimate for the covariance matrix, but ingeneral this is not well-suited for high-dimensional problems
and in fact is undefined for p > m. The ARMAX model
does capture some of the correlation across time via themoving average component and we see in Table II thatit performs significantly better than linear regression whenaggregating predictions across multiple times. However, the
SGCRF with the copula transform clearly achieves the best
results implying that it is more accurately capturing thedisperse nature of the actual joint distribution.
B. Application to ramp detection
One particularly appealing possibility for the probabilistic forecasting methods is in the area of predicting wind“ramps,” times when power experiences a sudden jump froma relatively low to a relatively high value. Because wind
power grows with the cube of wind speed in Region 2 of theturbine operating conditions (before the wind turbines reach
rated power), a small increase in wind speed can lead to alarge change in power, and predicting when these ramps occur is one of the primary open challenges in wind forecasting.
Indeed, a well-known issue with many forecasting methods
is that while they may accurately predict that a ramp will
occur, they are significantly limited in accurately capturing
the uncertainty over where the ramp will occur [4].
Although a detailed analysis of the ramp prediction capabilities of our approach is not the main focus on thispaper, we briefly highlight the potential of our approachin this task. In particular, because the model accurately

ble observations. From a larger perspective, these methods highlight a highly desirable direction for future work:
new machine learning approaches that can produce highdimensional probabilistic forecasts. But to fully exploit thesemodels, we also need stochastic dispatch and optimal power
flow methods than can exploit this uncertainty to mitigate
risks inherent with renewable intermittency. Integrating theseplanning approaches with probabilistic models thus seems
to be a key direction for future work in power systems
operation.

1
Mean
Sample
Truth

0.9
0.8
0.7

Energy

0.6
0.5
0.4
0.3
0.2

R EFERENCES

0.1[1] N. Amjady, F. Keynia, and H. Zareipour. Short-term load forecast of
microgrids by a new bilevel prediction strategy. Smart Grid, IEEE
Transactions on, 1(3):286–294, 2010.
[2] O. Banerjee, L. E. Ghaoui, and A. d’Aspremont. Model selection
through sparse maximum likelihood estimation for multivariate gaussian or binary data. Journal of Machine Learning Research, 9:485–
516, 2008.
[3] P. Brockwell. Time Series Analysis. Wiley Online Library, 2005.
[4] C. Ferreira, J. Gama, L. Matias, A. Botterud, and J. Wang. A survey
on wind power ramp forecasting. Technical report, Argonne National
Laboratory (ANL), 2011.
[5] J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covarianceestimation with the graphical lasso. Biostatistics, 9(3):432–441, 2008.
[6] T. Hong. Global energy forecasting competition, 2012.
[7] B. Lange, K. Rohrig, F. Schlögl, Ü. Cali, and R. Jursa. Wind power
forecasting. Renewable electricity and the grid, pages 95–120, 2008.
[8] H. Liu, J. Lafferty, and L. Wasserman. The nonparanormal: Semiparametric estimation of high dimensional undirected graphs. The Journalof Machine Learning Research, 10:2295–2328, 2009.
[9] J. Mendes, R. Bessa, H. Keko, J. Sumaili, V. Miranda, C. Ferreira,J. Gama, A. Botterud, Z. Zhou, and J. Wang. Development and testingof improved statistical wind power forecasting methods. Technicalreport, Argonne National Laboratory (ANL), 2011.
[10] M. Milligan, M. Schwartz, and Y. Wan. Statistical wind power
forecasting models: results for us wind farms. Technical report,
National Renewable Energy Laboratory, Golden, CO, 2003.
[11] C. Monteiro, R. Bessa, V. Miranda, A. Botterud, J. Wang, G. Conzelmann, et al. Wind power forecasting: state-of-the-art 2009. Technicalreport, Argonne National Laboratory (ANL), 2009.
[12] A. Y. Ng and M. I. Jordan. On discriminative vs. generative classifiers:a comparison of logistic regression and naive Bayes. 2002.
[13] J. Penm, J. Penm, and R. Terrell. The recursive fitting of subset varx
models. Journal of Time Series Analysis, 14(6):603–619, 2008.
[14] P. Pinson. Estimation of the uncertainty in wind power forecasting.
Centre Energétique et Procédés–Ecole des Mines de Paris Rue, 23,2006.
[15] G. Sideratos and N. Hatziargyriou. An advanced statistical method
for wind power forecasting. Power Systems, IEEE Transactions on,
22(1):258 –265, feb. 2007.
[16] K.-A. Sohn and S. Kim. Joint estimation of structured sparsity andoutput structure in multiple-output regression via inverse-covarianceregularization. In Proceedings of the Conference on Artificial Intelligence and Statistics, 2012.
[17] S. A. Soliman and A. M. Al-Kandari. Electrical Load Forecasting:Modeling and Model Construction. Elsevier, 2010.
[18] C. Sutton and A. McCallum. An introduction to conditional random
fields. Foundations and Trends in Machine Learning, 4(4):267–373,
2012.
[19] Various. PJM Manual 19: Load Forecasting and Analysis. PJM,
2012. Available at: http://www.pjm.com/planning/resource-adequacyplanning/ /media/documents/manuals/m19.ashx.
[20] M. Wytock and J. Z. Kolter. Sparse Gaussian conditional random
fields: Algorithms, theory, and application to energy forecasting. In
Proceedings of the International Conference on Machine Learning,
2013.
[21] X.-T. Yuan and T. Zhang. Partial gaussian graphical model estimation.
CoRR, abs/1209.6419, 2012.
[22] I. Žežula. On multivariate gaussian copulas. Journal of Statistical
Planning and Inference, 139(11):3942–3946, 2009.

0
5

10

15

20

15

20

Hour
0.9
0.8
0.7

Energy

0.6
0.5
0.4
0.3
0.2
0.1
0
5

10
Hour

Fig. 3. Samples drawn from the linear regression model (top), and SGCRF
model (bottom)

captures correlations in the predicted observations over time,
if we draw random samples from our model, then we expect
scenarios where the ramp occurs at different times; this is incontrast to most “independent” probabilistic methods, which
would assume a fairly tight distribution over possible futurescenarios. This situation is illustrated in Figure 3, where weshow the mean prediction along with 10 samples drawn fromeach model. Again, in a stochastic optimal control task, an
operator would be much better served by considering thepossible scenarios generated from our model than from theindependent probabilistic model, as they consist of several
different timings for the upcoming power ramp.
V. C ONCLUSION
In this paper, we present a probabilistic forecasting approach based on state-of-the-art methods in machine learning that can efficiently model high-dimensional and nonGaussian joint distributions over its predictions. Such probabilistic forecasting is capable of more accurately capturing
the true distribution and can be an indispensable tool formaking accurate predictions of the actual range of possi-