11. Probabilistic Forecast Methods

11.1 The Different Forms of Probabilistic Forecasts

As introduced in Sect. 5.​2 and Fig. 5.​4, there are three core forms of probabilistic forecasts which will be explored in this book: quantile forecasts, density forecasts and ensemble forecasts (not to be confused with ensemble machine learning models such as random forest in Sect. 10.​3.​2). These can be grouped into two core categories: univariate (quantile and density) and multivariate (the ensemble forecasts). To understand these types, consider the scenario of trying to estimate the distribution of the data for the time steps $t=N+1, N+2, \ldots , N+k$.

For a univariate forecast, the aim is to estimate the distribution of the demand at each time step. In other words, estimate a total of k univariate distributions. The distributions at different time steps are independent of each other, meaning that the spread of the variable at one time step is not influenced by information about the demand at other time steps. Another way to say this is there are no inter-dependencies modelled. The univariate Gaussian, as introduced in Sect. 3.​1 is the most famous example of a univariate distribution function. There are two main forms for estimates of univariate distributions, either a full continuous density function, or discrete values (quantiles) defining uniformly spaced levels of equal probability.¹ These two forms are illustrated in Fig. 11.1 for a standard Gaussian distribution. The density estimate is often preferable, as it describes the entire distribution, but requires either knowing/assuming that the distribution is from a particular parametric family (e.g. Gaussian in this case), or requires training relatively expensive methods, such as Kernel Density Estimates which will be described in Sect. 11.5. The quantile estimates (described in more detail in Sect. 3.​2) in Fig. 11.1, show the $5, 10, \ldots , 95\%$ quantiles for the Gaussian density, and are clearly less descriptive of the distribution, however they are a lot less expensive to compute and don’t rely on assuming a specific distribution of the data. Since a univariate distribution has to be estimated for each of the k time steps in the forecast horizon this reduction in computational cost can be particularly advantageous.

For multivariate forecasts the task is instead to estimate a single multivariate distribution for all k demand variables in the forecast horizon (see Sect. 3.​3 for more on multivariate distributions). The advantage of multivariate distributions is that they take into account the inter-dependencies over the entire forecast time horizon. To illustrate this, consider the example of household demand. This is mainly determined by the occupants behaviour. If a person gets into work late then they will likely get back from work later, hence their demand shift in the morning will correspond to a shift in the evening. In other words, there is an interdependency between the demand in the morning and the demand in the evening due to the link between these two activities. A multivariate forecast can therefore be sampled to produce demand profile scenarios which include these correlations. Thus more complicated and realistic interdependent behaviours can be simulated and utilised to optimise applications such as storage control (Sect. 15.​1).

Similar to the univariate case, the full multivariate density can be estimated but it is typically more complicated and difficult to model accurately. The methods are often more computationally expensive, there is fewer packages/resources for fitting them, and there are very few standard parametric multivariate distribution functions which can be used to fit to the data. Instead finite samples from the distribution are often estimated instead. These ensembles are realisations or ‘representatives’ from the distribution. To illustrate this consider a basic case with $k=2$ where the estimates at two consecutive time steps $t=1, 2$ are jointly described by a bivariate Gaussian distribution, shown in the left of Fig. 11.2. Drawing 30 random samples from this distribution gives the time-dependent correlated bivariate ensembles on the right. Notice in the language of distributions introduced in Sect. 3.​3 that the multivariate probabilistic forecast is a joint distribution and the univariate probabilistic forecasts are marginal distributions of the full joint distribution.

11.2 Estimating Future Distributions

As discussed in the previous section, the aim of a probabilistic forecast is to estimate the future distribution of the demand whether at a single time step (univariate) or multiple (multivariate). To estimate the uncertainty requires accurately modelling the variation. There are a few standard practical approaches, which will be outlined in the section, and are the basis for many of the techniques in the following sections.

