# Additive Covariance Matrix Models: Modelling Regional Electricity Net-Demand
in Great Britain

V. Gioia$\,{}^{1}$, M. Fasiolo$\,{}^{2}$, J. Browell$\,{}^{3}$, R.
Bellio$\,{}^{1}$
$\,{}^{1}\,$University of Udine, Department of Economics and Statistics
$\,{}^{2}\,$University of Bristol, School of Mathematics
$\,{}^{3}\,$University of Glasgow, School of Mathematics and Statistics
<EMAIL_ADDRESS>
###### Abstract

Forecasts of regional electricity net-demand, consumption minus embedded
generation, are an essential input for reliable and economic power system
operation, and energy trading. While such forecasts are typically performed
region by region, operations such as managing power flows require spatially
coherent joint forecasts, which account for cross-regional dependencies. Here
we forecast the joint distribution of net-demand across the 14 regions
constituting Great Britain’s electricity network. Joint modelling is
complicated by the fact that the net-demand variability within each region,
and the dependencies between regions, vary with temporal, socio-economical and
weather-related factors. We accommodate for these characteristics by proposing
a multivariate Gaussian model based on a modified Cholesky parametrisation,
which allows us to model each unconstrained parameter via an additive model.
Given that the number of model parameters and covariates is large, we adopt a
semi-automated approach to model selection, based on gradient boosting. In
addition to demonstrating that adopting a covariate-dependent covariance
matrix model leads to substantial forecasting performance improvements,
comparable to those obtained by using a full rather than a diagonal static
covariance matrix, we explore the model output via accumulated local effects
and other visual tools to get insights into how the covariates affect net-
demand variability and dependencies.

The code for reproducing the results in this paper is available at
https://doi.org/10.5281/zenodo.7315106.

_Keywords_ : Multivariate Electricity Net-Demand Forecasting; Covariance
Regression Modelling; Modified Cholesky Decomposition; Generalized Additive
Models.

##  1 Introduction

Electricity networks are changing from centralised systems, where power is
generated by large power plants connected to the transmission network and
consumed mostly on the distribution network, to decentralised networks where
significant generation and storage is connected directly to distribution
networks. The growth of this embedded generation means that the transmission
network now needs to serve the net-demand of customers, that is their demand
net of local production. In Great Britain (GB), embedded production comes
mostly from domestic and small-to-medium commercial solar and wind farms, as
well as small thermal power plants. The lack of visibility of these units at
the transmission level, combined with the weather-dependent nature of
renewable generation leads to considerable challenges in energy trading and
power system operation (Huxley et al., 2022).

The purpose of this work is to support such key operations by proposing an
_interpretable_ modelling approach that provides probabilistic, _spatially
coherent_ short-term net-demand forecasts. The energy industry is conservative
by nature due to the need to maintain security of supply (Von Meier, 2006, see
Chapter 8 of ). As a result, new processes will only be adopted if they are
trusted, and interpretability plays an important role in building trust.
Further, interpretability is of critical importance when extrapolation is
required, for example during exceptional events such as extreme temperatures,
or when the model’s predictions must be decomposed into the contribution of
several effects.

Predicting power flows on the electricity transmission network is a key
motivating application for probabilistic, spatially coherent modelling in
energy forecasting. This is important for both network operators, who are
responsible for system security, and traders who must be aware of spatial
variation in prices. Power flows are influenced by the injection and offtake
of power from the network, as well as network configuration. They are also
constrained by the physics of the network and must be forecasted to identify
and mitigate any risk of exceeding thermal or stability limits. Therefore,
spatial probabilistic forecasts of supply and demand are required to forecast
power flows, and quantify uncertainty and risk associated with these
constraints. Further, as the configuration of the network may change, any
forecasting system must be flexible enough to allow the aggregation of supply
and demand on the fly to calculate flows across relevant boundaries (Tuinema
et al., 2020).

Motivated by the need for probabilistic joint demand forecasts, we consider
joint modelling of net-demand across the 14 regions constituting GB’s
transmission network, which are shown in Figure 1. The net-demand in each
region is the aggregate of the net-demand across many Grid Supply Points
(GSPs), the latter being the interfaces between the transmission system and
either a distribution network or a high-voltage consumer. Correctly modelling
the dependency structure between regions is critical as an error in such a
structure would pollute downstream predictions of power flows. Further, joint
probabilistic forecasts of regional net-demand can be flexibly post-processed
to produce forecasts tailored to the needs of different analyses. For example,
when sustained winds blow in the North of the country, the boundary between
Scotland and the North of England is of particular interest. Indeed, the power
flows across this boundary can be substantial, their direction and size being
driven mostly by wind generation in Scotland, where embedded wind capacity far
exceeds regional demand. Similarly, when the sun shines across the country,
power flows from the South of England, where most of the country’s solar
generation units are installed, to London and the Midlands. During winter peak
hours, the boundaries between London and its surroundings are characterised by
heavy power flows directed toward the Capital, driven by high demand and low
generation capacity in the city.

Figure 1: A map of the GSP groups forming GB’s electricity grid. The letters
are the designation used by the electricity market in GB, while the numbers
correspond to the position of each GSP group in the response vector,
$\mathbb{y}_{i}$ (see Section 3.2). The colours represent five macro-regions
namely Scotland (red), Northern (green), Midlands (light blue), Southern
(blue) and London (yellow). The data on GSPs boundaries has been obtained from
data.nationalgrideso.com/system/gis-boundaries-for-gb-grid-supply-points.

While National Grid’s 2021 ten year statement (National Grid, 2021) gives a
detailed description of tens of transmission grid boundaries and explains
under which circumstances the power lines crossing them become heavily loaded,
the examples given above are meant to convey the fact that power flow analysis
on operational time scales needs regional net-demand forecasts that can be
aggregated, or more generally post-processed, to match the particular scenario
of interest. For concreteness, in this work we consider the aggregation into
the five GSP macro-regions shown in Figure 1, which are motivated by the
boundaries mentioned above and match closely the critical boundaries presented
in a Regional Trends and Insights report from the National Grid (National
Grid, 2018, see Figure 1 therein).

Figure 2: Conditional standard deviations (nodes) and correlations (edges)
across the 14 GSP regions (a–b) or macro-regions (c to e), predicted by the
model from Section 3.3. The plots correspond to 7am on 31/12/18 (a), midnight
on 20/08/18 (b) and 10am on 14/06/18 (c to e). Plot d is based on regional
wind and precipitation forecasts, while c and e correspond to, respectively, a
25% decrease and increase of such forecasts.

Forecasts of power flows, which are composite variables of demand and supply,
among many other factors, can only be calculated if a multivariate predictive
distribution of said quantities is available. Hence, joint forecasts of net-
demand across the GSP regions shown in Figure 1, or across an appropriate,
scenario-dependent aggregation of them, are essential to help the network
operator to take early action when managing the risk of breaching constraints.
However, the structure of spatial dependency in net-demand is complex, as it
is influenced by both socio-economic and weather effects, and is time-varying.
Figure 2a-b illustrates the issue. In particular, the nodes and edges of the
network represent, respectively, the conditional regional standard deviations
and inter-regional correlations of net-demand, predicted one day-ahead by one
of the models proposed in this paper. Figure 2a corresponds to New Year’s Eve,
a day where net-demand forecast uncertainty is particularly high and
correlation is strong between densely populated areas, such as London and the
West Midlands. In contrast, Figure 2b corresponds to the 20th of August and
show that net-demand is predicted to lead to a quieter day, from a network
management perspective, with weak spatial dependency in forecast uncertainty.

Figure 2a-b makes clear that capturing the time-varying nature of regional
net-demand dynamics is essential to produce operationally useful joint
forecast. However, several other factors affect the joint distribution of
regional net-demand, in addition to daily and yearly seasonalities. For
example, the right column of Figure 2 shows the joint configuration of the
macro-regional net-demand variabilities and dependencies during storm Hector.
In particular, Figure 2d shows the prediction obtained by conditioning on the
seasonal factors and weather forecasts corresponding to this time period.
Figures 2c–e have been obtained by respectively decreasing and increasing the
regional wind speed and precipitation forecasts by 25%. They show that weather
has a strong effect on the joint distribution of net-demand. Specifically, a
strengthening of the storm is predicted to lead to a higher variability in
Scotland and to stronger correlations between the latter and the other macro-
regions.

Having motivated the need for an interpretable covariance modelling framework
able to provide the spatially coherent net-demand forecasts required by power
flow analysis, we now outline the modelling approach proposed here. We jointly
model GB regional net-demand using a multivariate Gaussian model, based on a
covariance matrix parametrisation that allows us to model each of its
unconstrained parameters via a separate additive model, containing both
parametric and smooth spline-based effects. In particular, the covariance
matrix of the regional net-demand vector is parametrised via the modified
Cholesky decomposition (MCD) of Pourahmadi (1999). The wiggleness of the
smooth effects is controlled via smoothing penalties, the strength of the
latter being controlled via smoothing parameters. Model fitting is performed
via two nested iterations, the regression coefficients being estimated via
maximum a posteriori (MAP) methods, while the smoothing parameters are
selected by maximising a Laplace approximation to the marginal likelihood
(LAML).

The proposed model can be seen as a multi-parameter generalized additive model
(GAM, Hastie and Tibshirani, 1987) or as a generalized additive model for
location, scale and shape (GAMLSS, Rigby and Stasinopoulos, 2005). Given that
the distinction between scale and shape parameters is unclear under the chosen
distribution and parametrisation, henceforth we will refer to the model simply
as a GAM. Additive models are popular modelling tools in the electricity
demand forecasting (see, e.g., Fan and Hyndman, 2012), in part because they
strike a balance between predictive performance and interpretability, the
importance of the latter in this context having been discussed above. While
ensuring interpretability is challenging under the richly parametrised model
considered here, the MCD parametrisation provides some degree of
interpretability when the response vector has some, not necessarily unique,
intrinsic ordering, as is the case for regional net-demand. We further enhance
interpretable exploration of the model by summarising its output via the
accumulated local effects (ALEs) of Apley and Zhu (2020). Such effects are
aimed at quantifying and visualising the effect of one or several covariates
on a model output, such as the predicted correlation of net-demand between two
regions in our context. While Apley and Zhu (2020) consider the use of such
effects under black-box models, the asymptotic posterior approximation used
here allows us to quantify the uncertainty of the ALEs, by propagating
posterior parameter uncertainty via a standard asymptotic approximation.

The joint regional net-demand model proposed here has 119 distributional
parameters, controlling the mean vector and the covariance matrix of a
conditional multivariate Gaussian distribution. Each parameter can be modelled
via parametric and smooth effects of several covariates, hence the space of
possible models is large. While the effects controlling the mean vector can be
chosen on the basis of expert knowledge or previous research, manual selection
of an additive model for each of the remaining 105 parameters is unrealistic.
Here we leverage the interpretation of the MCD’s parameters to choose the set
of candidate effects that could be used to model each parameter. Then, we use
gradient boosting (Friedman, 2001) to order the effects on the basis of how
much they improve the fit, and we choose the number of effects modelling the
MCD on the basis of forecasting performance on a validation set. The results
show that the semi-automatic effect selection procedure just outlined leads to
satisfactory predictive performance and to model selection decisions that are
largely in agreement with intuition (e.g., wind speed and solar irradiance are
selected to model net-demand variability in, respectively, Scotland and the
South of England).

To our best knowledge, this is the first applied statistical paper to consider
full additive modelling of response vectors that are both more than two- or
three-dimensional and have heterogeneous elements (e.g., they are not lagged
values of the same variable). Additive modelling of multivariate responses has
been proposed by Klein et al. (2015), who consider bivariate Gaussian and
t-distributions based on a variance-correlation decomposition. Marra and
Radice (2017) propose a fitting framework where bivariate responses are
modelled via copulas with continuous margins, and all distributional
parameters are modelled additively. Copula-based models are considered also by
Hans et al. (2022) but, while the example considered therein features a
bivariate response, they propose model fitting via the gradient boosting
methods of Thomas et al. (2018), which can handle higher-dimensional
responses. The same fitting method is used by Strömer et al. (2022a), who
consider several two-dimensional response models, including the bivariate
Gaussian one.

In a generalized linear modelling (GLM) context, Pourahmadi (1999) uses the
MCD to parametrise a multivariate Gaussian model in eleven dimensions. They
are interested in capturing temporal dependencies in longitudinal data, which
allows them to impose a strong structure on the covariance model. Bonat and
Jørgensen (2016) propose methods to fit GLMs to multivariate response and
covariate-dependent covariance matrices, and apply their methods beyond the
bivariate case. They do not use an unconstrained parametrisation, but model
the covariance matrix directly and tune the optimiser to avoid generating
indefinite matrices. Browell et al. (2022) model covariance matrices beyond
two dimensions by using covariance functions whose parameters are allowed to
vary with covariates such as the forecast lead time. When positive
definiteness is not guaranteed by the model formulation, they do not enforce
any constraint but post-process the estimated covariance matrix to find the
nearest positive definite matrix.