The first approach tries to model the distribution directly by training on the observations. As with most point forecasts these models use the historical data to capture the variation and will typically make assumption about how the past distribution will relate to the future demand. The parameteric models (Sect. 11.3), kernel density estimation (Sect. 11.5) and the quantile regression (Sect. 11.4) all model the distribution in this way.

The aim is to train the parameters or hyperparameters of a distribution model directly (e.g. the Gaussian model) or use a model which will estimate the distribution (e.g. quantiles). The advantage of these approaches is that as long as the right model is used, and they are trained on sufficient data from the target distribution, then they can accurately capture the uncertainty. For example, if we are modeling demand for 2pm and we know that the historic 2pm data all come from the same distribution then this data can be used to estimate the true distribution. Unfortunately, it is often not known for certain which data comes from the same distribution and so certain assumptions will need to be made based on the analysis of the data. Another drawback of this approach is that its accuracy is correlated with the amount of available data. Small amounts of data will mean a potentially inaccurate estimate.

The second type of model doesn’t model the variation directly but instead inserts variation into the model either through adjusting the input variables and/or the model parameters. For example, assume there is a demand model which is dependent on temperature alone. Then the variation in the demand can be modelled by inserting different values of the temperature into the model. Usually these are formed by tweaks on an individual estimate of the temperature and simulates the sensitivity of the demand to the temperature.

This approach is used in numerical weather prediction to produce forecast ensembles/scenarios. Small deviations are applied to the most likely state of the atmosphere and the numerical weather prediction models are reapplied to the adjusted states to produce a range of weather scenarios. Analysis of closeness of the final ensembles can indicate confidence in the future weather states, and widely ranging ensembles may mean there the future weather is highly uncertain.

Alternatively small adjustments can be applied to the model parameters. This accounts for mis-specifications in the model and can generate other likely future states. Multiple adjustments can therefore produce a range of outputs allowing for an estimate of the future distribution. The difficultly with both of these adjustment approaches is that the correct deviations have to be applied to the inputs/parameters in order to produce an accurate distribution estimate. This can be aided in the input case by randomly sampling from the historical observations, or from estimating a distribution from which you can sample.

Another drawback of this model is that the demand variation is not being simulated directly but instead is estimating the sensitivity of the model to the inputs or parameters. Consider the temperature example above. The demand may change with the temperature but in fact it is the variation in the demand for a fixed temperature which is of primary interest (assuming the temperature can be accurately forecast). The key is to add adjustments to the temperature so that it captures this variation. Once again cross-validation is one approach which can be used to determine an appropriate adjustment to the inputs/parameters.

The following sections will mainly focus on the first approach for producing probabilistic forecasts and train the models directly on the historical observations.

11.3 Parametric Models

Parametric distribution models are desirable as they can give a full description of the spread of the data usually using only a few parameters. This section begins by discussing parametric models via a simple example of a univariate distribution (Sect. 11.3.1). Individual univariate parametric models are usually too inflexible to model the distributions accurately, but families of simple univariate distributions can be “mixed” to estimate much more general shapes and will be introduced in Sect. 11.3.2.

11.3.1 Simple Univariate Distributions

Some simple univariate distributions have already been introduced in Sect. 3.​1. The most common being the Gaussian (or Normal) distribution, but the lognormal, and the gamma distribution were also presented. A range of distributions can be modelled using these functions in addition to other similar ones. The advantage of such models is that they only require training a small number of parameters to fully estimate the distribution. However, the restriction to a specific functional form means simple parametric models cannot estimate more complex distributions. For example, the distribution functions mentioned above are all unimodal which means they describe distributions with a single modal (maximum) value. It would be impossible to model multimodal distributions (distributions with multiple distinct maximum values). An example of a unimodal, bimodal and trimodal distribution are shown in Fig. 11.3.

Although univariate models are unlikely to produce the most accurate univariate probabilistic forecasts they can be useful as benchmark models to compare to more sophisticated approaches described later in this chapter. Further since they are described by relatively few parameters they may be easier to train than nonparametric models. Training parametric models requires estimating each individual parameter which describes the chosen distribution family. For example, a Gaussian will require estimates for the mean and standard deviation, whereas the gamma distribution requires estimating the shape and scale parameters. In the case of the Gaussian distribution the mean and standard deviation can be found by maximum likelihood estimation (Sect. 8.​2.​1) and these values turn out to simply be the sample mean and sample standard deviation (Sect. 3.​5) respectively. To ensure the best possible estimate is produced requires carefully selecting the most appropriate input data to train the parameters (in contrast to data driven machine learning techniques which will learn from all the data). The data can be identified by the analysis methods outlined in Chap. 6. For example, suppose some hourly data is discovered to have strong daily periodicity then it may be appropriate to train 24 different models, each one using only the data from a specific hour of the day.

Parametric models also exist for multivariate models. In particular there is a multivariate version of the Gaussian distribution. As mentioned in Sect. 11.1 these parametric models can be used to produce ensemble probabilistic forecasts over the h-steps ahead (by estimating a h-dimensional parametric multivariate distribution). Unfortunately, there are much fewer well defined multivariate distributions which can accurately capture a wide variety of probabilistic forecast behaviours. This makes them less suitable compared to more versatile methods which will be introduced in Sects. 11.6 and 11.7 which can also capture interdependencies across time steps in the forecast horizon.

11.3.2 Mixture Models

More versatile univariate distribution can be modelled by combining mixtures of the simple parametric distributions discussed in Sect. 11.3.1. The general form of a PDF for a finite mixture model of a random variable ${\textbf {x}} \in \mathbb {R}^p$ is

$$\begin{aligned} f({\textbf {x}} ) = \sum _{k=1}^K \pi _k g_k({\textbf {x}} , \theta _k), \end{aligned}$$

(11.131)

where $g_k({\textbf {x}})$ are PDF’s usually from a single family (e.g. Gaussian’s) with their own corresponding parameter’s $\theta _k$ (e.g. mean and standard deviation for a Gaussian). The $\pi _k$ are weights which satisfy $\sum _{k=1}^K \pi _k = 1$, and are often called mixing probabilities. Mixture models are often used for clustering, and in this case each PDF defines a distribution of points from one of the clusters, and the weights signify what proportion of the observations are in each cluster.

The most common mixture models for continuous variables use Gaussian components, i.e.

$$\begin{aligned} g_k(\textbf{x}) = \frac{1}{(2\pi )^{N/2} det(\boldsymbol{\Sigma }_k)^{1/2}}exp \left( (\textbf{x}- \boldsymbol{\mu }_k )^T\boldsymbol{\Sigma }_k^{-1}(\textbf{x}- \boldsymbol{\mu }_k) \right) , \end{aligned}$$

(11.132)

with covariance $\mathbf {\Sigma }_k \in \mathbb {R}^{p \times p}$, and mean vector $\boldsymbol{\mu }_k \in \mathbb {R}^p$. This is called a Gaussian mixture model (GMM). A simple example of a GMM ($p=1$) with three clusters is shown in Fig. 11.4 with mixture probabilities of 0.5, 0.25 and 0.25, means of 1, 3, 6 and all the same standard deviation of 1. It is easy to see that more complicated distributions can be estimated by adding more groups/clusters.

Although GMMs have a lot more parameters to fit to the observations then a single Gaussian model they can be solved relatively efficiently via an iterative process called the expectation-maximisation algorithm (EM) which finds an optimal estimate² for the maximum likelihood function (See Sect. 8.​2.​1).

Consider observations, $\mathbf {x_1}, \mathbf {x_2}, \ldots , \mathbf {x_N} \in \mathbb {R}^p$. Without going into the details of the EM-algorithm, the process iterates between an expectation step (E-step), which calculates the expectation of the log-likelihood function with current estimates of the parameters, and the maximisation (M-step) which updates the parameters which maximises the current expected log-likelihood function. For a GMM this translates to the following steps (calculated for each iteration):