From a methodological point of view, our work is closely related to Muschinski
et al. (2022), who consider distributional regression with multivariate
Gaussian responses. They propose non-parametric modelling of the elements of
several Cholesky-based parametrisations, including MCD. However, they fit the
model using Markov chain Monte Carlo methods, rather than direct optimisation
methods as done here and, while they compare a set of manually-chosen
covariance matrix models on a weather forecasting application, here we
consider semi-automatic variable selection to handle a much larger set of
candidate covariates. Further, they aim at capturing temporal rather than
spatial dependencies, the latter being the focus of the present work.

The rest of the paper is structured as follows. Section 2 introduces, in a
general setting, the proposed multivariate Gaussian model structure and
fitting methodology. It also explains how ALEs can be used to quantify the
effect of the covariates on the estimated covariance matrix elements. Section
3 focuses on the regional net-demand modelling application. In particular, the
data is introduced in Section 3.1, while Section 3.2 describes the bespoke,
boosting-based model selection approach proposed here. The output of the final
model is explored in Section 3.3, while the forecasting performance of the
proposed model is assessed in Section 3.4. Section 4 summarises the main
results.

##  2 Multivariate Gaussian Additive Models

###  2.1 Model Structure

Let $\mathbb{y}_{i}=(y_{i1},\ldots,y_{id})^{\top}$, for $i=1,\ldots,n$, be
independent response vectors, normally distributed with mean
$\boldsymbol{\mu}_{i}$ and covariance matrix ${\boldsymbol{\Sigma}}_{i}$. The
$q=d+d(d+1)/2$ unique elements of $\boldsymbol{\mu}_{i}$ and of (a suitable
parametrisation of) $\boldsymbol{\Sigma}_{i}$ are modelled via
$\boldsymbol{\eta}_{i}$, a $q$-dimensional vector of linear predictors. The
$j$-th element of $\boldsymbol{\eta}_{i}$ is modelled via

$\eta_{ij}={\mathbf{Z}^{j}_{i}}^{\top}\boldsymbol{\psi}_{j}+\sum_{l}f_{jl}({\boldsymbol{x}}^{S_{jl}}_{i})\,\,,\;\;\;\text{for}\;\;\;j=1,\dots,q,$
(1)

where ${\mathbf{Z}^{j}_{i}}^{\top}$ is the $i$-th row of the design matrix
$\mathbf{Z}^{j}$, $\boldsymbol{\psi}_{j}$ is a vector of regression
coefficients, $\boldsymbol{x}_{i}$ is an $s$-dimensional vector of covariates
and $S_{jl}\subset\\{1,\dots,s\\}$. Hence, for example, if $S_{jl}=\\{2,4\\}$
then ${\boldsymbol{x}}_{i}^{S_{jl}}$ is a two dimensional vector formed by the
second and fourth element of $\boldsymbol{x}_{i}$. Each $f_{jl}$ is a smooth
function, built via

$f_{jl}(\boldsymbol{x}^{S_{jl}})=\sum_{k}b^{jl}_{k}(\boldsymbol{x}^{S_{jl}})\alpha_{k}^{jl}\,\,,$
(2)

where $b^{jl}_{k}$ are spline basis functions of dimension
$\text{card}(S_{jl})$, while $\alpha_{k}^{jl}$ are regression coefficients.
Denote with $\boldsymbol{\alpha}$ the vector of all such coefficients in the
model. The wiggliness of the effects is controlled by an improper multivariate
Gaussian prior on $\boldsymbol{\alpha}$. The prior is centered at the origin
and its precision matrix is ${\bf
S}^{\boldsymbol{\lambda}}=\sum_{u}\lambda_{u}{\bf S}_{u}$, where the ${\bf
S}_{u}$’s are positive semi-definite matrices and
$\boldsymbol{\lambda}=(\lambda_{1},\lambda_{2},\dots)^{\top}$ is a vector of
positive smoothing parameters. See Wood (2017) for a detailed introduction to
GAMs, smoothing splines bases and penalties.

Let us temporarily drop index $i$ to simplify the notation. In this work we
use the following parametrisation of $\boldsymbol{\mu}$ and
$\boldsymbol{\Sigma}$ in terms of $\boldsymbol{\eta}$: $\mu_{j}=\eta_{j}$ for
$j=1,\dots,d$, while the remaining elements of $\boldsymbol{\eta}$ parametrise
a MCD of $\boldsymbol{\Sigma}^{-1}$ (Pourahmadi, 1999). In particular,

$\boldsymbol{\Sigma}^{-1}=\mathbb{T}^{\top}\mathbb{D}^{-2}\mathbb{T}\,,$ (3)

where $\mathbb{D}^{2}$ is a diagonal matrix with ${\rm
D}^{2}_{jj}=\exp(\eta_{j+d})$, for $j=1,\dots,d$, and

$\mathbb{T}=\begin{pmatrix}1&0&0&\cdots&0\\\ \eta_{2d+1}&1&0&\cdots&0\\\
\eta_{2d+2}&\eta_{2d+3}&1&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots&\\\
\eta_{q-d+2}&\eta_{q-d+3}&\cdots&\eta_{q}&1\\\ \end{pmatrix}.$ (4)

Note that parametrisation (3) is unconstrained, that is, the resulting
covariance matrix ${\boldsymbol{\Sigma}}$ is positive definite for any finite
${\boldsymbol{\eta}}$, which facilitates the fitting process. Other
unconstrained parametrisations could have been used, such as those discussed
by Pinheiro and Bates (1996) and Pourahmadi (2011). However, the MCD approach
is particularly attractive in the context of this work. First, the fitting
methods described in Section 2.2 require the first two derivatives of the log-
likelihood w.r.t. ${\boldsymbol{\eta}}$ and, under the MCD parametrisation,
the multivariate Gaussian log-likelihood can be written directly in terms of
${\boldsymbol{\eta}}$, which facilitates the computation of such derivatives.
Second, the MCD has a regression-related interpretation, which can be
exploited when the response vector has some intrinsic ordering. In particular,
assume w.l.o.g. that $\mathbb{E}(\mathbb{y})=\mathbb{0}$ and that $\mathbb{y}$
follows the regression models

$y_{l}=\sum_{k=1}^{l-1}\phi_{lk}y_{k}+\epsilon_{l}\,\,,\quad\text{for}\quad
l=2,\dots,d,$

where $y_{1}=\epsilon_{1}$, $\text{var}(\epsilon_{l})=\sigma_{l}^{2}$ and
$\text{cov}(\epsilon_{l},\epsilon_{k})=0$ for $l\neq k$. Pourahmadi (1999)
shows that ${\rm T}_{lk}=-\phi_{lk}$ and ${\rm D}_{kk}^{2}=\sigma_{k}^{2}$,
for $k=1,\dots,d$, and $l=k+1,\dots,d$. Hence, the elements of $\mathbb{T}$
can be interpreted as the regression coefficients of the elements of
$\mathbb{y}$ on their predecessors. In Section 3.2 we will discuss how such an
interpretation can facilitate the development of a multivariate model for GB
regional net-demand.

###  2.2 Model Fitting

Let us indicate the set of all response vectors
$\mathbb{y}_{1},\dots,\mathbb{y}_{n}$ simply with $\mathbb{y}$ and with
${\boldsymbol{\beta}}$ the vector of all regression coefficients in the model,
which include ${\boldsymbol{\alpha}}$ and all the unpenalised coefficients
vectors $\boldsymbol{\psi}_{j}$. Let $\tilde{\bf S}^{\boldsymbol{\lambda}}$ be
the prior precision matrix of ${\boldsymbol{\beta}}$, that is an enlarged
version of ${\bf S}^{\boldsymbol{\lambda}}$ padded with zeros so that
${\boldsymbol{\alpha}}^{\top}{\bf
S}^{{\boldsymbol{\lambda}}}{\boldsymbol{\alpha}}={\boldsymbol{\beta}}^{\top}\tilde{\bf
S}^{\boldsymbol{\lambda}}{\boldsymbol{\beta}}$. Then, up to an additive
constant that does not depend on ${\boldsymbol{\beta}}$, the Bayesian
posterior log-density of the model from Section 2.1 is

$\mathcal{L}({\boldsymbol{\beta}})=\log
p({\boldsymbol{\beta}}|\mathbb{y},\boldsymbol{\lambda})=\sum_{i=1}^{n}\log
p(\mathbb{y}_{i}|{\boldsymbol{\beta}})-\frac{1}{2}{\boldsymbol{\beta}}^{\top}\tilde{\bf
S}^{\boldsymbol{\lambda}}{\boldsymbol{\beta}},$ (5)

where $\log p(\mathbb{y}_{i}|{\boldsymbol{\beta}})$ is the $i$-th log-
likelihood contribution.

For fixed smoothing parameters, $\boldsymbol{\lambda}$, we obtain MAP
estimates of the regression coefficients by maximising the log-posterior (5),
using Newton’s algorithm. The latter requires the gradient and Hessian of the
log-posterior w.r.t. ${\boldsymbol{\beta}}$, which are provided in the
Supplementary Material A (henceforth SM A). The real challenge is selecting
the smoothing parameters themselves. We do it by maximising an approximation
to the log marginal likelihood, $\mathcal{V}({\boldsymbol{\lambda}})=\log\int
p(\mathbb{y}|{\boldsymbol{\beta}})p({\boldsymbol{\beta}}|{\boldsymbol{\lambda}})d{\boldsymbol{\beta}}$.
In particular, we consider a LAML criterion

$\tilde{\mathcal{V}}({\boldsymbol{\lambda}})=\mathcal{L}(\hat{{\boldsymbol{\beta}}})+\frac{1}{2}\log|\tilde{\bf
S}^{{\boldsymbol{\lambda}}}|_{+}-\frac{1}{2}\log|\boldsymbol{\mathcal{H}}|+\frac{M_{p}}{2}\log(2\pi)\,\,,$
(6)

with $M_{p}$ being the dimension of the null space of $\tilde{\bf
S}^{\lambda}$, $|\tilde{\bf S}^{\lambda}|_{+}$ the product of its positive
eigenvalues, $\hat{{\boldsymbol{\beta}}}$ the maximiser of
$\mathcal{L}({\boldsymbol{\beta}})$ and $\boldsymbol{\mathcal{H}}$ its
negative Hessian, evaluated at $\hat{{\boldsymbol{\beta}}}$.

We maximise $\tilde{\mathcal{V}}({\boldsymbol{\lambda}})$ via the generalized
Fellner-Schall method of Wood and Fasiolo (2017), under which the $u$-th
smoothing parameter is updated using

$\lambda^{\text{new}}_{u}=\frac{{\rm
tr}\\{(\tilde{\mathbb{S}}^{\boldsymbol{\lambda}})^{-}\tilde{\mathbb{S}}_{u}\\}-{\rm
tr}(\boldsymbol{\mathcal{H}}^{-1}\tilde{\mathbb{S}}_{u})}{\hat{{\boldsymbol{\beta}}}^{\top}\tilde{\mathbb{S}}_{u}\hat{{\boldsymbol{\beta}}}}\lambda_{u}^{\text{old}}\,\,,$
(7)

where $(\tilde{\mathbb{S}}^{\boldsymbol{\lambda}})^{-}$ is the Moore-Penrose
pseudoinverse of $\tilde{\mathbb{S}}^{\boldsymbol{\lambda}}$ and
$\tilde{\mathbb{S}}_{u}$ is ${\mathbb{S}}_{u}$ after padding it with zeros.
That is, if we indicate with ${\boldsymbol{\beta}}_{u}$ the subvector of
${\boldsymbol{\beta}}$ that is penalised by ${\mathbb{S}}_{u}$, then
${\boldsymbol{\beta}}_{u}^{\top}{\mathbb{S}}_{u}{\boldsymbol{\beta}}_{u}={\boldsymbol{\beta}}^{\top}\tilde{\mathbb{S}}_{u}{\boldsymbol{\beta}}$.
An advantage of update (7) is that it does not require computing the
derivatives of $\tilde{\mathcal{V}}({\boldsymbol{\lambda}})$ w.r.t.
${\boldsymbol{\lambda}}$. In particular, as detailed in Wood et al. (2016),
computing the gradient of $\tilde{\mathcal{V}}({\boldsymbol{\lambda}})$
requires the third derivatives of log-likelihood w.r.t. each element of
${\boldsymbol{\eta}}$, which leads to computational effort of order
$O\left\\{n\binom{q+2}{3}\right\\}$ ($\approx 2\times 10^{10}$, if $q=119$ and
$n\approx 8\times 10^{4}$ as in the application considered here). Hence, for
moderately large dimension $d$ of the response vector, a quasi-Newton
iteration for maximising $\tilde{\mathcal{V}}({\boldsymbol{\lambda}})$ would
be too computational intensive, at least under naïve evaluation of the
likelihood derivatives.