1. Calculate the posterior probability $\tau _{ik}$ that each observations $\textbf{x}_i$ belongs to each group $k=1, \ldots , K$,

   $$\begin{aligned} \tau _{ik} = \frac{\pi _k g_k(\textbf{x}_i, \theta _k)}{\sum _{k=1}^K \pi _k g_k(\textbf{x}_i, \theta _k)}. \end{aligned}$$

   (11.133)

   This is the E-step.

2. Update the mixing probabilities $\pi _k^{new} = \sum _{i=1}^N \tau _{k, i}$, for each component $k=1, \ldots , K$.

3. Update the mean for each component $k=1, \ldots , K$,

   $$\begin{aligned} \boldsymbol{\mu }_k = \frac{\sum _{i=1}^N \tau _{k, i}{} {\textbf {x}}_i}{\sum _{i=1}^N \tau _{k, i}}. \end{aligned}$$

   (11.134)

   Note this is a weighted average, weighted based on the membership probabilities.
4. Update the Covariance matrix for each component $k=1, \ldots , K$,

   $$\begin{aligned} \boldsymbol{\Sigma }_k = \frac{\sum _{i=1}^N \tau _{k, i}({\textbf {x}}_i- \boldsymbol{\mu }_k)({\textbf {x}}_i- \boldsymbol{\mu }_k)^T}{\sum _{i=1}^N \tau _{k, i}}. \end{aligned}$$

   (11.135)

   I.e. a weighted version of the sample covariance.

The number of groups, K, is a hyperparameter that must be chosen. Although this could be picked during cross-validation, a likelihood function framing means that information criteria (as introduced in Sect. 8.​2.​2) can also be used to find the most appropriate number of clusters. For different sized clusters calculate the BIC (or AIC). Plotting the BIC against the number of clusters can be used to find the point where increasing the number of clusters shows diminishing returns in terms of the drop in the BIC. This point is the “elbow” point of the plot (and hence why this heuristic is called the“elbow method”) and indicates one choice for a suitable number of clusters. “Suitable” here is a relatively subjective term since there may be several other reasons why different numbers may be more appropriate or useful.

An example of the method is illustrated in Fig. 11.5. Here, the optimal number of clusters is around four since the tangential lines intersect around this value. Tangential lines are often used to make it easier to identify the elbow and hence the number of clusters.

Fitting a GMM is an easy way to estimate a distribution if all the data comes from the same distribution. For a time series this means that the data used for training forms a stationary series. Unfortunately this is unlikely to be the case in general. Different hours of the day may have different distributions and the time series may be dependent on weather, time of year, or a whole host of other variables. Hence, even though the EM-algorithm allows for relatively quick training of the GMM, there may be insufficient data to train several mixture models accurately.

11.4 Quantile Regression and Estimation

The majority of models used for probabilistic load forecasting are nonparametric and are popular because they allow more flexibility in what distributions are being modelled. One of the simplest and most common ways to generate univariate probabilistic forecasts is quantile regression, the subject of this section. One of the advantages of the method is that it is a simple adaption of standard least squares regression.

Consider estimating the q quantiles (See Sect. 3.​2 for introduction to quantiles) for the time steps $t=N+1, N+2, \ldots , N+k$. Popular choices are deciles (10-quantiles) or demi-deciles (20-quartiles) so that the distribution is split into 10 or 20 areas of equal probability respectively.

For quantile regression the principle is identical, except now instead of minimising a least squares cost function, for each quantile $\tau \in \{1, 2, \ldots , q \}$ a set of parameters, $\hat{\boldsymbol{\beta }}_\tau$ must be found which minimise the difference between the model and the observations according to the quantile loss function, i.e.