###  2.3 Inference and Effect Visualisation

The uncertainty of the fitted regression coefficients, ${\boldsymbol{\beta}}$,
can be quantified using approximate Bayesian methods. In particular, standard
Bayesian asymptotics justify approximating
$p({\boldsymbol{\beta}}|\mathbb{y},{\boldsymbol{\lambda}})$ with a Gaussian
distribution, $N(\hat{{\boldsymbol{\beta}}},{\bf V}_{{\boldsymbol{\beta}}})$,
centered at the MAP estimator and with covariance matrix ${\bf
V}_{{\boldsymbol{\beta}}}=-\boldsymbol{\mathcal{H}}^{-1}$. Such a posterior
approximation does not take into account the uncertainty of the smoothing
parameters estimates, which are considered fixed to the LAML maximiser. Wood
et al. (2016) use a Gaussian approximation to
$p({\boldsymbol{\lambda}}|\mathbb{y})$ and propagate forward the corresponding
smoothing parameter uncertainty to obtain an approximation to the
unconditional posterior, $p({\boldsymbol{\beta}}|\mathbb{y})$. In principle,
this approach could be adopted for the model class considered here, but the
formulae provided Wood et al. (2016) require the Hessian of
$\tilde{\mathcal{V}}$ w.r.t. ${\boldsymbol{\lambda}}$ which involves the
fourth derivative of log-likelihood w.r.t. each element of
${\boldsymbol{\eta}}$. Such derivatives are tedious to derive (see SM A.3)
and, as explained in Section 2.2, computationally expensive to evaluate.

Given that the smooth effects are linear combinations of the regression
coefficients, it is simple to derive pointwise Bayesian credible intervals for
the effects, the asymptotic frequentist properties of such intervals having
been studied by Nychka (1988). However, each effect acts directly on a linear
predictor, the latter being non-linearly related to one or more elements of
${\boldsymbol{\Sigma}}$. As explained in Section 2.1, the MCD parametrisation
is related to a set of regressions involving the elements of the response
vector. This fact aids interpretability only if the response vector has some
natural ordering. While this is to some extent the case in the application
considered here (see Section 3.2), communicating modelling results to non-
statisticians is more likely to be effective if framed in terms of widely-used
concepts such as covariances and correlations, rather than parametrisation-
specific quantities. Hence, we use the accumulated local effects (ALEs) of
Apley and Zhu (2020) to quantify the effect of a covariate on
${\boldsymbol{\Sigma}}$ or on the corresponding correlation matrix,
$\boldsymbol{\Gamma}$. Below we explain what ALEs are and how they are
computed, while in Section 3.3 we will use them to visualise a model for GB
regional net-demand.

Recall that ${\boldsymbol{\Sigma}}_{i}$ depends on the covariate vector
$\boldsymbol{x}_{i}$ via the linear predictor vector
${\boldsymbol{\eta}}_{i}$. In the following we drop index $i$, which is not
needed, and we indicate simply with $\sigma(\boldsymbol{x})$ an arbitrary
element of ${\boldsymbol{\Sigma}}$, its dependence on $\boldsymbol{x}$ having
been made explicit. Assuming that $\sigma(\boldsymbol{x})$ is differentiable
w.r.t. the $k$-th covariate, the main (first order) ALE of $x_{k}$ is

$\bar{\sigma}_{k}(x)=\int_{x_{k}^{\text{min}}}^{x}\mathbb{E}_{{\boldsymbol{x}}_{\backslash
k}}\left\\{\sigma^{k}(z,{\boldsymbol{x}}_{\backslash
k})|x_{k}=z\right\\}dz-\text{const}\,\,,$ (8)

where ${\boldsymbol{x}}_{\backslash k}$ is ${\boldsymbol{x}}$ with the $k$-th
element excluded, $\sigma^{k}=\partial\sigma/\partial x_{k}$ and
$\mathbb{E}_{{\boldsymbol{x}}_{\backslash k}}\\{\cdot|x_{k}=z\\}$ is a
conditional expectation w.r.t. $p({\boldsymbol{x}}_{\backslash k}|x_{k}=z)$.
The choice $x_{k}^{\text{min}}$ is unimportant, as changing it simply shifts
the effect vertically, so in practice $x_{k}^{\text{min}}$ is set to just
below the smallest observed value of $x_{k}$. Formula $(\ref{eq:ALE_1})$
applies to any scalar-valued, differentiable function of ${\boldsymbol{x}}$,
such as an element of $\boldsymbol{\Gamma}$. Let us indicate such a generic
model output with $\omega(\boldsymbol{x})$.

In contrast to the partial dependence plots of Friedman (2001), ALEs avoid
making an extrapolation error when the elements of ${\boldsymbol{x}}$ are
correlated. This is explained by Apley and Zhu (2020), who also provide
formulas for obtaining estimated effects $\hat{\bar{\omega}}_{k}(x)$, by
approximating the integral in (8). While Apley and Zhu (2020) quantify the
uncertainty of ALEs via bootstrapping, Capezza et al. (2021) show that, in a
GAM context, the delta method can be used to approximate their variance via
$\text{var}\\{\hat{\bar{\omega}}_{k}(x)\\}\approx\nabla_{{\boldsymbol{\beta}}}^{\top}\hat{\bar{\omega}}_{k}{\bf
V}_{{\boldsymbol{\beta}}}\nabla_{{\boldsymbol{\beta}}}\hat{\bar{\omega}}_{k}$.
They provide formulas for the Jacobian
$\nabla_{{\boldsymbol{\beta}}}\hat{\bar{\omega}}_{k}(x)$ that apply to any GAM
with multiple linear predictors, covering also the case where $x_{k}$ is a
categorical variable. The only output-specific component is the Jacobian of
the parametrisation linking the output of interest, $\omega(\boldsymbol{x})$,
with ${\boldsymbol{\eta}}$, hence in SM A.4 we provide the Jacobians needed
when $\omega(\boldsymbol{x})$ is an element of ${\boldsymbol{\Sigma}}$ or of
$\boldsymbol{\Gamma}$.

##  3 Joint Multivariate Regional Net-Demand Modelling

###  3.1 Data Description and Modelling Setting

We consider data on regional net-demand in Great Britain spanning five years,
2014 to 2018. We define net-demand as the load measured at the interface
between transmission and distribution networks. In GB these interfaces are
called Grid Supply Points (GSPs), and are grouped into 14 regions called GSP
groups. Let $y_{ij}$, for $i=1,\dots,n$, be the standardised net-demand of GSP
group $j$ measured at a 30min resolution.

In addition to net-demand, the data contains the following covariates or
transformations thereof. Some covariates are common to all GSP groups, these
are: time of day, $\text{tod}_{i}\in\\{0,\dots,23.5\\}$; day of the year,
$\text{doy}_{i}\in\\{1,\dots,366\\}$; day of the week, $\text{dow}_{i}$; time
since the 1st Jan 2014, $\text{t}_{i}$, re-scaled to take value in $[0,1]$; GB
embedded wind generation capacity, $\text{wcap}_{i}$ (MW) and N2EX day-ahead
electricity price, $\text{n2ex}_{i}$ (€/MWh). Some of the covariates that are
specific to the $j$-region are: net-demand at a 24 hours lag, $y_{ij}^{24}$;
school holidays, $\text{shol}_{ij}$, a factor with three levels to distinguish
Christmas from other holidays, and day of the week with additional factor
levels for public holidays (e.g., boxing day), $\text{dow}^{+}_{ij}$. The last
two covariates are region-specific because some holidays do not coincide
between England, Wales and Scotland. Further, region-specific, covariates are
derived from the day-ahead weather forecasts from the operational ECMWF-HRES
model. The model produces predictions on a regular mesh and a GSP region
covers more than one grid element. For each region, we summarise the gridded
output via the following features: mean of forecasted solar irradiance
(W/$\text{m}^{2}$) times embedded solar generation capacity (MW),
$\text{irr}_{ij}$; forecasted temperature at the cell corresponding to highest
population density in a region, $\text{temp}_{ij}$ (K); 48 hours rolling mean
of $\text{temp}_{ij}$, $\text{temp}_{ij}^{S}$; mean of forecasted wind speed
at 10 and 100 meters, $\text{wsp}^{10}_{ij}$ and $\text{wsp}^{100}_{ij}$
(m/s), and mean of forecasted precipitation, $\text{rain}_{ij}$ (mm/h).

Browell and Fasiolo (2021) use the data just described to model the
conditional distribution of $y_{ij}$, separately for each of the $d=14$
regions. They do so by proposing a composite modelling approach, where the raw
residuals of a Gaussian GAM are modelled via linear quantile regression to
obtain non-parametric quantile estimates, as done by Gaillard et al. (2016).
Extreme conditional quantiles are modelled using a GAMLSS model based on the
generalized Pareto distribution. In this work, we are interested in modelling
the joint distribution of the $d$-dimensional response vector
$\mathbb{y}_{i}$. We consider a multivariate Gaussian model
$\mathbb{y}_{i}\sim N({\boldsymbol{\mu}}_{i},{\boldsymbol{\Sigma}}_{i})$,
where ${\boldsymbol{\mu}}_{i}$ and ${\boldsymbol{\Sigma}}_{i}$ are controlled
by the linear predictors vector ${\boldsymbol{\eta}}_{i}$, as described in
Section 2.1, and each element of ${\boldsymbol{\eta}}_{i}$ is modelled via
(1).

The proposed model has $q=119$ linear predictors and each of them could be
modelled via any of the covariates described above. Hence, model selection is
challenging. As will be explained in the next section, we use data from 2014
to 2016 to generate a large list of candidate covariate effects, ordered in
terms of decreasing importance. We choose the number of effects to add to the
final multivariate Gaussian model, that is where to stop along the ordered
effect list, by maximising the out-of-sample predictive performance on 2017
net-demand. Having chosen the model structure, in Section 3.4 we evaluate the
accuracy of the resulting forecasts on 2018 data.

###  3.2 Semi-Automatic Model Selection

Browell and Fasiolo (2021) consider a progression of univariate GAMs based on
an increasingly rich set of covariates and assess their performance on a day-
ahead forecasting task. We use their results to choose a model for the first
$d=14$ elements of ${\boldsymbol{\eta}}_{i}$ or, equivalently,
${\boldsymbol{\mu}}_{i}$. In particular, we adopt the model formula

$\displaystyle\eta_{ij}$
$\displaystyle=g_{j1}(t_{i})+g_{j2}(t_{i}^{2})+g_{j3}(\text{dow}^{+}_{ij})+g_{j4}(\text{shol}_{ij})+g_{j5}(y_{ij}^{24})+g_{j6}(\text{wsp}_{ij}^{10})$
$\displaystyle+f_{j1}^{20}(\text{doy}_{i})+f_{j2}^{35}(\text{tod}_{i})+f_{j3}^{10}(\text{n2ex}_{i})+f_{j4}^{35}(\text{temp}_{ij})+f_{j5}^{35}(\text{temp}_{ij}^{S})+f_{j6}^{10}(\text{rain}_{ij})$
$\displaystyle+\text{wcap}_{i}\times
f_{j7}^{20}(\text{wsp}_{ij}^{100})+f_{j8}^{5}(\text{irr}_{ij})+f_{j9}^{30}(\text{tod}_{i},\text{dow}^{+}_{ij})+f_{j10}^{20}(\text{tod}_{i},\text{shol}_{ij})$
$\displaystyle+f_{j11}^{5,5}(\text{n2ex}_{i},\text{tod}_{i})+f_{j12}^{5,5}(\text{temp}_{ij},\text{tod}_{i})+f_{j13}^{5,5}(\text{rain}_{ij},\text{tod}_{i})\,\,,$
(9)

for $j=1,\dots,d$. Here $g_{j1}$ to $g_{j6}$ are parametric (linear) effects,
while $f_{j1}$ to $f_{j13}$ are smooth effects. In particular, $f_{j1}$ to
$f_{j8}$ are univariate smooth effects, the spline bases dimensions being
indicated by the superscripts. Effects $f_{j9}$ and $f_{j10}$ are smooth-
factor interactions, where a different univariate smooth is defined for each
level of the $\text{dow}^{+}_{ij}$ or $\text{shol}_{ij}$ factor variables. The
last three effects in (3.2) are bivariate tensor-product smooths, where the
dimension of each marginal basis is indicated by the superscripts. All smooth
effects are built using cubic regression spline bases, with the exception of
$f_{j1}^{20}(\text{doy}_{i})$, which uses a B-spline basis with an adaptive
P-spline penalty. The latter allows the smoothness of the effect to vary with
$\text{doy}_{i}$, see Section 5.3.5 of Wood (2017) for details. Model (3.2) is
almost identical to the model called “GAM-point” in Browell and Fasiolo (2021)
but they use a parametric effect, based on basic trigonometric functions of
$\text{doy}_{i}$, to model the annual seasonality. The more flexible approach
proposed here is meant to model seasonality more flexibly, which is
particularly important around year-end and in densely populated regions, such
as London (see the results in Section 3.3).

It is challenging to develop a model for the remaining 105 elements of
${\boldsymbol{\eta}}_{i}$. As explained in Section 2.1, the elements of the
$\mathbb{T}_{i}$ and $\mathbb{D}_{i}$ matrices correspond to the parameters of
a set of linear models, where the $j$-th element of $\mathbb{y}_{i}$ is
regressed on its predecessors $y_{ij-1}$, …, $y_{i1}$. Hence, the parameters
of the decomposition depend on the ordering of the elements of
$\mathbb{y}_{i}$. For the GSP net-demand data, a sensible ordering can be
chosen on the basis of the location of the GSP regions. As Figure 1 shows, we
order the regions North to South hence $y_{i1}$ and $y_{i14}$ are,
respectively, net-demand in the North of Scotland and in the South West of
England. Under such an ordering, neighbouring regions, which are more likely
to be affected by similar weather and socio-economical events, occupy nearby
positions in $\mathbb{y}_{i}$. Of course, variations on the proposed ordering
could be considered, for example one could think about swapping the order of
regions 12 and 13, which are at a similar latitude.

Given that the complexity of the model considered here would make the search
for some ‘optimal’ ordering burdensome, we design a semi-automatic model
section procedure that does not explicitly take the ordering into account. In
fact, while Pourahmadi (1999) proposed highly structured models for
$\mathbb{D}_{i}$ and $\mathbb{T}_{i}$ which rely on the interpretation of
their elements, and thus on their ordering, we use gradient boosting to chose
which matrix elements should be modelled, and the effects that should be used
to do so. The proposed approach is related to the method of Strömer et al.
(2022b) who use non-cyclical component-wise gradient boosting (Thomas et al.,
2018) to determine the importance of the effects, and then run a further
boosting procedure based on a subset of selected effects, chosen via on a
user-defined importance threshold. In contrast, here we use the out-of-sample
predictive performance to determine the number of effects to include in the
final model and we fit the latter using the methods from Section 2.2, rather
than boosting.

In the following we provide the details of the proposed model selection
approach, part of which is summarised by Algorithm 1. For $j=1,\dots,d$, we
fit a univariate Gaussian GAM $y_{ij}\sim N(\mu_{ij},\sigma_{j}^{2})$, using
net-demand data from the 1st of January 2014 to the 31st of December 2016,
with $\mu_{ij}=\eta_{ij}$ modelled via (3.2). Then, define
$\epsilon_{ij}=y_{ij}-\eta_{ij}$ and let $\tilde{{\boldsymbol{\Sigma}}}$ be
the empirical covariance matrix of such residuals. Initialise each element of
${\boldsymbol{\eta}}_{j}$, for $j=d+1,\dots,q$, using the corresponding
element of the MCD of $\tilde{{\boldsymbol{\Sigma}}}^{-1}$. Let
$\mathcal{R}_{j}$ be a vector containing the indices of the candidate effects
for modelling ${\boldsymbol{\eta}}_{j}$, for $j=d+1,\dots,q$.

The effects in $\mathcal{R}_{j}$ are those appearing in

$\displaystyle\eta_{ij}$
$\displaystyle=g_{j1}(t_{i})+g_{j2}(t_{i}^{2})+g_{j3}(\text{dow}_{i})+f_{j1}^{10}(\text{doy}_{i})+f_{j2}^{10}(\text{tod}_{i})+$
$\displaystyle+f_{j3}^{5}(\text{wsp}_{il_{j}}^{100})+f_{j4}^{5}(\text{irr}_{il_{j}})+f_{j5}^{5}(\text{temp}_{il_{j}})+f_{j6}^{5}(\text{rain}_{il_{j}})+f_{j7}^{5}(\text{n2ex}_{i})\,\,.$
(10)

We assume that the first element of each $\mathcal{R}_{j}$ refers to both
$g_{j1}$ and $g_{j2}$, so that these two terms effectively form a single
effect in the model selection process. For $j=d+1,\dots,q$, we indicate with
$l_{j}$ the row of $\mathbb{D}_{i}$ or $\mathbb{T}_{i}$ on which the $j$-th
linear predictor appears, and with $\text{wsp}_{il_{j}}^{100}$,
$\text{irr}_{il_{j}}$ and so on the weather forecasts corresponding to the
$l_{j}$-th region. Hence, the non-trivial elements of the $l_{j}$-th row of
$\mathbb{D}_{i}$ and $\mathbb{T}_{i}$ are modelled via the weather forecasts
for the $l_{j}$-th region. That is, the forecasts for the North of Scotland
are used to model $(\mathbb{D}_{i})_{11}$ but none of the elements of
$\mathbb{T}_{i}$, those for the South of Scotland are used to model
$(\mathbb{D}_{i})_{22}$ and second row of $\mathbb{T}_{i}$, and so on. To see
the reasoning behind this choice, recall from Section 2.1 that the linear
predictors appearing on the $l_{j}$-row of $\mathbb{T}_{i}$ or
$\mathbb{D}_{i}$ are related to, respectively, the coefficients or the
residual variance of the regression of the $l_{j}$-th element of
$\boldsymbol{y}$ on its predecessors. Hence, it seems reasonable to use the
weather forecasts for the $l_{j}$-th region to model the effect of the
preceding regions on $l_{j}$.

Note that $\mathcal{R}_{j}$ contains only a subset of the effects appearing in
(3.2). In particular, no bivariate tensor-product smooth effect is used to
model the linear predictors controlling the covariance matrix. This choice is
motivated by the fact that we are performing model selection across 105 linear
predictors, so it is important to limit the number of candidate effects to
ensure computational feasibility and statistical parsimony. For the same
reason, the number of basis functions used to construct the effects in (3.2)
is kept low.

Algorithm 1 Quantifying the effects’ importance via gradient boosting

1: Let ${\boldsymbol{\eta}}_{j}$ be the $j$-th linear predictor and indicate
with ${\boldsymbol{\eta}}$ the $n\times q$ matrix containing all the linear
predictors. Initialise ${\boldsymbol{\eta}}_{j}$, for $j=1,\dots,q$, as
explained in Section 3.2. Let $\mathcal{R}_{j}$ be the indices of the
candidate effects for modelling ${\boldsymbol{\eta}}_{j}$ and let
$\mathbb{X}^{r}$, with $r\in\mathcal{R}_{j}$, be the model matrix of the
$r$-th effect. Let $\boldsymbol{\Delta}$ be a list of length $q$, its $j$-th
element $\boldsymbol{\Delta}_{j}$ being a vector of dimension
$\text{card}(\mathcal{R}_{j})$ with all elements, $\Delta_{rj}$, initialised
at 0.

2: For $j=d+1,\dots,q,$

  1. I.

Compute the gradient $\mathbb{u}_{j}$ of the log-likelihood w.r.t.
${\boldsymbol{\eta}}_{j}$, that is

$u_{ij}=\frac{\partial\log
p(\mathbb{y}_{i}|{\boldsymbol{\eta}}_{i})}{\partial\eta_{ij}}\,\,,\quad\text{for}\quad
i=1,\dots,n,$

where ${\boldsymbol{\eta}}_{i}$ is the $i$-th row of ${\boldsymbol{\eta}}$.

  2. II.

For $r\in\mathcal{R}_{j}$

    1. (a)

Project $\mathbb{u}_{j}$ on $\mathbb{X}^{r}$ via
$\hat{\mathbb{u}}_{j}^{r}=\mathbb{X}^{r}({\mathbb{X}^{r}}^{\top}\mathbb{X}^{r})^{-1}{{\mathbb{X}}^{r}}^{\top}\mathbb{u}_{j}$.

    2. (b)

Update the corresponding linear predictor via
$\tilde{{\boldsymbol{\eta}}}_{j}^{r}={\boldsymbol{\eta}}_{j}+\nu\hat{\mathbb{u}}_{j}^{r}$,
where $0<\nu<1$ is the learning rate. Let $\tilde{{\boldsymbol{\eta}}}^{rj}$
be the same as ${\boldsymbol{\eta}}$ matrix, but with the $j$-th column set to
$\tilde{{\boldsymbol{\eta}}}_{j}^{r}$.

    3. (c)

Compute the corresponding change in log-likelihood

$\delta_{rj}=\sum_{i=1}^{n}\left\\{\log
p(\mathbb{y}_{i}|\tilde{{\boldsymbol{\eta}}}_{i}^{rj})-\log
p(\mathbb{y}_{i}|{\boldsymbol{\eta}}_{i})\right\\}.$

  3. III.

Let $j^{*}$ and $r^{*}$ be the indices corresponding to the largest
$\delta_{jr}$. Update the relevant linear predictor and the cumulative gain
vector by doing

${\boldsymbol{\eta}}_{j}\leftarrow\tilde{{\boldsymbol{\eta}}}_{j^{*}}^{r^{*}},\quad\text{and}\quad\Delta_{r^{*}j^{*}}\leftarrow\Delta_{r^{*}j^{*}}+\delta_{r^{*}j^{*}}\,\,.$

3: Unless the maximum number of iterations $M$ has been reached, go back to
step 2.

Having fixed ${\boldsymbol{\eta}}_{j}$, for $j=1,\dots,d$, and initialised
${\boldsymbol{\eta}}_{j}$, for $j=d+1,\dots,q$, to a constant as explained
above, Algorithm 1 fits by least squares the gradient of the log-likelihood
w.r.t. ${\boldsymbol{\eta}}_{j}$, for $j=d+1,\dots,q$, with each candidate
effect in $\mathcal{R}_{j}$. The effect-linear predictor pair that leads to
the largest log-likelihood increase is used to update the model. The
corresponding increase in log-likelihood is stored, and the output of the
procedure is a list of vectors $\boldsymbol{\Delta}_{j}$, for $j=1,\dots,q$.
The $r$-th element of $\boldsymbol{\Delta}_{j}$ contains the cumulative log-
likelihood gain achieved by the $r$-th effect on the $j$-th linear predictor.
Note that $\boldsymbol{\Delta}_{j}$, for $j=1,\dots,d$, are empty vectors
because the model for the first $d$ linear predictors is given by (3.2), hence
Step 2 of Algorithm 1 skips these elements. We select the model for
${\boldsymbol{\eta}}_{j}$, with $j=1,\dots,q$, by running the Algorithm 1 for
$M$ iterations. We sort all the elements of the cumulative gain list,
$\boldsymbol{\Delta}$, from the largest to the smallest to obtain a single
vector. Assuming that the priority with which an effect-linear predictor pair
should be added to the final model is proportional to its cumulative
likelihood gain, the $L$ pairs corresponding to the $L$-largest gains should
be included in the final model, for some $L\geq 0$.

Let $L_{1}=0,L_{2}=5,L_{3}=10,\dots$, be a grid of potential values for the
total number of effects, $L$, used to model ${\boldsymbol{\eta}}_{j}$, for
$j=d+1,\dots,q$. To determine $L$, we optimise the predictive performance of
the full multivariate Gaussian model on net-demand data from year 2017. In
particular, we start by using the univariate Gaussian GAMs fitted to the first
three years of data to predict the value of $\eta_{ij}$, for $j=1,\dots,d$,
with $i$ covering January 2017. Then, we refit the models using data from the
1st of January 2014 to the 31st of January 2017, and obtain predictions for
February 2017. By iterating this rolling forecasting origin procedure twelve
times, we obtain predicted values for $\eta_{ij}$, for $j=1,\dots,d$, covering
the whole of 2017. For each candidate value of $L_{j}$ we fit, using the
rolling procedure described above, the multivariate Gaussian model
$\mathbb{y}_{i}\sim N({\boldsymbol{\mu}}_{i},{\boldsymbol{\Sigma}}_{i})$ on
2017 net-demand data and obtain out-of-sample predictions for the remaining
linear predictors (i.e., those used to model ${\boldsymbol{\Sigma}}_{i}$).
Model fitting is performed via the methods from Section 2.2 and, for each
candidate $L_{j}$, the day-ahead predictions for ${\boldsymbol{\mu}}_{i}$ and
${\boldsymbol{\Sigma}}_{i}$ are used to compute the out-of-sample log-
likelihood.

Here we described a model selection approach which builds on Browell and
Fasiolo (2021) for modelling the mean vector and, given a pool of candidate
effects, uses gradient boosting and a validation set to choose the effects for
modelling the covariance matrix. In the next section we discuss the results
obtained via this model selection procedure.