$$\begin{aligned} \hat{\boldsymbol{\beta }}_\tau =\arg \min _{\boldsymbol{\beta }\in \mathcal {B}} \left( \min \sum _{t=1}^N c_{\tau }(L_t, f_t(\textbf{Z}, \boldsymbol{\beta })) \right) , \end{aligned}$$

(11.137)

where the cost function $c_{\tau }(x, y)$ is defined by

$$ c_{\tau }(x,y) = {\left\{ \begin{array}{ll} \tau (x-y) &{} x \ge y \\ (1-\tau )(y-x) &{} x < y \end{array}\right. }, $$

This is repeated for all quantiles. Recall this is the same pinball loss score introduced in Chap. 7. The process of quantile regression is slightly more complicated than for least squares regression as the cost function isn’t differentiable. However the problem is easily reformulated as a linear programming problem and can be solved very quickly. Quantile regression is only applicable for models, $f_t(\textbf{Z}, \beta )$, which are linear combination of the parameters. This still allows a lot of versatility in the types of relationships that can be modelled.

To illustrate the process consider a very simple example. Generate 400 points from a Gaussian distribution (See Sect. 3.​1) with mean $\mu =2$ and standard deviation $\sigma =3$ to represent a time series of 400 points. Here the y-axis values are the random points and the order in time is simply the order in which they were sampled. The time series is shown in black in Fig. 11.6. Now consider a simple linear model of the form $f_t(\boldsymbol{\beta })=at+b$, (i.e. the parameters $\boldsymbol{\beta }= (a, b)^T$) where there is no other inputs $\textbf{Z}$ since there is no dependencies in this particular model. A quantile regression for this linear model is applied to the time series as in Eq. (11.137) for each deciles (or 10-quantiles). These are shown in Fig. 11.6 as the red dashed lines. Notice in theory, in the limit of increasing numbers of samples, the final quantiles should be flat horizontal lines, and should describe quantiles for a Gaussian distribution with mean $\mu =2$ and standard deviation $\sigma =3$. However, in this case there is a slight gradient since there is only a relatively small amount of data and skews in the sampling can have a large effect on the model fit.

Recall in Chap. 7 that the Probability Integral Transform (PIT) can be used to assess the calibration of a probabilistic forecast. The quantile regression lines should split the data into equal probabilities of observations which would mean 40 (400/10) observations are expected between each of the consecutive deciles. This is shown to be the case in the PIT in Fig. 11.7. Notice in some quantiles there is actually 39 or 41 observations due to the relatively small number of samples. Having a uniform PIT on the training set should be expected when an appropriate model is chosen. The true assessment of the model is, as always, determined by evaluating it on an unseen test set rather than the training set. In addition, for a probabilistic forecast both calibration and sharpness are important properties and therefore the proper scoring functions introduced in Sect. 7.​2 should be used to evaluate forecasts rather than the PIT alone.

Finally, it should be noted that since each quantile is trained independently, sometimes the quantiles may cross over with each other which would obviously be inconsistent since say the $5\%$ quantile may end up higher than the $10\%$ quantile. To prevent this, it is valid to reorder each quantile at each time step to ensure the p percentile is lower than the q percentile, when $p < q$.

11.5 Kernel Density Estimation Methods

Often to help accelerate the optimisation, the variables are normalised (e.g. to [0, 1]) in order to reduce the search space (see Sect. 6.​1.​3). The forecast can be rescaled after the training is complete. As mentioned in Sect. 3.​4 there are options for the different kernels, and different ones can be tested, although often the choice has minimal impact on the accuracy of the forecasts [1].

The different modifications presented here can obviously be combined to create other models. For example the conditional kernel density form shown in Eq. (11.143) can be extended to include a decay factor like in Eq. (11.141) or restrictions can be applied on the inputs like in Eq. (11.142). As with many KDE methods, the drawback is that each modification often increases the training complexity and computational cost.

11.6 Ensemble Methods

This section introduces ensemble forecasts, by which we mean a set of point forecasts from the same forecast origin, estimating each time step with the same forecast horizon (of length h time steps). The point forecasts are samples of equal probability from a h-dimensional multivariate distribution representing the joint distribution over the forecast horizon (see Sect. 3.​3 for more on joint distributions). In other words, each ensemble represents an equally likely load trajectory. The methods described in this section produce these ensembles without needing to produce the full joint distribution.

11.6.1 Residual Bootstrap Ensembles (Homoscedasticity)

This section describes a method for generating ensemble forecast which can be viewed as realisations from a full multivariate probabilistic forecast in the case where the time series is assumed to have fixed variance. To begin consider a 1-step ahead point forecast model, for example, this could be the exponential smoothing model in Sect. 9.​2 or the ARIMA models in Sect. 9.​4. Denote this as $f(L_t| \textbf{Z}, \boldsymbol{\beta })$ which may use previous historical data as well as any explanatory inputs (all described by the set of variables $\textbf{Z}$) to create an estimate, $\hat{L}_{t+1|t}$, for the true value, $L_{t+1}$, at $t+1$. The $\boldsymbol{\beta }$ are the parameters for the model. The next time step ahead can be forecast by iteratively applying the model and including the forecast from the previous time step as a pseudo-historical input to the model. Hence for the next time step

$$\begin{aligned} \hat{L}_{t+2|t} = f(\hat{L}_{t+1|t}| L_t, \textbf{Z}, \beta ) = f(f(L_t| \textbf{Z}, \beta )| L_t, \textbf{Z}, \beta ). \end{aligned}$$

(11.144)

The process can obviously be repeated to produce k-step ahead forecasts. Now, recall that there is an error process describing the difference between the observations and the 1-step ahead forecasts described by the residual $\epsilon _{t+1} = L_{t+1} -\hat{L}_{t+1|t}$. By including the small deviations, described by the residual series, into the forecast, different trajectories can be created which represent different, but equally likely outcomes.

To describe the algorithm in more detail, consider a k-step ahead forecast generated from a one step ahead model, e.g. $\hat{L}_{t+1} = f(L_t| \textbf{Z}, \boldsymbol{\beta })$. Assume the forecast origin is at time step $t=N$. Hence the aim is to produce ensemble forecasts which cover the period $N+1, \ldots , N+k$. The residual series, $\epsilon _t = L_{t+1} - \hat{L}_{t+1}$ (which is calculated for the entire training set), is assumed to have a fixed variance (the series is said to have homoscedasticity) and are uncorrelated with each other. Using this residual series the process of generating a new ensemble for a k-step ahead forecast, using a 1-step ahead forecast model f, is relatively simple. For each ensemble, b, the procedure is as follows:

1. Randomly sample with replacement (this is called a bootstrap sample) a residual, $\hat{e}^{(b)}_1$, from the set of all residuals, $\{\epsilon _1, \epsilon _2, \ldots , \epsilon _N \}$.

2. Add this residual to the current 1-step ahead forecast value $\tilde{L}_{N+1|N}$ to produce a new value $\hat{L}^{(b)}_{N+1|N} = \tilde{L}_{N+1|N}+\hat{e}^{(b)}_1$.

3. Include $\hat{L}^{(b)}_{N+1|N}$ in the forecast model to generate an estimate for the next time step, $\tilde{L}^{(b)}_{N+2|N+1} = f(\hat{L}^{(b)}_{N+1|N}| L_N, \textbf{Z}, \beta )$.

4. Update this value using another bootstrap sample from the residual series to give $\hat{L}^{(b)}_{N+2|N+1} = \tilde{L}^{(b)}_{N+2|N+1}+\hat{e}^{(b)}_2$.

5. Continue this procedure until the kth step is reached.

6. The final series, $\hat{L}^{(b)}_{N+1|N}, \hat{L}^{(b)}_{N+2|N+1}, \ldots , \hat{L}^{(b)}_{N+k|N+k-1}$ is the bth bootstrap ensemble.

This is also known as a residual bootstrap forecast. The process can be repeated to produce as many ensembles as desired. The more ensembles generated, the more load diversity is captured. Generating more samples increases the computational cost, but since each ensemble is independent of the others they can be generated in parallel. Notice that this method strongly assumes that the 1-step ahead errors in the future will be similar to the past 1-step ahead errors. An example of a simple periodic series with homoscedasticity is shown in Fig. 11.9a.

If instead of sampling from the actual residuals you sample from an assumed or fitted distribution then the method can be referred to as a Monte Carlo forecast. For example, it is often assumed that residuals are Gaussian distributed with zero mean and therefore instead of sampling from the set of residuals, the values can be sampled from a Gaussian distribution trained on the residuals. An example of an ensemble forecast generated from the Monte Carlo simulations, for 100 ensembles is shown in Fig. 11.10 for a simple ARIMA(4, 1, 1) model. Notice that the errors get wider (have larger variation) with forecast horizon length. This is due to the accumulation of the errors from one step to the next. This intuitively makes sense as the uncertainty should increase the further ahead the prediction.

Notice that the forecasts at each time step can be used to estimate a univariate estimate. This can be done by either fitting quantiles or a density estimate to the collection of ensemble points at each time step.

11.7 Copula Models for Multivariate Forecasts

This section introduces Copulas, another popular method for generating multivariate probabilistic forecasts, widely used in quantitative finance, but now used extensively in energy forecasts. Here only the basics will be described. Suggested further reading can be found in Appendix D.

A copula is simply a function, C, from an N-dimensional unit box to a 1-dimensional unit box, i.e. $C: [0, 1]^N \longrightarrow [0,1]$, and describes a cumulative distribution function on variables $U_1, \ldots , U_N$ where each variable $U_i$ is uniformly distribution over [0, 1], i.e. $C (u_1, u_2, \ldots , u_N) = P (U_1 \le u_1, U_2 \le u_2, \ldots ,U_N \le u_N)$. In other words a copula models the joint distribution on variables $U_1, \ldots , U_N$ with uniform marginal distributions (See Sect. 3.​3 for more details on joint and marginal distributions). A copula is focused on the correlation/inter-dependence structure of the variables. An advantage of these methods is the wide range of different possible copula functions which can be used to model the inter-dependence.

The process in Eq. (11.147) is reversible which means once a copula model has been trained on observations then multivariate samples can be easily produced with the relevant dependency structure. First samples $u_1, u_2, \ldots , u_N$ are generated from the copula distribution and then each variable is transformed into the original random variable space using the inverse CDF of the marginal $F^{-1}_i(u_i)$ for each corresponding component of the sample point.

Suppose that $\textbf{X} =(X_1, X_2, \ldots , X_N)^T \in \mathbb {R}^N$ has a multivariate Gaussian distribution, then the corresponding copula is defined purely by the correlation matrix since this explains the entire dependency structure. This also means that each correlation matrix defines a specific Gaussian copula. Of course just because a multivariate random variable has a Gaussian copula doesn’t mean that it follows a multivariate Gaussian distribution since the marginals need not be Gaussian. If the correlation matrix is the identity then the copula is called the independence copula defined by

$$\begin{aligned} C_0(u_1, u_2, \ldots , u_N) = u_1\ldots u_N, \end{aligned}$$

(11.148)

where each component is independent of the other components. In general, there is no simple analytic form for a Gaussian copula but it can be expressed as

$$\begin{aligned} C^{Gauss}_R(u_1, u_2, \ldots , u_N) = \Psi _R(\Psi ^{-1}(u_1), \ldots , \Psi ^{-1}(u_N)), \end{aligned}$$

(11.149)

where $\Psi$ is the univariate CDF of the standard Gaussian (mean zero and unit standard deviation) and $\Psi _R$ is the multivariate Gaussian with zero mean and correlation matrix R.

It is clear that different copula’s are useful for different dependency structures and that not all copulas are useful for all types of data. For example, Gumbel copula’s shouldn’t be used with data with negative correlations. How to choose and fit a copula will be briefly considered later in this section.

The value of the Pearson’s correlation coefficient depends on the marginals and the copula, which means random variables will have different Pearson values when transformed using the marginal CDFs. A more convenient measure for the correlation structure for a copula are Rank correlation coefficients which only depend on the rank of the data (see Sect. 3.​5). Since the rank of data is unchanged by the application of a monotonically increasing function, it means the rank correlation coefficients won’t change when the marginal CDFs are applied. An example of a rank correlation coefficient was given in Sect. 3.​5: the Spearman’s rank correlation coefficient.

To better understand how copula’s work and how they can be used to generate new samples, consider a simple bivariate distribution for the random variables (X1, X2), whose dependency structure is described by a Gaussian copula with covariance

$$ R= \begin{bmatrix} 1 &{} 0.8 \\ 0.8 &{} 1 \end{bmatrix} $$

i.e. it has linear Pearson correlation $\rho =0.8$. Also suppose that the marginal of the first variable, X1, has a distribution described by a Gamma distribution, where

$$\begin{aligned} Gamma(X, \alpha , \beta ) = \frac{1}{\beta ^\alpha \Gamma (\alpha )}\int _{0}^{X} t^{\alpha -1} e^{-t/\beta } dt \end{aligned}$$

(11.152)

where $\alpha = 2$ is the shape parameter (determines the shape of the distribution), and $\beta = 1$ is the scale parameter (determines how spread the distribution is). $\Gamma (.)$ is the so-called gamma function. The second variable, X2, is described by the standard Gaussian distribution (mean zero and unit standard deviation). An example of 1000 observations from this distribution is shown in Fig. 11.12. On the horizontal and vertical of the plots are marginal histograms of each variable, which shows the Gamma and Gaussian distributions respectively.

Given the copula, new samples can be generated and then transformed to the original space (with the same linear correlation) by using the inverse CDF of the marginals. An example of 1000 new points generated using the copula (and transformed using the inverse CDFs) is shown in Fig. 11.14. Notice how the final distribution successfully resembles the original distribution in Fig. 11.12.

A similar process can be used to generate ensemble demand forecasts. In this case, consider a demand time series $L_1, L_2, \ldots ,$ where the aim is to generate a day-ahead multivariate forecast with, say, forecast origin $t=N$. Further for simplicity, suppose the data is hourly and hence a 24-step ahead forecast is being considered. The aim here will be to generate a multivariate distribution for the day, $F_{N+1, \ldots , N+24}(L_{N+1}, \ldots , L_{N+24})$, and hence model the inter-dependencies between different times of the day. It is assumed that CDFs for the marginals at different times of the day are already known, i.e. $F_1(L_{N+1}), \ldots , F_{24}(L_{N+24})$ are known. These could be estimated, e.g. by the univariate probabilistic models described earlier in this chapter. In this case a copula can be used to model the intra-day dependency structure by training on the daily profiles transformed by the marginals.

Choosing the correct copula’s depends on many factors and a detailed investigation is beyond the scope of this book. Further reading is suggested in Appendix D. In summary, the choice depends on the type of correlation being modelled as well as the dependencies within the tails/extremes of the distribution. One possible approach for choosing an appropriate copula(s) model can be based on comparison on a validation set as described in Sect. 8.​1.​3.

11.8 Questions

For the questions which require using real demand data, try using some of the data as listed in Appendix D.4. Preferably choose data with at least a year of hourly or half hourly data. In all the cases using this data, split into training, validation and testing in a $60, 20, 20\%$ split (Sect. 8.​1.​3).