Figure 3: Model selection results. The diagonal corresponds to the elements of
$\mathbb{D}_{i}$, the rest to those of $\mathbb{T}_{i}$. The symbols represent
different effects and their colour is determined by the magnitude of the
element of $\boldsymbol{\Delta}$ corresponding to each effect.

###  3.3 Model Selection Results

The output of Algorithm 1 is a list $\boldsymbol{\Delta}$ containing the log-
likelihood gains, accumulated up to step $M$, achieved by each effect-linear
predictor pair. As $M$ grows, the linear predictors fit the 2014-16 training
data increasingly well, which eventually leads to over-fitting. To stop the
gains accumulation before over-fitting, we choose $M$ to maximise the log-
likelihood on a validation set comprising data from 2017. By running Algorithm
1 up to $M=2000$, we find that $M=1000$ approximately maximises the out-of-
sample log-likelihood. Running the procedure in parallel on a 12-core
workstation takes several hours.

As explained in Section 3.2, we use the predictive performance on 2017 data to
choose the number of effects, $L$, to be included in the model. The procedure
takes around a day to run on a 12-core workstation and suggests including
$L=60$ effects, when the out-of-sample log-likelihood is used as performance
metric. Figure 3 shows the effects selected to model each element of
$\mathbb{D}_{i}$ and $\mathbb{T}_{i}$, while Figure 4 shows a set of ALEs
obtained by fitting a model containing the effects shown in Figure 3 to all
the data available (2014-2018).

First, note that most of the effects in Figure 3 act on the diagonal elements
of $\mathbb{D}_{i}$ which, under the interpretation detailed in Section 2.1,
are related to the residual variance of the regression of $y_{il}$ on
$y_{il-1}$, $\dots$ $y_{i1}$. This suggests that the effects acting on
$\mathbb{D}_{i}$ are modelling the residual variance of the net-demand in each
region, after having conditioned on the preceding regions. With the exception
of the North of Scotland (Region 1, henceforth R. 1), all the non-zero
elements of $\mathbb{D}_{i}$ depend on the time of day, $\text{tod}_{i}$. It
is not surprising to see that the cumulative log-likelihood gain of the
$\text{tod}_{i}$ effect is particular large in highly urbanised areas, such as
the Midlands (R. 8 and 9) and London (R. 11). The red ALE in Figure 4a shows
the effect of daily seasonality in London, which is characterised by high net-
demand variance during peak hours. The same effect has a similar, but flatter,
shape in the South of Scotland (R. 2). In the South of England (R. 13) the
effect has a single peak and it is even stronger than in London. As we discuss
later, this is likely related to the high capacity of embedded solar
generation relative to electricity consumption.

The time of day is used to model several elements of $\mathbb{T}_{i}$ as well.
The strongest such effect, in terms of cumulative log-likelihood gain, acts on
the element corresponding to London and the West Midlands (R. 11 and 8). The
ALE of $\text{tod}_{i}$ on the correlation between these regions is shown in
green in Figure 4e. Together with Figure 4a, it shows that prediction errors
are larger and more correlated at peak times than during the night or in the
middle of the day in these highly urbanised regions. The second-strongest
effect of $\text{tod}_{i}$ on an element of $\mathbb{T}_{i}$ controls the
dependency between the West Midlands (R. 8) and the South West of England (R.
14). Interestingly, the blue curve in Figure 4e shows that the corresponding
correlation is not maximal at peak times, but at midnight.

Figure 4: Top row: ALEs of the time of day (a), day of year (b), forecasted
wind speed (c) and solar irradiation (d) on the standard deviation of net-
demand in a selected group of regions. Bottom row: ALEs of same covariates on
a selected group of pairwise correlations.

The effect of the day of the year, $\text{doy}_{i}$, is used to model many
elements of $\mathbb{D}_{i}$ and $\mathbb{T}_{i}$. It is not surprising to see
this effect appearing on the 8th and 11th row of Figure 3, which correspond to
the highly urbanised West Midlands (R. 8) and London (R. 11). As the green ALE
in Figure 4b shows, net-demand forecast uncertainty is very high in London at
the end of the year, due to holidays that have a sizeable, hard-to-model
effect on demand patterns. Furthermore, as the effects in Figure 4f show,
uncertainty between regions is also highly correlated during this period,
meaning that forecast errors are likely to have the same sign across regions
as they are driven by the same underlying behavioural effects.

In accordance with intuition, wind speed, $\text{wsp}_{ij}^{100}$, is selected
to model the elements of $\mathbb{D}_{i}$ corresponding to regions with a high
penetration of embedded wind generation, such as the South East of England (R.
12), and the North (R. 1) and South (R. 2) of Scotland. The ALEs of
$\text{wsp}_{ij}^{100}$ in these regions are shown in Figure 4c and could be
interpreted as follows. At low wind, the variability of wind production is low
because little or no generation is occurring, but it increases at modest wind
speeds that are sufficient for power generation to occur, while being in a
range where generation is highly sensitive to wind speed. Then variability
decreases for high wind speeds, where power production is less variable as
turbines self-regulate to maintain their maximum power production. At very
high wind speeds, wind turbines may automatically shut down, and small
differences in wind speed may result in large differences in production,
leading to greater forecast uncertainty.

It is perhaps surprising that wind speed is not selected to control the
element of $\mathbb{T}_{i}$ controlling the dependency between the two
Scottish GSPs, which both contain large amounts of embedded wind. However,
capacity as a fraction of peak load is considerably higher in North than in
South of Scotland. Further, the fact that $(\mathbb{T}_{i})_{21}$ is constant,
does not mean that the correlation between the two Scottish regions is
constant, as illustrated by the green curve in Figure 4g. The plot shows that
the correlation is proportional to the variance in these regions, hence wind
speed controls both the size and correlation between prediction errors.
Interestingly, the blue curve in Figure 4g shows that the net-demand in the
South East of England is less correlated with that in London (R. 11) as wind
speed in the former region increases, which suggests that this covariate
effects the net-demand patterns of these two regions quite differently.

The time of day and solar irradiance, $\text{irr}_{ij}$, are both strongly
related to solar energy production, hence it is interesting to see that the
effects of both variables are selected to model the elements of
$\mathbb{D}_{i}$ corresponding to several Southern regions, which have high
embedded solar generation capacity. The ALEs of $\text{irr}_{ij}$ on the net-
demand variability in South Wales (R. 7), South (R. 13) and South-West (R. 14)
England are shown in Figure 4d. Note that the horizontal scales are different
because the installed solar capacity differs between regions. The shape of
these effects is similar and could be interpreted as follows. Variability is
low at low or high levels of irradiance which correspond to, respectively,
heavily clouded (or night) and to clear sky conditions. Variability is highest
at intermediate levels of irradiance, which might correspond to partial or
broken cloud conditions. However, the shape of the effects may also be
affected by the correlation between irradiance and temperature and by changes
in installed solar capacity over the study period, hence is it important not
to over-interpret them. Solar irradiance is selected to control the dependency
between the South of Wales and the South-West of England via the corresponding
element of $\mathbb{T}_{i}$. This is interesting because, while the two
regions are separated by the Bristol channel, they are geographically close,
hence likely to be affected by similar weather patterns, and they both feature
very high solar penetration relative to peak load.

As explained above, the set of effects selected by the semi-automated
procedure proposed in Section 3.2 matches intuition in many respects. However,
looking again at Figure 3, note that more than half of the selected effects
are related to calendars variables, namely progressive time, time of day and
day of year. Further, some of the candidate variables appearing in (3.2), such
as the day of the week and the N2EX day-ahead price, have not been selected to
model any element of $\mathbb{T}_{i}$ or $\mathbb{D}_{i}$. Hence, it is
interesting to analyse how the predictive performance of the model depends on
the set of candidate covariates that are considered by the model selection
procedure. This is the focus of Section 3.4.

###  3.4 Model Validation

Here we assess the predictive performance of several alternative models
obtained via the selection procedure proposed in Section 3.2. We consider
three sets of candidate effects and two benchmarks with no covariates. The
first set (Full) includes all the effects appearing in (3.2), hence it
corresponds to the model analysed in Section 3.3. Then we consider a calendar-
only model (Cal), obtained by using only the first five effects in (3.2), that
is from $g_{j1}(t_{i})$ to $f_{j2}^{10}(\text{tod}_{i})$, and a larger model
(Cal+Ren) which also includes the effects of wind speed and solar irradiance.
Finally, we consider two baseline models, Static, where the $\mathbb{D}$ and
$\mathbb{T}$ matrices are static, that is they do not depend on any covariate,
and Indep where we use a static, diagonal covariance matrix.

For each set of candidate effects we use data from 2014 to 2017 to perform
model selection, as described in Section 3.2 and 3.3. Having selected the
model structure, we assess performance on 2018 data. In particular, we start
by fitting each model to data up to the end of 2017 and we forecast net-demand
during January 2018. Then, we refit the model using data up to the 31st of
January 2018 and we forecast net-demand for February. By iterating this
rolling forecasting origin procedure, we obtain day-ahead predictions covering
the whole of 2018, for each model. To speed up computation, we first fit model
(3.2) to the net-demand from each region using separate univariate Gaussian
GAMs and then we fit the corresponding residuals vectors using each of the
covariance matrix models described above, using the methods described in
Section 2.2. For each of the three dynamic covariance matrix models, the
procedure takes around a day to run on a 12-core workstation.

We use the day-ahead multivariate predictions obtained as just described to
compute the performance metrics reported in the first four columns of Table 1.
We consider the log score (i.e., the negative log-likelihood, also know as the
ignorance score), the multivariate energy score (Gneiting and Raftery, 2007)
and the $p$-variogram score (Var $p$, Scheuerer and Hamill, 2015) with $p=0.5$
and $p=1$. The results in Table 1 show that enabling the covariance matrix to
vary with the covariates leads to a substantial improvement in performance.
Remarkably, improvement is of the same magnitude as the difference between the
two static models, where the simplest of the two assumes independence between
regions.

| GSP regions |  | GSP macro-regions
---|---|---|---
| Log | Energy | Var 0.5 | Var 1.0 |  | Log | Energy | Var 0.5 | Var 1.0
Indep | 27496 | 3547 | 43403 | 175076 |  | 15534 | 4706 | 9849 | 88731
Static | 24214 | 3519 | 39928 | 161166 |  | 13461 | 4596 | 9609 | 85902
Cal | 23139 | 3475 | 37411 | 150253 |  | 12975 | 4523 | 8816 | 78183
Cal+Ren | 22787 | 3470 | 36772 | 147251 |  | 12882 | 4517 | 8737 | 77353
Full | 22773 | 3470 | 36693 | 146843 |  | 12882 | 4518 | 8734 | 77365
Table 1: Day-ahead performance scores on 2018 test data, when forecasting the
joint distribution of net-demand across the 14 GSP groups (left) and across
the five macro-regions (right) shown in Figure 1, for each of the three sets
of candidate effects and two benchmark Indep and Static models. The best score
in each column has been underlined.

Note that the Cal model is much better that the Static model under all scores,
demonstrating that regional net-demand covariance is characterised by strong
daily and yearly seasonalities. Predictive performance is further improved by
including covariates that are strongly related to embedded renewable
generation, while the gains obtained by further extending the set of candidate
effects are modest in comparison. The performance differences are least
apparent under the energy score, which is known to be insensitive to
differences in dependency structures (Scheuerer and Hamill, 2015).

We further analyse the relative performance of the models by computing
differences in skill scores
($\text{skill}_{\text{model}}=1-\text{score}_{\text{model}}/\text{score}_{\text{baseline}}$)
for four pairs of model. We calculate such differences to quantify the gains
from adding complexity, beginning with Static relative to Indep, followed by
Cal relative to Static, and so on. The skill scores are computed under the log
score and using the Indep model as the baseline. Non-parametric bootstrapping
with week-long blocks is used to quantify the variability of differences in
skill scores between the four pairs of models. The resulting boxplots are
shown in Figure 5a. They show that modelling inter-regional dependencies and
then expanding the set of candidate effects improves performance at each
stage. However, the gains are diminishing and, while the benefits of capturing
calendar and renewable-generation-related effects on the covariance matrix are
substantial, extending the candidate effects pool to include temperature,
rainfall and N2EX electricity price leads to improvements of marginal
statistical significance.

While Figure 5a summarises the improvements in forecasting skill scores over
the whole test set (year 2018), the curves in Figure 5b-c have been obtained
by smoothing the differences in skill scores along the day of the year and the
time of day. They show that the Static, Cal and Cal+Ren models perform
consistently better than their predecessors in the model sequence.
Interestingly, the performance gain of the Cal model relative to Static varies
widely with the two covariates. In particular, it is surprising to see that,
in Figure 5c, the advantage of moving from Static to Cal is minimal at peak
times and maximal at night, the latter being the time when net-demand is least
variable in many regions (see Figure 4a). A tentative explanation is that,
while the time at which net-demand variability is maximal is different
depending on the regions’ characteristics, net-demand variability tends to be
minimal at night across all regions. Hence, the improvements brought by
capturing the effect of the time of day on the covariance matrix are most
sizeable at night.

Figure 5: Left: bootstrapped differences in skill scores between four pairs of
models, computed under the log score. Center and right: smoothed differences
in skill scores, under the same loss and model pairs, as a function of the day
of the year and time of day.

The results discussed so far demonstrate the importance of taking into account
the covariate-dependent nature of the covariance matrix when modelling the
joint distribution of regional net-demand. However, this work is motivated by
the need for spatially coherent probabilistic net-demand forecasts in power
flow studies and, as explained in Section 1, the transmission grid boundaries
of interest in such studies vary depending on, for example, the status of the
network. Hence, it would interesting to verify whether the performance gains
discussed above still hold when the forecast is post-processed to match the
needs of an operationally relevant scenario. While considering realistic
scenarios would require covering engineering aspects that are well beyond the
scope of this work, in Section 1 we proposed aggregating the GSP regions into
the five macro-regions shown in Figure 1, which were motivated by some of the
boundaries of interest described in the 2021 National Grid’s ten year report
(National Grid, 2021).

Joint macro-regional net-demand forecasts are easy to obtain, because they are
linear transformations of the regional forecasts. The last four columns of
Table 1 show the performance of each model when forecasting the joint
distribution of macro-regional net-demand. The results are similar to those
for regional net-demand, but on a different scale due to the aggregation of
net-demand (5 macro-regions vs 14 regions). As in the regional case, the Indep
model leads to the poorest performance, while Cal+Ren and Full produce the
best forecasts, delivering similar improvements to that of Static over Indep.

It is also possible to linearly transform the joint regional forecasts to
obtain marginal probabilistic forecasts of differences in net-demand between
regions or macro-regions. Such forecasts could be of particular interest in
the context of power flow analyses focused on specific transmission grid
boundaries. Here we consider the difference in net-demand between the Scottish
or South macro-regions and the rest of the country, or between London and its
neighbouring regions. These boundaries are of particular interest to the
network operator due to the strong influence of wind/solar generation on power
flows and to constraints related to network capacity and stability. Table 2
reports the performance of each model, when forecasting the marginal
distribution of each net-demand difference.

Given that providing probabilistic forecasts of net-demand differences between
(aggregations of) regions requires accurate quantification of inter-
dependencies, it is no surprise that Indep performs relatively poorly in the
Scotland- and South-minus-rest scenarios. In the London-minus-neighbours
scenario, the number of regions over which net-demand is aggregated is lower
than in the other two scenarios, which might explain why Indep performs
similarly to Static in this case. Note that, while the Indep and Static models
are still the worst, the performance gains obtained by explicitly modelling
the covariate-dependent nature of net-demand covariance matrix are less
substantial, relative to Table 1. Arguably, assessing the quality of a joint
net-demand forecast should be more difficult when the performance is evaluated
on a one-dimensional summary of the net-demand vector, rather than on the
whole vector. This might also explain why the Full model does not seem to
offer any benefit here but, considering the very moderate benefits relative to
Cal+Ren reported in Table 1, one might argue that this model might be slightly
over-fit. Nevertheless, considering the need for probabilistic forecasts that
can be flexibly aggregated (or more generally transformed) to match the needs
of specific power flow (or other) analyses, it is a notable result that all
the covariate-dependent matrix models considered here outperform the Static
model under all the forecast transformations and performance scores considered
in Tables 1 and 2. Due to the high cost of operating power systems (and volume
of energy traded in wholesale markets), marginal improvement in forecast
performance and associated decision-making can yield substantial economic and
operational benefits.

| Sco - Rest | Sou - Rest | Lon - Nei
---|---|---|---
| Log | Energy | Log | Energy | Log | Energy
Indep | 5879 | 6169 | 4495 | 4310 | 2842 | 1491
Static | 4645 | 5790 | 4206 | 4221 | 2850 | 1489
Cal | 4543 | 5701 | 4117 | 4150 | 2715 | 1454
Cal+Ren | 4541 | 5698 | 4121 | 4150 | 2695 | 1450
Full | 4545 | 5703 | 4122 | 4153 | 2703 | 1452
Table 2: Day-ahead performance scores on 2018 test data, when forecasting the
marginal distribution of differences in net-demand defined in Section 3.4. The
best score in each column has been underlined. Note that, in one dimension,
the energy score is equivalent to the continuous ranked probability score
(CRPS).

##  4 Conclusion

Forecasts of supply and demand are essential inputs to predict and manage
power flows on electricity networks, as well as prices and other important
variables. Given the imperative to maintain a reliable electricity supply,
these predictions must enable risk to be quantified and managed. As the
complexity of energy systems increases, the heuristic approaches widely used
today are becoming inadequate and will have to be replaced by explicit
probabilistic forecasts of power flows (Morales et al., 2014).

Motivated by the need for spatially coherent, probabilistic net-demand
forecasts to support energy system operations, we have focused on joint day-
ahead forecasting of net-demand across the GSP regions comprising GB’s
transmission system. To accommodate for the dynamic nature of the net-demand
covariance matrix, we let the elements of its MCD parametrisation vary with a
number of temporal and weather-related covariates. To perform effect selection
for a model comprising more than one hundred linear predictors, we leverage
the interpretability of the chosen parametrisation and we combine it with a
semi-automatic effect selection method, based on gradient-boosting. The
results on the test set show that additive covariance matrix models
significantly outperform, in terms of forecasting performance, a model where
the covariance matrix is static.

A promising direction for future work would be to extend the model presented
here to capture temporal, in addition to spatial, dependencies. In particular,
the covariance matrix models used here implicitly assume that regional net-
demand residual vectors are uncorrelated in time. While the mean vector model
(3.2) contains the effect of lagged net-demand, which is meant to capture part
of the intra-regional temporal dependencies, more complex temporal effects
could be captured by extending the covariance matrix to explicitly model the
longitudinal nature of the data considered here. Such an extension should lead
to models able to generate multivariate net-demand trajectories that are
coherent both in space and in time, thus supporting important operations (e.g.
determining the schedules for power generating units) that consider both
spatial and temporal constraints.

## References

  * Apley and Zhu (2020) Apley, D. W. and J. Zhu (2020).  Visualizing the effects of predictor variables in black box supervised learning models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 82(4), 1059–1086.
  * Bonat and Jørgensen (2016) Bonat, W. H. and B. Jørgensen (2016).  Multivariate covariance generalized linear models.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 65(5), 649–675.
  * Browell and Fasiolo (2021) Browell, J. and M. Fasiolo (2021).  Probabilistic forecasting of regional net-load with conditional extremes and gridded NWP.  IEEE Transactions on Smart Grid 12(6), 5011–5019.
  * Browell et al. (2022) Browell, J., C. Gilbert, and M. Fasiolo (2022).  Covariance structures for high-dimensional energy forecasting.  Electric Power Systems Research 211, 108446.
  * Capezza et al. (2021) Capezza, C., B. Palumbo, Y. Goude, S. N. Wood, and M. Fasiolo (2021).  Additive stacking for disaggregate electricity demand forecasting.  The Annals of Applied Statistics 15(2), 727–746.
  * Fan and Hyndman (2012) Fan, S. and R. J. Hyndman (2012).  Short-term load forecasting based on a semi-parametric additive model.  IEEE Transactions on Power Systems 27(1), 134–141.
  * Friedman (2001) Friedman, J. H. (2001).  Greedy function approximation: A gradient boosting machine.  The Annals of Statistics 29(5), 1189–1232.
  * Gaillard et al. (2016) Gaillard, P., Y. Goude, and R. Nedellec (2016).  Additive models and robust aggregation for GEFCom2014 probabilistic electric load and electricity price forecasting.  International Journal of Forecasting 32(3), 1038–1050.
  * Gneiting and Raftery (2007) Gneiting, T. and A. E. Raftery (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association 102(477), 359–378.
  * Hans et al. (2022) Hans, N., N. Klein, F. Faschingbauer, M. Schneider, and A. Mayr (2022).  Boosting distributional copula regression.  arXiv preprint arXiv:2202.12851v1.
  * Hastie and Tibshirani (1987) Hastie, T. and R. Tibshirani (1987).  Generalized additive models: Some applications.  Journal of the American Statistical Association 82(398), 371–386.
  * Huxley et al. (2022) Huxley, O., J. Taylor, A. Everard, J. Briggs, K. Tilley, J. Harwood, and A. Buckley (2022).  The uncertainties involved in measuring national solar photovoltaic electricity generation.  Renewable and Sustainable Energy Reviews 156, 112000\.
  * Klein et al. (2015) Klein, N., T. Kneib, S. Klasen, and S. Lang (2015).  Bayesian structured additive distributional regression for multivariate responses.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 64(4), 569–591.
  * Marra and Radice (2017) Marra, G. and R. Radice (2017).  Bivariate copula additive models for location, scale and shape.  Computational Statistics & Data Analysis 112, 99–113.
  * Morales et al. (2014) Morales, J. M., A. J. Conejo, H. Madsen, P. Pinson, and M. Zugno (2014).  Integrating Renewables in Electricity Markets: Operational Problems, Volume 205.  New York, NY, USA: Springer.
  * Muschinski et al. (2022) Muschinski, T., J. G. Mayr, T. Simon, N. Umlauf, and A. Zeileis (2022).  Cholesky-based multivariate Gaussian regression.  Econometrics and Statistics (to appear).
  * National Grid (2018) National Grid (2018).  Regional trends and insights.  nationalgrideso.com/sites/eso/files/documents/SOF%20Report-%20Regional%20Trends%20and%20Insights.pdf.
  * National Grid (2021) National Grid (2021).  Electricity ten year statement.  nationalgrideso.com/document/223046/download.
  * Nychka (1988) Nychka, D. (1988).  Bayesian confidence intervals for smoothing splines.  Journal of the American Statistical Association 83(404), 1134–1143.
  * Pinheiro and Bates (1996) Pinheiro, J. C. and D. M. Bates (1996).  Unconstrained parametrizations for variance-covariance matrices.  Statistics and Computing 6(3), 289–296.
  * Pourahmadi (1999) Pourahmadi, M. (1999).  Joint mean-covariance models with applications to longitudinal data: Unconstrained parameterisation.  Biometrika 86(3), 677–690.
  * Pourahmadi (2011) Pourahmadi, M. (2011).  Covariance estimation: The GLM and regularization perspectives.  Statistical Science 26(3), 369–387.
  * Rigby and Stasinopoulos (2005) Rigby, R. A. and D. M. Stasinopoulos (2005).  Generalized additive models for location, scale and shape.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 54(3), 507–554.
  * Scheuerer and Hamill (2015) Scheuerer, M. and T. M. Hamill (2015).  Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities.  Monthly Weather Review 143(4), 1321–1334.
  * Strömer et al. (022a) Strömer, A., N. Klein, C. Staerk, H. Klinkhammer, and A. Mayr (2022a).  Boosting multivariate structured additive distributional regression models.  arXiv preprint arXiv:2207.08470v1.
  * Strömer et al. (022b) Strömer, A., C. Staerk, N. Klein, L. Weinhold, S. Titze, and A. Mayr (2022b).  Deselection of base-learners for statistical boosting—with an application to distributional regression.  Statistical Methods in Medical Research 31(2), 207–224.
  * Thomas et al. (2018) Thomas, J., A. Mayr, B. Bischl, M. Schmid, A. Smith, and B. Hofner (2018).  Gradient boosting for distributional regression: Faster tuning and improved variable selection via noncyclical updates.  Statistics and Computing 28(3), 673–687.
  * Tuinema et al. (2020) Tuinema, B. W., J. L. Rueda Torres, A. I. Stefanov, F. M. Gonzalez-Longatt, and M. A. M. M. van der Meijden (2020).  Probabilistic power flow analysis.  In Probabilistic Reliability Analysis of Power Systems, pp. 179–208. Cham, CH: Springer.
  * Von Meier (2006) Von Meier, A. (2006).  Electric Power Systems: A Conceptual Introduction.  Hoboken, NJ, USA: John Wiley & Sons.
  * Wood (2017) Wood, S. N. (2017).  Generalized Additive Models: An Introduction with R. (2 ed.).  Boca Raton, FL, USA: Chapman & Hall/CRC.
  * Wood and Fasiolo (2017) Wood, S. N. and M. Fasiolo (2017).  A generalized Fellner-Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models.  Biometrics 73(4), 1071–1081.
  * Wood et al. (2016) Wood, S. N., N. Pya, and B. Säfken (2016).  Smoothing parameter and model selection for general smooth models.  Journal of the American Statistical Association 111(516), 1548–1563.

## Supplementary Material to “Additive Covariance Matrix Models: Modelling
Regional Electricity Net-Demand in Great Britain”

##  Appendix A Derivatives of the Log-Likelihood

###  A.1 Setting up the Notation

Consider a scalar-valued function $f$ of the $n$-dimensional vectors
${\boldsymbol{\eta}}_{1},\dots,{\boldsymbol{\eta}}_{q}$. We indicate with
$f^{{\boldsymbol{\eta}}_{k}}$ and
$f^{{\boldsymbol{\eta}}_{k}{\boldsymbol{\eta}}_{j}}$ the vectors with $i$-th
elements

$f^{\eta_{ik}}=\frac{\partial f}{\partial\eta_{ik}}\quad\text{and}\quad
f^{\eta_{ik}\eta_{ij}}=\frac{\partial^{2}f}{\partial\eta_{ik}\partial\eta_{ij}}\,\,,$

where $\eta_{ik}$ indicates the $i$-th element of ${\boldsymbol{\eta}}_{k}$.
Each ${\boldsymbol{\eta}}_{k}$ is a function of a corresponding
$p_{k}$-dimensional vector ${\boldsymbol{\beta}}_{k}$. For the derivatives of
$f$ w.r.t. the elements of ${\boldsymbol{\beta}}_{k}$, we use the compact
notation

$f^{\beta_{kr}}=\frac{\partial f}{\partial\beta_{kr}}\quad\text{and}\quad
f^{\beta_{kr}\beta_{js}}=\frac{\partial^{2}f}{\partial\beta_{kr}\partial\beta_{js}}\,\,,$

where $\beta_{kr}$ indicates the $r$-th element of ${\boldsymbol{\beta}}_{k}$.
Finally, we denote with
$f^{{\boldsymbol{\beta}}_{k}}=\nabla_{{\boldsymbol{\beta}}_{k}}f$ the gradient
of $f$ w.r.t. ${\boldsymbol{\beta}}_{k}$ and with
$f^{{\boldsymbol{\beta}}_{k}{\boldsymbol{\beta}}_{j}}=\nabla_{{\boldsymbol{\beta}}_{j}}^{\top}\nabla_{{\boldsymbol{\beta}}_{k}}f$
the matrix of second derivatives.

###  A.2 Gradient and Hessian w.r.t. ${\boldsymbol{\beta}}$

To simplify the notation, let us indicate with $\mathbb{y}$ the collection of
all response vectors $\mathbb{y}_{1},\dots,\mathbb{y}_{n}$ and define
$\mathcal{L}({\boldsymbol{\beta}})=\log
p({\boldsymbol{\beta}}|\mathbb{y},{\boldsymbol{\lambda}})$. Recall that

$\mathcal{L}({\boldsymbol{\beta}})=\sum_{i=1}^{n}\ell_{i}-\frac{1}{2}{\boldsymbol{\beta}}^{\top}\tilde{\bf
S}^{\boldsymbol{\lambda}}{\boldsymbol{\beta}}$

where $\ell_{i}=\log p(\mathbb{y}_{i}|{\boldsymbol{\beta}})$. The gradient and
Hessian of the log-posterior w.r.t ${\boldsymbol{\beta}}$ are

$\mathcal{L}^{{\boldsymbol{\beta}}}({\boldsymbol{\beta}})=\sum_{i=1}^{n}\ell_{i}^{{\boldsymbol{\beta}}}-\tilde{\mathbf{S}}^{{\boldsymbol{\lambda}}}{\boldsymbol{\beta}}\quad\text{and}\quad\mathcal{L}^{{\boldsymbol{\beta}}{\boldsymbol{\beta}}}({\boldsymbol{\beta}})=\sum_{i=1}^{n}\ell_{i}^{{\boldsymbol{\beta}}{\boldsymbol{\beta}}}-\tilde{\mathbf{S}}^{{\boldsymbol{\lambda}}}\,\,.$

Let us define $\bar{\ell}=\sum_{i=1}^{n}\ell_{i}$ . To provide formulas for
$\bar{\ell}^{{\boldsymbol{\beta}}}$ and
$\bar{\ell}^{{\boldsymbol{\beta}}{\boldsymbol{\beta}}}$, let us assume that
${\boldsymbol{\beta}}=({{\boldsymbol{\beta}}_{1}}^{\top},\dots,{{\boldsymbol{\beta}}_{q}}^{\top})^{\top}$,
where ${\boldsymbol{\beta}}_{j}$ is the vector of regression coefficients
specific to the $j$-th linear predictor, that is
${\boldsymbol{\eta}}_{j}=\mathbf{X}^{j}{\boldsymbol{\beta}}_{j}$ where
$\mathbf{X}^{j}$ is an $n\times p_{j}$ model matrix. With this notation, the
$j$-th sub-vector of $\bar{\ell}^{{\boldsymbol{\beta}}}$ is

$\bar{\ell}^{{\boldsymbol{\beta}}_{j}}={\mathbf{X}^{j}}^{\top}\bar{\ell}^{{\boldsymbol{\eta}}_{j}}\,\,,$

while the $j,k$-th block of the Hessian is

$\bar{\ell}^{{\boldsymbol{\beta}}_{j}{\boldsymbol{\beta}}_{k}}=(\bar{\ell}^{{\boldsymbol{\beta}}_{k}{\boldsymbol{\beta}}_{j}})^{\top}={\mathbf{X}^{k}}^{\top}\text{diag}(\bar{\ell}^{{\boldsymbol{\eta}}_{j}{\boldsymbol{\eta}}_{k}}){\mathbf{X}^{j}}\,\,,$

where ${\rm diag}(\cdot)$ is the vector-to-matrix diagonal operator. The
formulas provided so far are applicable to any GAM with multiple linear
predictors and independent response vectors. In contrast, the expressions for
$\bar{\ell}^{{\boldsymbol{\eta}}_{j}}$ and
$\bar{\ell}^{{\boldsymbol{\eta}}_{j}{\boldsymbol{\eta}}_{k}}$ are model-
specific and are provided in the following section for a multivariate Gaussian
distribution, with covariance matrix parametrised via the MCD.

###  A.3 Derivatives w.r.t. ${\boldsymbol{\eta}}$

Let us start by defining a few useful quantities. Let $\mathbb{G}$ be a
$(d-1)\times(d-1)$ lower triangular matrix such that
$G_{jk}=C_{jk}+2d\mathbbm{1}_{\\{k\leq j\\}}$, where

${\rm C}_{jk}=\begin{cases}\quad\binom{j+1}{2}\quad\quad\quad\,\,\,k=j\\\ {\rm
C}_{j(k+1)}-1\quad\quad\,k<j\\\ \quad\quad
0\quad\quad\quad\quad\,\,k>j\,\,,\end{cases}$

and $\mathbbm{1}$ is the indicator function. Define the $(d-1)\times(d-1)$
lower triangular matrices $\mathbb{Z}$ and $\mathbb{W}$ such that ${\rm
Z}_{jk}=k\mathbbm{1}_{\\{k\leq j\\}}$ and ${\rm
W}_{jk}=(j+1)\mathbbm{1}_{\\{k\leq j\\}}$. Let
$\mathbb{z}=\operatorname{rvech}(\mathbb{Z})$ and
$\mathbb{w}=\operatorname{rvech}(\mathbb{W})$, where
$\operatorname{rvech}(\cdot)$ is the row-wise half-vectorisation operator,
that is
$\operatorname{rvech}(\mathbb{Z})=(Z_{11},Z_{21},Z_{22},Z_{31},Z_{32},Z_{33},\dots,Z_{(d-1)(d-1)})^{\top}$.
Let $\mathbb{Q}_{l}$, for $l=1,\ldots,d$, and $\mathbb{P}_{l}$, for
$l=1,\ldots,d(d-1)/2$, by $d\times d$ matrices such that
$(\mathbb{Q}_{l})_{ll}=e^{-\eta_{l+d}}$ and
$(\mathbb{P}_{l})_{{z}_{l}{w}_{l}}=1$, while all other elements are equal to
zero.

Here index $i$ is not needed, hence we drop it and we indicate the $i$-th log-
likelihood component $\ell_{i}$ simply with $\ell$. Note that, given that we
are focusing on an individual $i$, here ${\boldsymbol{\eta}}$ is a
$q$-dimensional vector and $q=d+d(d+1)/2$. If we omit the constants that do
not depend on ${\boldsymbol{\eta}}$ and we indicate with $r_{k}$ the $k$-th
element of the residual vector, $\mathbb{r}=\mathbb{y}-\boldsymbol{\mu}$, the
Gaussian log-density can be written

$\displaystyle\ell$ $\displaystyle=-\frac{1}{2}\left\\{{\rm
tr}(\log\mathbb{D}^{2})+\mathbb{r}^{\top}\mathbb{T}^{\top}\mathbb{D}^{-2}\mathbb{T}\mathbb{r}\right\\}$
$\displaystyle=-\frac{1}{2}\sum_{j=1}^{d}\Big{\\{}\eta_{j+d}+e^{-\eta_{j+d}}\Big{(}\sum_{k=1}^{j-1}\eta_{{\rm
G}_{(j-1)k}}r_{k}+r_{j}\Big{)}^{2}\Big{\\}}\,\,,$

where we used $\log|\boldsymbol{\Sigma}|={\rm
tr}(\log\mathbb{D}^{2})=\sum_{j=1}^{d}\eta_{j+d}$ and we implicitly assumed
that the sum $\sum_{k=1}^{j-1}$ should not be computed when $j=1$ (we will use
the same convention in several places below). Similarly, below we assume that
$\sum_{j=l+1}^{d}$ will not be computed when $l=d$. Here we provide the first
and second derivatives of $\ell$ w.r.t. ${{\boldsymbol{\eta}}}$ both in
compact matrix form and in an extended format, the latter being more useful
for efficient numerical implementation.

With notation above, the elements of
$\ell^{{\boldsymbol{\eta}}}=(\ell_{1}^{{\boldsymbol{\eta}}},\dots,\ell_{q}^{{\boldsymbol{\eta}}})^{\top}=(\partial\ell/\partial\eta_{1},\dots,\partial\ell/\partial\eta_{q})^{\top}$
are

$\displaystyle\ell^{{\boldsymbol{\eta}}}_{l}$
$\displaystyle=\left(\mathbb{T}^{\top}\mathbb{D}^{-2}\mathbb{T}\mathbb{r}\right)_{l}$
$\displaystyle=e^{-\eta_{d+l}}\Big{(}\sum_{k=1}^{l-1}\eta_{{\rm
G}_{(l-1)k}}r_{k}+r_{l}\Big{)}+\sum_{j=l+1}^{d}e^{-\eta_{j+d}}\Big{(}\sum_{k=1}^{j-1}\eta_{{\rm
G}_{(j-1)k}}r_{k}+r_{j}\Big{)}\eta_{{\rm G}_{(j-1)l}}\,\,,$

for $l=1,\ldots,d$,

$\displaystyle\ell^{{\boldsymbol{\eta}}}_{l}$
$\displaystyle=\frac{1}{2}\mathbb{r}^{\top}\mathbb{T}^{\top}\mathbb{Q}_{l-d}\mathbb{T}\mathbb{r}-\frac{1}{2}$
$\displaystyle=\frac{1}{2}e^{-\eta_{l}}\Big{(}\sum_{k=1}^{l-d-1}\eta_{{\rm
G}_{(l-d-1)k}}r_{k}+r_{l-d}\Big{)}^{2}-\frac{1}{2}\,\,,$

for $l=d+1,\ldots,2d$, and

$\displaystyle\ell^{{\boldsymbol{\eta}}}_{l}$
$\displaystyle=-\mathbb{r}^{\top}\mathbb{P}_{l-2d}\mathbb{D}^{-2}\mathbb{T}\mathbb{r}$
$\displaystyle=-e^{\eta_{w_{l-2d}+d}}\Big{(}\sum_{k=1}^{w_{l-2d}-1}\eta_{{\rm
G}_{(w_{l-2d}-1)k}}r_{k}+r_{w_{l-2d}}\Big{)}r_{z_{l-2d}}\,\,,$

for $l=2d+1,\ldots,q$.

The elements forming the upper triangle of
$\ell^{{\boldsymbol{\eta}}{\boldsymbol{\eta}}}$ (here
$\ell^{{\boldsymbol{\eta}}{\boldsymbol{\eta}}}_{lm}=\partial^{2}\ell/\partial\eta_{l}\partial\eta_{m}$),
are

$\displaystyle\ell^{{\boldsymbol{\eta}}{\boldsymbol{\eta}}}_{lm}$
$\displaystyle=-\left(\mathbb{T}^{\top}\mathbb{D}^{-2}\mathbb{T}\right)_{lm}$
$\displaystyle=-\Big{\\{}e^{-{\eta}_{l+d}}+\sum_{k=l+1}^{d}e^{-{\eta}_{k+d}}\Big{(}{\eta}_{{\rm
G}_{(k-1)l}}\Big{)}^{2}\Big{\\}}\mathbbm{1}_{\\{m=l\\}}$ $\displaystyle\hskip
14.22636pt-\Big{(}e^{-{\eta}_{m+d}}{\eta}_{{\rm
G}_{(m-1)l}}+\sum_{k=m+1}^{d}e^{-{\eta}_{k+d}}{\eta}_{{\rm
G}_{(k-1)l}}{\eta}_{{\rm G}_{(k-1)m}}\Big{)}\mathbbm{1}_{\\{m>l\\}}\,\,,$

for $l=1,\ldots,d$ and $m=l,\ldots,d$,

$\displaystyle\ell^{{\boldsymbol{\eta}}{\boldsymbol{\eta}}}_{lm}$
$\displaystyle=-\left(\mathbb{T}^{\top}\mathbb{Q}_{m-d}\mathbb{T}\mathbb{r}\right)_{l}$
$\displaystyle=-e^{-{\eta}_{m}}\Big{\\{}\Big{(}\sum_{k=1}^{l-1}\eta_{G_{(l-1)k}}r_{k}+r_{l}\Big{)}\mathbbm{1}_{\\{m-d=l\\}}$
$\displaystyle\hskip
14.22636pt+\Big{(}\sum_{k=1}^{m-d-1}\eta_{G_{(m-d-1)k}}r_{k}+r_{m-d}\Big{)}\eta_{G_{(m-d-1)l}}\mathbbm{1}_{\\{m-d>l\\}}\Big{\\}}\,\,,$

for $l=1,\ldots,d$ and $m=d+1,\ldots,2d$,

$\displaystyle\ell^{{\boldsymbol{\eta}}{\boldsymbol{\eta}}}_{lm}$
$\displaystyle=\left(\mathbb{P}_{m-2d}\mathbb{D}^{-2}\mathbb{T}\mathbb{r}+\mathbb{T}^{\top}\mathbb{D}^{-2}\mathbb{P}^{\top}_{m-2d}\mathbb{r}\right)_{l}\,\,$
$\displaystyle=e^{-\eta_{w_{m-2d}+d}}\Big{\\{}r_{z_{m-2d}}\Big{(}\mathbbm{1}_{\\{w_{m-2d}=l\\}}+\eta_{{\rm
G}_{(w_{m-2d}-1)l}}\mathbbm{1}_{\\{w_{m-2d}>l\\}}\Big{)}$ $\displaystyle\hskip
14.22636pt+\Big{(}\sum_{k=1}^{w_{m-2d}-1}\eta_{{\rm
G}_{(w_{m-2d}-1)k}}r_{k}+r_{w_{m-2d}}\Big{)}\mathbbm{1}_{\\{z_{m-2d}=l\\}}\Big{\\}}\,\,,$

for $l=1,\ldots,d$ and $m=2d+1,\ldots,q$,

$\displaystyle\ell^{{\boldsymbol{\eta}}{\boldsymbol{\eta}}}_{lm}$
$\displaystyle=-\frac{1}{2}\mathbb{r}^{\top}\mathbb{T}^{\top}\mathbb{Q}_{l-d}\mathbb{T}\mathbb{r}$
$\displaystyle=-\frac{1}{2}e^{-\eta_{l}}\Big{(}\sum_{k=1}^{l-d-1}\eta_{{\rm
G}_{(l-d-1)k}}r_{k}+r_{l-d}\Big{)}^{2}\mathbbm{1}_{\\{m=l\\}}\,\,,$

for $l=d+1,\ldots,2d$ and $m=l,\ldots,2d$,

$\displaystyle\ell^{{\boldsymbol{\eta}}{\boldsymbol{\eta}}}_{lm}$
$\displaystyle=\mathbb{r}^{\top}\mathbb{P}_{m-2d}\mathbb{Q}_{l-d}\mathbb{T}\mathbb{r}$
$\displaystyle=e^{-\eta_{l}}\Big{(}\sum_{k=1}^{l-d-1}\eta_{{\rm
G}_{(l-d-1)k}}r_{k}+r_{l-d}\Big{)}r_{z_{m-2d}}\mathbbm{1}_{\\{w_{m-2d}=l-d\\}}\,\,,$

for $l=d+1,\ldots,2d$ and $m=2d+1,\ldots,q$, and finally

$\displaystyle\ell^{{\boldsymbol{\eta}}{\boldsymbol{\eta}}}_{lm}$
$\displaystyle=-\mathbb{r}^{\top}\mathbb{P}_{l-2d}\mathbb{D}^{-2}\mathbb{P}^{\top}_{m-2d}\mathbb{r}$
$\displaystyle=-e^{-{\eta}_{w_{l-2d}+d}}r_{z_{l-2d}}r_{z_{m-2d}}\mathbbm{1}_{\\{w_{m-2d}=w_{l-2d}\\}}\,\,,$

for $l=2d+1,\ldots,q$ and $m=l,\ldots,q$.

###  A.4 Jacobian Formulas

Recall that ${{\boldsymbol{\eta}}}_{j}=(\eta_{1j},\dots,\eta_{nj})^{\top}$ is
the vector containing the values of the $j$-th linear predictor at each
observation. For $j=1,\dots,q$, we have that
${\boldsymbol{\eta}}_{j}={\mathbb{X}}^{j}{\boldsymbol{\beta}}_{j}$ where
${\mathbb{X}}^{j}$ and ${\boldsymbol{\beta}}_{j}$ are, respectively, the
$n\times p_{j}$ model matrix and the $p_{j}$-dimensional vector of regression
coefficients belonging to the $j$-th linear predictor. Define the vector
$\boldsymbol{\sigma}_{jk}=\\{({\boldsymbol{\Sigma}}_{1})_{jk},\dots,({\boldsymbol{\Sigma}}_{n})_{jk}\\}^{\top}$.
Then, the $n\times p$ Jacobian matrix of $\boldsymbol{\sigma}_{jk}$ w.r.t.
${\boldsymbol{\beta}}$ is

$\mathbf{J}^{jk}=\nabla_{{\boldsymbol{\beta}}}^{\top}\boldsymbol{\sigma}_{jk}=(\nabla_{{\boldsymbol{\beta}}_{1}}^{\top}\boldsymbol{\sigma}_{jk},\cdots,\nabla_{{\boldsymbol{\beta}}_{q}}^{\top}\boldsymbol{\sigma}_{jk})\,\,,$

where $p=\sum_{j=1}^{q}p_{j}$. The $a$-th block of $\mathbf{J}^{jk}$ is

$\nabla_{{\boldsymbol{\beta}}_{a}}^{\top}\boldsymbol{\sigma}_{jk}=\nabla_{{\boldsymbol{\eta}}_{a}}^{\top}\boldsymbol{\sigma}_{jk}\nabla^{\top}_{{\boldsymbol{\beta}}_{a}}{{\boldsymbol{\eta}}_{a}}=\nabla_{{\boldsymbol{\eta}}_{a}}^{\top}\boldsymbol{\sigma}_{jk}{\bf
X}^{a}\,\,,$

for $a,k=1,\dots,q$, and where
$\nabla_{{\boldsymbol{\eta}}_{a}}^{\top}\boldsymbol{\sigma}_{jk}$ is an
$n\times n$ diagonal matrix with non-zero elements

$\left(\nabla_{{\boldsymbol{\eta}}_{a}}^{\top}\boldsymbol{\sigma}_{jk}\right)_{ii}=\frac{\partial({\boldsymbol{\Sigma}}_{i})_{jk}}{\partial\eta_{a}}\,\,.$

Note $\partial({\boldsymbol{\Sigma}}_{i})_{jk}/\partial\eta_{a}$ is the only
parametrisation-dependent component of the Jacobian. In Section A.4.1 we
provide formulas for
$\nabla_{{\boldsymbol{\eta}}_{a}}^{\top}\boldsymbol{\sigma}_{jk}$ under the
MCD parametrisation. The Jacobian of
$\\{(\boldsymbol{\Gamma}_{1})_{jk},\dots,(\boldsymbol{\Gamma}_{n})_{jk}\\}$
w.r.t. ${\boldsymbol{\beta}}$ is computed analougously, but with
$\partial(\boldsymbol{\Gamma}_{i})_{jk}/\partial\eta_{a}$ in place of
$\partial({\boldsymbol{\Sigma}}_{i})_{jk}/\partial\eta_{a}$. Formulas for the
former quantity are provided in Section A.4.2.

####  A.4.1 Derivatives of $\boldsymbol{\Sigma}$ w.r.t. ${\boldsymbol{\eta}}$

Here the index $i$ is not needed, hence we drop it. Consider the factorisation
$\boldsymbol{\Sigma}=\mathbb{R}\mathbb{R}^{\top}$ where
$\mathbb{R}=\mathbb{L}\mathbb{D}$ and $\mathbb{L}=\mathbb{T}^{-1}$. The
partial derivative of the $(l,m)$ element of $\boldsymbol{\Sigma}$ w.r.t.
$\eta_{j}$, is

$\frac{\partial\Sigma_{lm}}{\partial\eta_{j}}=\sum_{k=1}^{d}\left(\frac{\partial{\rm
R}_{lk}}{\partial\eta_{j}}{\rm R}_{mk}+{\rm R}_{lk}\frac{\partial{\rm
R}_{mk}}{\partial\eta_{j}}\right),$

where

$\frac{\partial{\rm R}_{lk}}{\partial\eta_{j}}=0\,\,,\quad\text{for}\quad
j=1,\ldots,d,$ $\frac{\partial{\rm
R}_{lk}}{\partial\eta_{j}}=\frac{\partial{\rm L_{lk}}{\rm
D_{kk}}}{\partial\eta_{j}}=\frac{1}{2}{\rm L_{l(j-d)}}{\rm
D_{(j-d)(j-d)}}\mathbbm{1}_{\\{j-d=k\\}}\,\,,\quad\text{for}\quad
j=d+1,\ldots,2d,$

and

$\frac{\partial{\rm R}_{lk}}{\partial\eta_{j}}=\frac{\partial{\rm L_{lk}}{\rm
D_{kk}}}{\partial\eta_{j}}=-{\rm L_{ls}}{\rm L_{tk}}{\rm
D_{kk}}\,\,,\quad\text{for}\quad j=2d+1,\ldots,q,$

with $s$ and $t$ being the indices of the only element of $\mathbb{T}$ that
depends on $\eta_{j}$. Hence, we obtain

$\frac{\partial\Sigma_{lm}}{\partial\eta_{j}}=0\,\,,\quad\text{for}\quad
j=1,\ldots,d,$ $\frac{\partial\Sigma_{lm}}{\partial\eta_{j}}={\rm
L_{l(j-d)}}{\rm L_{m(j-d)}}{\rm D_{(j-d)(j-d)}}\,\,,\quad\text{for}\quad
j=d+1,\ldots,2d,$

and

$\frac{\partial\Sigma_{lm}}{\partial\eta_{j}}=-{\rm L_{ls}}\sum_{k=1}^{d}{\rm
L_{tk}}{\rm D^{2}_{kk}}{\rm L_{mk}}-{\rm L_{ms}}\sum_{k=1}^{d}{\rm L_{lk}}{\rm
D^{2}_{kk}}{\rm L_{tk}}\,\,,\quad\text{for}\quad j=2d+1,\ldots,q,$

where $s$ and $t$ are defined as above.

####  A.4.2 Derivatives of $\boldsymbol{\Gamma}$ w.r.t. ${\boldsymbol{\eta}}$

To simplify the notation indicate $(\Sigma_{ll})^{-1/2}$ with
$\Sigma^{-1/2}_{ll}$. The $(l,m)$ element of $\boldsymbol{\Gamma}$ is

$\Gamma_{lm}=\Sigma^{-1/2}_{ll}\,\Sigma_{lm}\,\Sigma^{-1/2}_{mm}\,\,,$

and its partial derivative w.r.t. $\eta_{j}$, for $j=1,\dots,q$, is

$\frac{\partial\Gamma_{lm}}{\partial\eta_{j}}=\Sigma^{-1/2}_{ll}\,\frac{\partial\Sigma_{lm}}{\partial\eta_{j}}\,\Sigma^{-1/2}_{mm}-\frac{1}{2}\Sigma_{lm}\Bigg{\\{}\Sigma^{-3/2}_{ll}\,\Sigma^{-1/2}_{mm}\,\frac{\partial\Sigma_{ll}}{\partial\eta_{j}}+\Sigma^{-1/2}_{ll}\,\Sigma^{-3/2}_{mm}\,\frac{\partial\Sigma_{mm}}{\partial\eta_{j}}\Bigg{\\}}\,\,,$

and where the derivatives of the elements of ${\boldsymbol{\Sigma}}$ w.r.t.
$\eta_{j}$ are provided in Section A.4.1.